diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3e3a8c0d0..221fac574 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Change the structure of block diagonalization functions, using `BlockDiagonalForm` struct and changing the function name from `bdf` to `block_diagonal_form`. ([#349])
 - Add **GPUArrays** compatibility for `ptrace` function, by using **KernelAbstractions.jl**. ([#350])
+- Introduce `Space`, `Dimensions`, `GeneralDimensions` structures to support wider definitions and operations of `Qobj/QobjEvo`, and potential functionalities in the future. ([#271], [#353], [#360])
 
 ## [v0.24.0]
 Release date: 2024-12-13
@@ -55,6 +56,7 @@ Release date: 2024-11-13
 [v0.24.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.24.0
 [#86]: https://github.com/qutip/QuantumToolbox.jl/issues/86
 [#139]: https://github.com/qutip/QuantumToolbox.jl/issues/139
+[#271]: https://github.com/qutip/QuantumToolbox.jl/issues/271
 [#292]: https://github.com/qutip/QuantumToolbox.jl/issues/292
 [#306]: https://github.com/qutip/QuantumToolbox.jl/issues/306
 [#309]: https://github.com/qutip/QuantumToolbox.jl/issues/309
@@ -71,3 +73,5 @@ Release date: 2024-11-13
 [#347]: https://github.com/qutip/QuantumToolbox.jl/issues/347
 [#349]: https://github.com/qutip/QuantumToolbox.jl/issues/349
 [#350]: https://github.com/qutip/QuantumToolbox.jl/issues/350
+[#353]: https://github.com/qutip/QuantumToolbox.jl/issues/353
+[#360]: https://github.com/qutip/QuantumToolbox.jl/issues/360
diff --git a/docs/src/resources/api.md b/docs/src/resources/api.md
index ff890fea9..1063f21f3 100644
--- a/docs/src/resources/api.md
+++ b/docs/src/resources/api.md
@@ -11,6 +11,9 @@ CurrentModule = QuantumToolbox
 ## [Quantum object (Qobj) and type](@id doc-API:Quantum-object-and-type)
 
 ```@docs
+Space
+Dimensions
+GeneralDimensions
 AbstractQuantumObject
 BraQuantumObject
 Bra
@@ -73,7 +76,7 @@ LinearAlgebra.diag
 proj
 ptrace
 purity
-permute
+SparseArrays.permute
 tidyup
 tidyup!
 get_data
diff --git a/ext/QuantumToolboxCUDAExt.jl b/ext/QuantumToolboxCUDAExt.jl
index 7d379c3a1..459356b05 100644
--- a/ext/QuantumToolboxCUDAExt.jl
+++ b/ext/QuantumToolboxCUDAExt.jl
@@ -10,35 +10,37 @@ import SparseArrays: SparseVector, SparseMatrixCSC
 
 If `A.data` is a dense array, return a new [`QuantumObject`](@ref) where `A.data` is in the type of `CUDA.CuArray` for gpu calculations.
 """
-CuArray(A::QuantumObject{Tq}) where {Tq<:Union{Vector,Matrix}} = QuantumObject(CuArray(A.data), A.type, A.dims)
+CuArray(A::QuantumObject{Tq}) where {Tq<:Union{Vector,Matrix}} = QuantumObject(CuArray(A.data), A.type, A.dimensions)
 
 @doc raw"""
     CuArray{T}(A::QuantumObject)
 
 If `A.data` is a dense array, return a new [`QuantumObject`](@ref) where `A.data` is in the type of `CUDA.CuArray` with element type `T` for gpu calculations.
 """
-CuArray{T}(A::QuantumObject{Tq}) where {T,Tq<:Union{Vector,Matrix}} = QuantumObject(CuArray{T}(A.data), A.type, A.dims)
+CuArray{T}(A::QuantumObject{Tq}) where {T,Tq<:Union{Vector,Matrix}} =
+    QuantumObject(CuArray{T}(A.data), A.type, A.dimensions)
 
 @doc raw"""
     CuSparseVector(A::QuantumObject)
 
 If `A.data` is a sparse vector, return a new [`QuantumObject`](@ref) where `A.data` is in the type of `CUDA.CUSPARSE.CuSparseVector` for gpu calculations.
 """
-CuSparseVector(A::QuantumObject{<:SparseVector}) = QuantumObject(CuSparseVector(A.data), A.type, A.dims)
+CuSparseVector(A::QuantumObject{<:SparseVector}) = QuantumObject(CuSparseVector(A.data), A.type, A.dimensions)
 
 @doc raw"""
     CuSparseVector{T}(A::QuantumObject)
 
 If `A.data` is a sparse vector, return a new [`QuantumObject`](@ref) where `A.data` is in the type of `CUDA.CUSPARSE.CuSparseVector` with element type `T` for gpu calculations.
 """
-CuSparseVector{T}(A::QuantumObject{<:SparseVector}) where {T} = QuantumObject(CuSparseVector{T}(A.data), A.type, A.dims)
+CuSparseVector{T}(A::QuantumObject{<:SparseVector}) where {T} =
+    QuantumObject(CuSparseVector{T}(A.data), A.type, A.dimensions)
 
 @doc raw"""
     CuSparseMatrixCSC(A::QuantumObject)
 
 If `A.data` is in the type of `SparseMatrixCSC`, return a new [`QuantumObject`](@ref) where `A.data` is in the type of `CUDA.CUSPARSE.CuSparseMatrixCSC` for gpu calculations.
 """
-CuSparseMatrixCSC(A::QuantumObject{<:SparseMatrixCSC}) = QuantumObject(CuSparseMatrixCSC(A.data), A.type, A.dims)
+CuSparseMatrixCSC(A::QuantumObject{<:SparseMatrixCSC}) = QuantumObject(CuSparseMatrixCSC(A.data), A.type, A.dimensions)
 
 @doc raw"""
     CuSparseMatrixCSC{T}(A::QuantumObject)
@@ -46,14 +48,14 @@ CuSparseMatrixCSC(A::QuantumObject{<:SparseMatrixCSC}) = QuantumObject(CuSparseM
 If `A.data` is in the type of `SparseMatrixCSC`, return a new [`QuantumObject`](@ref) where `A.data` is in the type of `CUDA.CUSPARSE.CuSparseMatrixCSC` with element type `T` for gpu calculations.
 """
 CuSparseMatrixCSC{T}(A::QuantumObject{<:SparseMatrixCSC}) where {T} =
-    QuantumObject(CuSparseMatrixCSC{T}(A.data), A.type, A.dims)
+    QuantumObject(CuSparseMatrixCSC{T}(A.data), A.type, A.dimensions)
 
 @doc raw"""
     CuSparseMatrixCSR(A::QuantumObject)
 
 If `A.data` is in the type of `SparseMatrixCSC`, return a new [`QuantumObject`](@ref) where `A.data` is in the type of `CUDA.CUSPARSE.CuSparseMatrixCSR` for gpu calculations.
 """
-CuSparseMatrixCSR(A::QuantumObject{<:SparseMatrixCSC}) = QuantumObject(CuSparseMatrixCSR(A.data), A.type, A.dims)
+CuSparseMatrixCSR(A::QuantumObject{<:SparseMatrixCSC}) = QuantumObject(CuSparseMatrixCSR(A.data), A.type, A.dimensions)
 
 @doc raw"""
     CuSparseMatrixCSR(A::QuantumObject)
@@ -61,7 +63,7 @@ CuSparseMatrixCSR(A::QuantumObject{<:SparseMatrixCSC}) = QuantumObject(CuSparseM
 If `A.data` is in the type of `SparseMatrixCSC`, return a new [`QuantumObject`](@ref) where `A.data` is in the type of `CUDA.CUSPARSE.CuSparseMatrixCSR` with element type `T` for gpu calculations.
 """
 CuSparseMatrixCSR{T}(A::QuantumObject{<:SparseMatrixCSC}) where {T} =
-    QuantumObject(CuSparseMatrixCSR{T}(A.data), A.type, A.dims)
+    QuantumObject(CuSparseMatrixCSR{T}(A.data), A.type, A.dimensions)
 
 @doc raw"""
     cu(A::QuantumObject; word_size::Int=64)
diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index 28a50f251..27055e097 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -79,6 +79,8 @@ include("progress_bar.jl")
 include("linear_maps.jl")
 
 # Quantum Object
+include("qobj/space.jl")
+include("qobj/dimensions.jl")
 include("qobj/quantum_object_base.jl")
 include("qobj/quantum_object.jl")
 include("qobj/quantum_object_evo.jl")
diff --git a/src/correlations.jl b/src/correlations.jl
index 5b7068d9e..ada51d3da 100644
--- a/src/correlations.jl
+++ b/src/correlations.jl
@@ -51,8 +51,7 @@ function correlation_3op_2t(
         ψ0 = steadystate(L)
     end
 
-    allequal((L.dims, ψ0.dims, A.dims, B.dims, C.dims)) ||
-        throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
+    check_dimensions(L, ψ0, A, B, C)
 
     kwargs2 = merge((saveat = collect(tlist),), (; kwargs...))
     ρt_list = mesolve(L, ψ0, tlist; kwargs2...).states
@@ -137,7 +136,7 @@ function correlation_2op_2t(
     HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
 }
-    C = eye(prod(H.dims), dims = H.dims)
+    C = eye(prod(H.dimensions), dims = H.dimensions)
     if reverse
         corr = correlation_3op_2t(H, ψ0, tlist, τlist, c_ops, A, B, C; kwargs...)
     else
diff --git a/src/negativity.jl b/src/negativity.jl
index d91b76108..f9fa2a215 100644
--- a/src/negativity.jl
+++ b/src/negativity.jl
@@ -41,7 +41,7 @@ julia> round(negativity(ρ, 2), digits=2)
 ```
 """
 function negativity(ρ::QuantumObject, subsys::Int; logarithmic::Bool = false)
-    mask = fill(false, length(ρ.dims))
+    mask = fill(false, length(ρ.dimensions))
     try
         mask[subsys] = true
     catch
@@ -70,37 +70,46 @@ Return the partial transpose of a density matrix ``\rho``, where `mask` is an ar
 # Returns
 - `ρ_pt::QuantumObject`: The density matrix with the selected subsystems transposed.
 """
-function partial_transpose(ρ::QuantumObject{T,OperatorQuantumObject}, mask::Vector{Bool}) where {T}
-    if length(mask) != length(ρ.dims)
+function partial_transpose(ρ::QuantumObject{DT,OperatorQuantumObject}, mask::Vector{Bool}) where {DT}
+    if length(mask) != length(ρ.dimensions)
         throw(ArgumentError("The length of \`mask\` should be equal to the length of \`ρ.dims\`."))
     end
     return _partial_transpose(ρ, mask)
 end
 
 # for dense matrices
-function _partial_transpose(ρ::QuantumObject{<:AbstractArray,OperatorQuantumObject}, mask::Vector{Bool})
+function _partial_transpose(ρ::QuantumObject{DT,OperatorQuantumObject}, mask::Vector{Bool}) where {DT<:AbstractArray}
+    isa(ρ.dimensions, GeneralDimensions) &&
+        (get_dimensions_to(ρ) != get_dimensions_from(ρ)) &&
+        throw(ArgumentError("Invalid partial transpose for dims = $(_get_dims_string(ρ.dimensions))"))
+
     mask2 = [1 + Int(i) for i in mask]
     # mask2 has elements with values equal to 1 or 2
     #   1 - the subsystem don't need to be transposed
     #   2 - the subsystem need be transposed
 
     nsys = length(mask2)
+    dims = dimensions_to_dims(get_dimensions_to(ρ))
     pt_dims = reshape(Vector(1:(2*nsys)), (nsys, 2))
     pt_idx = [
-        [pt_dims[n, mask2[n]] for n in 1:nsys] # origin   value in mask2
-        [pt_dims[n, 3-mask2[n]] for n in 1:nsys]  # opposite value in mask2 (1 -> 2, and 2 -> 1)
+        [pt_dims[n, mask2[n]] for n in 1:nsys]   # origin   value in mask2
+        [pt_dims[n, 3-mask2[n]] for n in 1:nsys] # opposite value in mask2 (1 -> 2, and 2 -> 1)
     ]
     return QuantumObject(
-        reshape(permutedims(reshape(ρ.data, (ρ.dims..., ρ.dims...)), pt_idx), size(ρ)),
+        reshape(permutedims(reshape(ρ.data, (dims..., dims...)), pt_idx), size(ρ)),
         Operator,
-        ρ.dims,
+        Dimensions(ρ.dimensions.to),
     )
 end
 
 # for sparse matrices
 function _partial_transpose(ρ::QuantumObject{<:AbstractSparseArray,OperatorQuantumObject}, mask::Vector{Bool})
+    isa(ρ.dimensions, GeneralDimensions) &&
+        (get_dimensions_to(ρ) != get_dimensions_from(ρ)) &&
+        throw(ArgumentError("Invalid partial transpose for dims = $(_get_dims_string(ρ.dimensions))"))
+
     M, N = size(ρ)
-    dimsTuple = Tuple(ρ.dims)
+    dimsTuple = Tuple(dimensions_to_dims(get_dimensions_to(ρ)))
     colptr = ρ.data.colptr
     rowval = ρ.data.rowval
     nzval = ρ.data.nzval
@@ -132,5 +141,5 @@ function _partial_transpose(ρ::QuantumObject{<:AbstractSparseArray,OperatorQuan
         end
     end
 
-    return QuantumObject(sparse(I_pt, J_pt, V_pt, M, N), Operator, ρ.dims)
+    return QuantumObject(sparse(I_pt, J_pt, V_pt, M, N), Operator, ρ.dimensions)
 end
diff --git a/src/qobj/arithmetic_and_attributes.jl b/src/qobj/arithmetic_and_attributes.jl
index f9af5aeeb..0ae908b6f 100644
--- a/src/qobj/arithmetic_and_attributes.jl
+++ b/src/qobj/arithmetic_and_attributes.jl
@@ -4,7 +4,7 @@ Arithmetic and Attributes for QuantumObject
     - export most of the attribute functions in "Python Qobj class"
 =#
 
-export proj, ptrace, purity, permute
+export proj, ptrace, purity
 export tidyup, tidyup!
 export get_data, get_coherence
 
@@ -13,23 +13,23 @@ Base.broadcastable(x::QuantumObject) = x.data
 for op in (:(+), :(-), :(*), :(/), :(^))
     @eval begin
         function Base.Broadcast.broadcasted(::typeof($op), x::QuantumObject, y::QuantumObject)
-            return QuantumObject(broadcast($op, x.data, y.data), x.type, x.dims)
+            return QuantumObject(broadcast($op, x.data, y.data), x.type, x.dimensions)
         end
 
         function Base.Broadcast.broadcasted(::typeof($op), x::QuantumObject, y::Number)
-            return QuantumObject(broadcast($op, x.data, y), x.type, x.dims)
+            return QuantumObject(broadcast($op, x.data, y), x.type, x.dimensions)
         end
 
         function Base.Broadcast.broadcasted(::typeof($op), x::Number, y::QuantumObject)
-            return QuantumObject(broadcast($op, x, y.data), y.type, y.dims)
+            return QuantumObject(broadcast($op, x, y.data), y.type, y.dimensions)
         end
 
         function Base.Broadcast.broadcasted(::typeof($op), x::QuantumObject, y::AbstractArray)
-            return QuantumObject(broadcast($op, x.data, y), x.type, x.dims)
+            return QuantumObject(broadcast($op, x.data, y), x.type, x.dimensions)
         end
 
         function Base.Broadcast.broadcasted(::typeof($op), x::AbstractArray, y::QuantumObject)
-            return QuantumObject(broadcast($op, x, y.data), y.type, y.dims)
+            return QuantumObject(broadcast($op, x, y.data), y.type, y.dimensions)
         end
     end
 end
@@ -37,79 +37,120 @@ end
 for op in (:(+), :(-), :(*))
     @eval begin
         function LinearAlgebra.$op(A::AbstractQuantumObject, B::AbstractQuantumObject)
-            check_dims(A, B)
+            check_dimensions(A, B)
             QType = promote_op_type(A, B)
-            return QType($(op)(A.data, B.data), A.type, A.dims)
+            return QType($(op)(A.data, B.data), A.type, A.dimensions)
         end
-        LinearAlgebra.$op(A::AbstractQuantumObject) = get_typename_wrapper(A)($(op)(A.data), A.type, A.dims)
+        LinearAlgebra.$op(A::AbstractQuantumObject) = get_typename_wrapper(A)($(op)(A.data), A.type, A.dimensions)
 
         LinearAlgebra.$op(n::T, A::AbstractQuantumObject) where {T<:Number} =
-            get_typename_wrapper(A)($(op)(n * I, A.data), A.type, A.dims)
+            get_typename_wrapper(A)($(op)(n * I, A.data), A.type, A.dimensions)
         LinearAlgebra.$op(A::AbstractQuantumObject, n::T) where {T<:Number} =
-            get_typename_wrapper(A)($(op)(A.data, n * I), A.type, A.dims)
+            get_typename_wrapper(A)($(op)(A.data, n * I), A.type, A.dimensions)
+    end
+end
+
+function check_mul_dimensions(from::NTuple{NA,AbstractSpace}, to::NTuple{NB,AbstractSpace}) where {NA,NB}
+    (from != to) && throw(
+        DimensionMismatch(
+            "The quantum object with (right) dims = $(dimensions_to_dims(from)) can not multiply a quantum object with (left) dims = $(dimensions_to_dims(to)) on the right-hand side.",
+        ),
+    )
+    return nothing
+end
+
+for ADimType in (:Dimensions, :GeneralDimensions)
+    for BDimType in (:Dimensions, :GeneralDimensions)
+        if ADimType == BDimType == :Dimensions
+            @eval begin
+                function LinearAlgebra.:(*)(
+                    A::AbstractQuantumObject{DT1,OperatorQuantumObject,<:$ADimType},
+                    B::AbstractQuantumObject{DT2,OperatorQuantumObject,<:$BDimType},
+                ) where {DT1,DT2}
+                    check_dimensions(A, B)
+                    QType = promote_op_type(A, B)
+                    return QType(A.data * B.data, Operator, A.dimensions)
+                end
+            end
+        else
+            @eval begin
+                function LinearAlgebra.:(*)(
+                    A::AbstractQuantumObject{DT1,OperatorQuantumObject,<:$ADimType},
+                    B::AbstractQuantumObject{DT2,OperatorQuantumObject,<:$BDimType},
+                ) where {DT1,DT2}
+                    check_mul_dimensions(get_dimensions_from(A), get_dimensions_to(B))
+                    QType = promote_op_type(A, B)
+                    return QType(
+                        A.data * B.data,
+                        Operator,
+                        GeneralDimensions(get_dimensions_to(A), get_dimensions_from(B)),
+                    )
+                end
+            end
+        end
     end
 end
 
 function LinearAlgebra.:(*)(
     A::AbstractQuantumObject{DT1,OperatorQuantumObject},
-    B::QuantumObject{DT2,KetQuantumObject},
+    B::QuantumObject{DT2,KetQuantumObject,<:Dimensions},
 ) where {DT1,DT2}
-    check_dims(A, B)
-    return QuantumObject(A.data * B.data, Ket, A.dims)
+    check_mul_dimensions(get_dimensions_from(A), get_dimensions_to(B))
+    return QuantumObject(A.data * B.data, Ket, Dimensions(get_dimensions_to(A)))
 end
 function LinearAlgebra.:(*)(
-    A::QuantumObject{DT1,BraQuantumObject},
+    A::QuantumObject{DT1,BraQuantumObject,<:Dimensions},
     B::AbstractQuantumObject{DT2,OperatorQuantumObject},
 ) where {DT1,DT2}
-    check_dims(A, B)
-    return QuantumObject(A.data * B.data, Bra, A.dims)
+    check_mul_dimensions(get_dimensions_from(A), get_dimensions_to(B))
+    return QuantumObject(A.data * B.data, Bra, Dimensions(get_dimensions_from(B)))
 end
 function LinearAlgebra.:(*)(
     A::QuantumObject{DT1,KetQuantumObject},
     B::QuantumObject{DT2,BraQuantumObject},
 ) where {DT1,DT2}
-    check_dims(A, B)
-    return QuantumObject(A.data * B.data, Operator, A.dims)
+    check_dimensions(A, B)
+    return QuantumObject(A.data * B.data, Operator, A.dimensions) # to align with QuTiP, don't use kron(A, B) to do it.
 end
 function LinearAlgebra.:(*)(
     A::QuantumObject{DT1,BraQuantumObject},
     B::QuantumObject{DT2,KetQuantumObject},
 ) where {DT1,DT2}
-    check_dims(A, B)
+    check_dimensions(A, B)
     return A.data * B.data
 end
 function LinearAlgebra.:(*)(
     A::AbstractQuantumObject{DT1,SuperOperatorQuantumObject},
     B::QuantumObject{DT2,OperatorQuantumObject},
 ) where {DT1,DT2}
-    check_dims(A, B)
-    return QuantumObject(vec2mat(A.data * mat2vec(B.data)), Operator, A.dims)
+    check_dimensions(A, B)
+    return QuantumObject(vec2mat(A.data * mat2vec(B.data)), Operator, A.dimensions)
 end
 function LinearAlgebra.:(*)(
     A::QuantumObject{DT1,OperatorBraQuantumObject},
     B::QuantumObject{DT2,OperatorKetQuantumObject},
 ) where {DT1,DT2}
-    check_dims(A, B)
+    check_dimensions(A, B)
     return A.data * B.data
 end
 function LinearAlgebra.:(*)(
     A::AbstractQuantumObject{DT1,SuperOperatorQuantumObject},
     B::QuantumObject{DT2,OperatorKetQuantumObject},
 ) where {DT1,DT2}
-    check_dims(A, B)
-    return QuantumObject(A.data * B.data, OperatorKet, A.dims)
+    check_dimensions(A, B)
+    return QuantumObject(A.data * B.data, OperatorKet, A.dimensions)
 end
 function LinearAlgebra.:(*)(
     A::QuantumObject{<:AbstractArray{T1},OperatorBraQuantumObject},
     B::AbstractQuantumObject{<:AbstractArray{T2},SuperOperatorQuantumObject},
 ) where {T1,T2}
-    check_dims(A, B)
-    return QuantumObject(A.data * B.data, OperatorBra, A.dims)
+    check_dimensions(A, B)
+    return QuantumObject(A.data * B.data, OperatorBra, A.dimensions)
 end
 
-LinearAlgebra.:(^)(A::QuantumObject{DT}, n::T) where {DT,T<:Number} = QuantumObject(^(A.data, n), A.type, A.dims)
+LinearAlgebra.:(^)(A::QuantumObject{DT}, n::T) where {DT,T<:Number} = QuantumObject(^(A.data, n), A.type, A.dimensions)
 LinearAlgebra.:(/)(A::AbstractQuantumObject{DT}, n::T) where {DT,T<:Number} =
-    get_typename_wrapper(A)(A.data / n, A.type, A.dims)
+    get_typename_wrapper(A)(A.data / n, A.type, A.dimensions)
 
 @doc raw"""
     A ⋅ B
@@ -126,7 +167,7 @@ function LinearAlgebra.dot(
     A::QuantumObject{DT1,OpType},
     B::QuantumObject{DT2,OpType},
 ) where {DT1,DT2,OpType<:Union{KetQuantumObject,OperatorKetQuantumObject}}
-    A.dims != B.dims && throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
+    check_dimensions(A, B)
     return LinearAlgebra.dot(A.data, B.data)
 end
 
@@ -148,8 +189,7 @@ function LinearAlgebra.dot(
     A::AbstractQuantumObject{DT2,OperatorQuantumObject},
     j::QuantumObject{DT3,KetQuantumObject},
 ) where {DT1,DT2,DT3}
-    ((i.dims != A.dims) || (A.dims != j.dims)) &&
-        throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
+    check_dimensions(i, A, j)
     return LinearAlgebra.dot(i.data, A.data, j.data)
 end
 function LinearAlgebra.dot(
@@ -157,8 +197,7 @@ function LinearAlgebra.dot(
     A::AbstractQuantumObject{DT2,SuperOperatorQuantumObject},
     j::QuantumObject{DT3,OperatorKetQuantumObject},
 ) where {DT1,DT2,DT3}
-    ((i.dims != A.dims) || (A.dims != j.dims)) &&
-        throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
+    check_dimensions(i, A, j)
     return LinearAlgebra.dot(i.data, A.data, j.data)
 end
 
@@ -167,7 +206,7 @@ end
 
 Return a similar [`AbstractQuantumObject`](@ref) with `dims` and `type` are same as `A`, but `data` is a zero-array.
 """
-Base.zero(A::AbstractQuantumObject) = get_typename_wrapper(A)(zero(A.data), A.type, A.dims)
+Base.zero(A::AbstractQuantumObject) = get_typename_wrapper(A)(zero(A.data), A.type, A.dimensions)
 
 @doc raw"""
     one(A::AbstractQuantumObject)
@@ -179,14 +218,14 @@ Note that `A` must be [`Operator`](@ref) or [`SuperOperator`](@ref).
 Base.one(
     A::AbstractQuantumObject{DT,OpType},
 ) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    get_typename_wrapper(A)(one(A.data), A.type, A.dims)
+    get_typename_wrapper(A)(one(A.data), A.type, A.dimensions)
 
 @doc raw"""
     conj(A::AbstractQuantumObject)
 
 Return the element-wise complex conjugation of the [`AbstractQuantumObject`](@ref).
 """
-Base.conj(A::AbstractQuantumObject) = get_typename_wrapper(A)(conj(A.data), A.type, A.dims)
+Base.conj(A::AbstractQuantumObject) = get_typename_wrapper(A)(conj(A.data), A.type, A.dimensions)
 
 @doc raw"""
     transpose(A::AbstractQuantumObject)
@@ -196,7 +235,7 @@ Lazy matrix transpose of the [`AbstractQuantumObject`](@ref).
 LinearAlgebra.transpose(
     A::AbstractQuantumObject{DT,OpType},
 ) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    get_typename_wrapper(A)(transpose(A.data), A.type, A.dims)
+    get_typename_wrapper(A)(transpose(A.data), A.type, transpose(A.dimensions))
 
 @doc raw"""
     A'
@@ -211,13 +250,15 @@ Lazy adjoint (conjugate transposition) of the [`AbstractQuantumObject`](@ref)
 LinearAlgebra.adjoint(
     A::AbstractQuantumObject{DT,OpType},
 ) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    get_typename_wrapper(A)(adjoint(A.data), A.type, A.dims)
-LinearAlgebra.adjoint(A::QuantumObject{DT,KetQuantumObject}) where {DT} = QuantumObject(adjoint(A.data), Bra, A.dims)
-LinearAlgebra.adjoint(A::QuantumObject{DT,BraQuantumObject}) where {DT} = QuantumObject(adjoint(A.data), Ket, A.dims)
+    get_typename_wrapper(A)(adjoint(A.data), A.type, adjoint(A.dimensions))
+LinearAlgebra.adjoint(A::QuantumObject{DT,KetQuantumObject}) where {DT} =
+    QuantumObject(adjoint(A.data), Bra, adjoint(A.dimensions))
+LinearAlgebra.adjoint(A::QuantumObject{DT,BraQuantumObject}) where {DT} =
+    QuantumObject(adjoint(A.data), Ket, adjoint(A.dimensions))
 LinearAlgebra.adjoint(A::QuantumObject{DT,OperatorKetQuantumObject}) where {DT} =
-    QuantumObject(adjoint(A.data), OperatorBra, A.dims)
+    QuantumObject(adjoint(A.data), OperatorBra, adjoint(A.dimensions))
 LinearAlgebra.adjoint(A::QuantumObject{DT,OperatorBraQuantumObject}) where {DT} =
-    QuantumObject(adjoint(A.data), OperatorKet, A.dims)
+    QuantumObject(adjoint(A.data), OperatorKet, adjoint(A.dimensions))
 
 @doc raw"""
     inv(A::AbstractQuantumObject)
@@ -227,13 +268,13 @@ Matrix inverse of the [`AbstractQuantumObject`](@ref). If `A` is a [`QuantumObje
 LinearAlgebra.inv(
     A::AbstractQuantumObject{DT,OpType},
 ) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    QuantumObject(sparse(inv(Matrix(A.data))), A.type, A.dims)
+    QuantumObject(sparse(inv(Matrix(A.data))), A.type, A.dimensions)
 
 LinearAlgebra.Hermitian(
     A::QuantumObject{DT,OpType},
     uplo::Symbol = :U,
 ) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    QuantumObject(Hermitian(A.data, uplo), A.type, A.dims)
+    QuantumObject(Hermitian(A.data, uplo), A.type, A.dimensions)
 
 @doc raw"""
     tr(A::QuantumObject)
@@ -336,9 +377,9 @@ Also, see [`norm`](@ref) about its definition for different types of [`QuantumOb
 LinearAlgebra.normalize(
     A::QuantumObject{<:AbstractArray{T},ObjType},
     p::Real = 2,
-) where {T,ObjType<:Union{KetQuantumObject,BraQuantumObject}} = QuantumObject(A.data / norm(A, p), A.type, A.dims)
+) where {T,ObjType<:Union{KetQuantumObject,BraQuantumObject}} = QuantumObject(A.data / norm(A, p), A.type, A.dimensions)
 LinearAlgebra.normalize(A::QuantumObject{<:AbstractArray{T},OperatorQuantumObject}, p::Real = 1) where {T} =
-    QuantumObject(A.data / norm(A, p), A.type, A.dims)
+    QuantumObject(A.data / norm(A, p), A.type, A.dimensions)
 
 @doc raw"""
     normalize!(A::QuantumObject, p::Real)
@@ -370,12 +411,12 @@ LinearAlgebra.triu(
     A::QuantumObject{<:AbstractArray{T},OpType},
     k::Integer = 0,
 ) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    QuantumObject(triu(A.data, k), A.type, A.dims)
+    QuantumObject(triu(A.data, k), A.type, A.dimensions)
 LinearAlgebra.tril(
     A::QuantumObject{<:AbstractArray{T},OpType},
     k::Integer = 0,
 ) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    QuantumObject(tril(A.data, k), A.type, A.dims)
+    QuantumObject(tril(A.data, k), A.type, A.dimensions)
 
 LinearAlgebra.lmul!(a::Number, B::QuantumObject{<:AbstractArray}) = (lmul!(a, B.data); B)
 LinearAlgebra.rmul!(B::QuantumObject{<:AbstractArray}, a::Number) = (rmul!(B.data, a); B)
@@ -393,7 +434,7 @@ Matrix square root of [`QuantumObject`](@ref)
     `√(A)` (where `√` can be typed by tab-completing `\sqrt` in the REPL) is a synonym of `sqrt(A)`.
 """
 LinearAlgebra.sqrt(A::QuantumObject{<:AbstractArray{T}}) where {T} =
-    QuantumObject(sqrt(sparse_to_dense(A.data)), A.type, A.dims)
+    QuantumObject(sqrt(sparse_to_dense(A.data)), A.type, A.dimensions)
 
 @doc raw"""
     log(A::QuantumObject)
@@ -405,7 +446,7 @@ Note that this function only supports for [`Operator`](@ref) and [`SuperOperator
 LinearAlgebra.log(
     A::QuantumObject{<:AbstractMatrix{T},ObjType},
 ) where {T,ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    QuantumObject(log(sparse_to_dense(A.data)), A.type, A.dims)
+    QuantumObject(log(sparse_to_dense(A.data)), A.type, A.dimensions)
 
 @doc raw"""
     exp(A::QuantumObject)
@@ -417,11 +458,11 @@ Note that this function only supports for [`Operator`](@ref) and [`SuperOperator
 LinearAlgebra.exp(
     A::QuantumObject{<:AbstractMatrix{T},ObjType},
 ) where {T,ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    QuantumObject(dense_to_sparse(exp(A.data)), A.type, A.dims)
+    QuantumObject(dense_to_sparse(exp(A.data)), A.type, A.dimensions)
 LinearAlgebra.exp(
     A::QuantumObject{<:AbstractSparseMatrix{T},ObjType},
 ) where {T,ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    QuantumObject(_spexp(A.data), A.type, A.dims)
+    QuantumObject(_spexp(A.data), A.type, A.dimensions)
 
 function _spexp(A::SparseMatrixCSC{T,M}; threshold = 1e-14, nonzero_tol = 1e-20) where {T,M}
     m = checksquare(A) # Throws exception if not square
@@ -554,7 +595,7 @@ function ptrace(QO::QuantumObject{<:AbstractArray,KetQuantumObject}, sel::Union{
     if n_s == 0 # return full trace for empty sel
         return tr(ket2dm(QO))
     else
-        n_d = length(QO.dims)
+        n_d = length(QO.dimensions)
 
         (any(>(n_d), sel) || any(<(1), sel)) && throw(
             ArgumentError("Invalid indices in `sel`: $(sel), the given QuantumObject only have $(n_d) sub-systems"),
@@ -565,19 +606,24 @@ function ptrace(QO::QuantumObject{<:AbstractArray,KetQuantumObject}, sel::Union{
 
     _sort_sel = sort(SVector{length(sel),Int}(sel))
     ρtr, dkeep = _ptrace_ket(QO.data, QO.dims, _sort_sel)
-    return QuantumObject(ρtr, type = Operator, dims = dkeep)
+    return QuantumObject(ρtr, type = Operator, dims = Dimensions(dkeep))
 end
 
 ptrace(QO::QuantumObject{<:AbstractArray,BraQuantumObject}, sel::Union{AbstractVector{Int},Tuple}) = ptrace(QO', sel)
 
 function ptrace(QO::QuantumObject{<:AbstractArray,OperatorQuantumObject}, sel::Union{AbstractVector{Int},Tuple})
+    # TODO: support for special cases when some of the subsystems have same `to` and `from` space
+    isa(QO.dimensions, GeneralDimensions) &&
+        (get_dimensions_to(QO) != get_dimensions_from(QO)) &&
+        throw(ArgumentError("Invalid partial trace for dims = $(_get_dims_string(QO.dimensions))"))
+
     _non_static_array_warning("sel", sel)
 
     n_s = length(sel)
     if n_s == 0 # return full trace for empty sel
         return tr(QO)
     else
-        n_d = length(QO.dims)
+        n_d = length(QO.dimensions)
 
         (any(>(n_d), sel) || any(<(1), sel)) && throw(
             ArgumentError("Invalid indices in `sel`: $(sel), the given QuantumObject only have $(n_d) sub-systems"),
@@ -586,9 +632,10 @@ function ptrace(QO::QuantumObject{<:AbstractArray,OperatorQuantumObject}, sel::U
         (n_d == 1) && return QO
     end
 
+    dims = dimensions_to_dims(get_dimensions_to(QO))
     _sort_sel = sort(SVector{length(sel),Int}(sel))
-    ρtr, dkeep = _ptrace_oper(QO.data, QO.dims, _sort_sel)
-    return QuantumObject(ρtr, type = Operator, dims = dkeep)
+    ρtr, dkeep = _ptrace_oper(QO.data, dims, _sort_sel)
+    return QuantumObject(ρtr, type = Operator, dims = Dimensions(dkeep))
 end
 ptrace(QO::QuantumObject, sel::Int) = ptrace(QO, SVector(sel))
 
@@ -674,7 +721,7 @@ purity(ρ::QuantumObject{<:AbstractArray{T},OperatorQuantumObject}) where {T} =
 Given a [`QuantumObject`](@ref) `A`, check the real and imaginary parts of each element separately. Remove the real or imaginary value if its absolute value is less than `tol`.
 """
 tidyup(A::QuantumObject{<:AbstractArray{T}}, tol::T2 = 1e-14) where {T,T2<:Real} =
-    QuantumObject(tidyup(A.data, tol), A.type, A.dims)
+    QuantumObject(tidyup(A.data, tol), A.type, A.dimensions)
 tidyup(A::AbstractArray{T}, tol::T2 = 1e-14) where {T,T2<:Real} = tidyup!(copy(A), tol)
 
 @doc raw"""
@@ -698,7 +745,7 @@ tidyup!(A::AbstractArray{T}, tol::T2 = 1e-14) where {T,T2<:Real} =
 
 Returns the data of a [`AbstractQuantumObject`](@ref).
 """
-get_data(A::AbstractQuantumObject) = A.data
+get_data(A::AbstractQuantumObject) = getfield(A, :data)
 
 @doc raw"""
     get_coherence(ψ::QuantumObject)
@@ -708,7 +755,7 @@ Get the coherence value ``\alpha`` by measuring the expectation value of the des
 It returns both ``\alpha`` and the corresponding state with the coherence removed: ``\ket{\delta_\alpha} = \exp ( \alpha^* \hat{a} - \alpha \hat{a}^\dagger ) \ket{\psi}`` for a pure state, and ``\hat{\rho_\alpha} = \exp ( \alpha^* \hat{a} - \alpha \hat{a}^\dagger ) \hat{\rho} \exp ( -\bar{\alpha} \hat{a} + \alpha \hat{a}^\dagger )`` for a density matrix. These states correspond to the quantum fluctuations around the coherent state ``\ket{\alpha}`` or ``|\alpha\rangle\langle\alpha|``.
 """
 function get_coherence(ψ::QuantumObject{<:AbstractArray,KetQuantumObject})
-    a = destroy(prod(ψ.dims))
+    a = destroy(prod(ψ.dimensions))
     α = expect(a, ψ)
     D = exp(α * a' - conj(α) * a)
 
@@ -716,7 +763,7 @@ function get_coherence(ψ::QuantumObject{<:AbstractArray,KetQuantumObject})
 end
 
 function get_coherence(ρ::QuantumObject{<:AbstractArray,OperatorQuantumObject})
-    a = destroy(prod(ρ.dims))
+    a = destroy(prod(ρ.dimensions))
     α = expect(a, ρ)
     D = exp(α * a' - conj(α) * a)
 
@@ -750,11 +797,11 @@ true
 !!! warning "Beware of type-stability!"
     It is highly recommended to use `permute(A, order)` with `order` as `Tuple` or `SVector` to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
-function permute(
+function SparseArrays.permute(
     A::QuantumObject{<:AbstractArray{T},ObjType},
     order::Union{AbstractVector{Int},Tuple},
 ) where {T,ObjType<:Union{KetQuantumObject,BraQuantumObject,OperatorQuantumObject}}
-    (length(order) != length(A.dims)) &&
+    (length(order) != length(A.dimensions)) &&
         throw(ArgumentError("The order list must have the same length as the number of subsystems (A.dims)"))
 
     !isperm(order) && throw(ArgumentError("$(order) is not a valid permutation of the subsystems (A.dims)"))
@@ -766,22 +813,26 @@ function permute(
     # obtain the arguments: dims for reshape; perm for PermutedDimsArray
     dims, perm = _dims_and_perm(A.type, A.dims, order_svector, length(order_svector))
 
-    return QuantumObject(
-        reshape(permutedims(reshape(A.data, dims...), Tuple(perm)), size(A)),
-        A.type,
-        A.dims[order_svector],
-    )
+    order_dimensions = _order_dimensions(A.dimensions, order_svector)
+
+    return QuantumObject(reshape(permutedims(reshape(A.data, dims...), Tuple(perm)), size(A)), A.type, order_dimensions)
 end
 
-function _dims_and_perm(
+_dims_and_perm(
     ::ObjType,
     dims::SVector{N,Int},
     order::AbstractVector{Int},
     L::Int,
-) where {ObjType<:Union{KetQuantumObject,BraQuantumObject},N}
-    return reverse(dims), reverse((L + 1) .- order)
-end
+) where {ObjType<:Union{KetQuantumObject,BraQuantumObject},N} = reverse(dims), reverse((L + 1) .- order)
 
-function _dims_and_perm(::OperatorQuantumObject, dims::SVector{N,Int}, order::AbstractVector{Int}, L::Int) where {N}
-    return reverse(vcat(dims, dims)), reverse((2 * L + 1) .- vcat(order, order .+ L))
-end
+# if dims originates from Dimensions
+_dims_and_perm(::OperatorQuantumObject, dims::SVector{N,Int}, order::AbstractVector{Int}, L::Int) where {N} =
+    reverse(vcat(dims, dims)), reverse((2 * L + 1) .- vcat(order, order .+ L))
+
+# if dims originates from GeneralDimensions
+_dims_and_perm(::OperatorQuantumObject, dims::SVector{2,SVector{N,Int}}, order::AbstractVector{Int}, L::Int) where {N} =
+    reverse(vcat(dims[2], dims[1])), reverse((2 * L + 1) .- vcat(order, order .+ L))
+
+_order_dimensions(dimensions::Dimensions, order::AbstractVector{Int}) = Dimensions(dimensions.to[order])
+_order_dimensions(dimensions::GeneralDimensions, order::AbstractVector{Int}) =
+    GeneralDimensions(dimensions.to[order], dimensions.from[order])
diff --git a/src/qobj/block_diagonal_form.jl b/src/qobj/block_diagonal_form.jl
index cc1f97eaa..ca3e89a68 100644
--- a/src/qobj/block_diagonal_form.jl
+++ b/src/qobj/block_diagonal_form.jl
@@ -54,7 +54,7 @@ function block_diagonal_form(
     A::QuantumObject{DT,OpType},
 ) where {DT<:AbstractSparseMatrix,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
     bdf = block_diagonal_form(A.data)
-    B = QuantumObject(bdf.B, type = A.type, dims = A.dims)
-    P = QuantumObject(bdf.P, type = A.type, dims = A.dims)
+    B = QuantumObject(bdf.B, type = A.type, dims = A.dimensions)
+    P = QuantumObject(bdf.P, type = A.type, dims = A.dimensions)
     return BlockDiagonalForm(B, P, bdf.blocks, bdf.block_sizes)
 end
diff --git a/src/qobj/dimensions.jl b/src/qobj/dimensions.jl
new file mode 100644
index 000000000..0fd19f3a9
--- /dev/null
+++ b/src/qobj/dimensions.jl
@@ -0,0 +1,104 @@
+#=
+This file defines the Dimensions structures, which can describe composite Hilbert spaces.
+=#
+
+export AbstractDimensions, Dimensions, GeneralDimensions
+
+abstract type AbstractDimensions{N} end
+
+@doc raw"""
+    struct Dimensions{N,T<:Tuple} <: AbstractDimensions{N}
+        to::T
+    end
+
+A structure that describes the Hilbert [`Space`](@ref) of each subsystems.
+"""
+struct Dimensions{N,T<:Tuple} <: AbstractDimensions{N}
+    to::T
+
+    # make sure the elements in the tuple are all AbstractSpace
+    Dimensions(to::NTuple{N,T}) where {N,T<:AbstractSpace} = new{N,typeof(to)}(to)
+end
+function Dimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Integer,N}
+    _non_static_array_warning("dims", dims)
+    L = length(dims)
+    (L > 0) || throw(DomainError(dims, "The argument dims must be of non-zero length"))
+
+    return Dimensions(NTuple{L,Space}(Space.(dims)))
+end
+Dimensions(dims::Int) = Dimensions(Space(dims))
+Dimensions(dims::DimType) where {DimType<:AbstractSpace} = Dimensions((dims,))
+Dimensions(dims::Any) = throw(
+    ArgumentError(
+        "The argument dims must be a Tuple or a StaticVector of non-zero length and contain only positive integers.",
+    ),
+)
+
+@doc raw"""
+    struct GeneralDimensions{N,T1<:Tuple,T2<:Tuple} <: AbstractDimensions{N}
+        to::T1
+        from::T2
+    end
+
+A structure that describes the left-hand side (`to`) and right-hand side (`from`) Hilbert [`Space`](@ref) of an [`Operator`](@ref).
+"""
+struct GeneralDimensions{N,T1<:Tuple,T2<:Tuple} <: AbstractDimensions{N}
+    # note that the number `N` should be the same for both `to` and `from`
+    to::T1   # space acting on the left
+    from::T2 # space acting on the right
+
+    # make sure the elements in the tuple are all AbstractSpace
+    GeneralDimensions(to::NTuple{N,T1}, from::NTuple{N,T2}) where {N,T1<:AbstractSpace,T2<:AbstractSpace} =
+        new{N,typeof(to),typeof(from)}(to, from)
+end
+function GeneralDimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Union{AbstractVector,NTuple},N}
+    (length(dims) != 2) && throw(ArgumentError("Invalid dims = $dims"))
+
+    _non_static_array_warning("dims[1]", dims[1])
+    _non_static_array_warning("dims[2]", dims[2])
+
+    L1 = length(dims[1])
+    L2 = length(dims[2])
+    ((L1 > 0) && (L1 == L2)) || throw(
+        DomainError(
+            (L1, L2),
+            "The length of the arguments `dims[1]` and `dims[2]` must be in the same length and have at least one element.",
+        ),
+    )
+
+    return GeneralDimensions(NTuple{L1,Space}(Space.(dims[1])), NTuple{L1,Space}(Space.(dims[2])))
+end
+
+_gen_dimensions(dims::AbstractDimensions) = dims
+_gen_dimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Integer,N} = Dimensions(dims)
+_gen_dimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Union{AbstractVector,NTuple},N} =
+    GeneralDimensions(dims)
+_gen_dimensions(dims::Any) = Dimensions(dims)
+
+# obtain dims in the type of SVector with integers
+dimensions_to_dims(dimensions::NTuple{N,AbstractSpace}) where {N} = vcat(map(dimensions_to_dims, dimensions)...)
+dimensions_to_dims(dimensions::Dimensions) = dimensions_to_dims(dimensions.to)
+dimensions_to_dims(dimensions::GeneralDimensions) =
+    SVector{2}(dimensions_to_dims(dimensions.to), dimensions_to_dims(dimensions.from))
+
+dimensions_to_dims(::Nothing) = nothing # for EigsolveResult.dimensions = nothing
+
+Base.length(::AbstractDimensions{N}) where {N} = N
+
+# need to specify return type `Int` for `_get_space_size`
+# otherwise the type of `prod(::Dimensions)` will be unstable
+_get_space_size(s::AbstractSpace)::Int = s.size
+Base.prod(dims::Dimensions) = prod(dims.to)
+Base.prod(spaces::NTuple{N,AbstractSpace}) where {N} = prod(_get_space_size, spaces)
+
+LinearAlgebra.transpose(dimensions::Dimensions) = dimensions
+LinearAlgebra.transpose(dimensions::GeneralDimensions) = GeneralDimensions(dimensions.from, dimensions.to) # switch `to` and `from`
+LinearAlgebra.adjoint(dimensions::AbstractDimensions) = transpose(dimensions)
+
+# this is used to show `dims` for Qobj and QobjEvo
+_get_dims_string(dimensions::Dimensions) = string(dimensions_to_dims(dimensions))
+function _get_dims_string(dimensions::GeneralDimensions)
+    dims = dimensions_to_dims(dimensions)
+    return "[$(string(dims[1])), $(string(dims[2]))]"
+end
+_get_dims_string(::Nothing) = "nothing" # for EigsolveResult.dimensions = nothing
diff --git a/src/qobj/eigsolve.jl b/src/qobj/eigsolve.jl
index 3728ffc5d..cbccbed2d 100644
--- a/src/qobj/eigsolve.jl
+++ b/src/qobj/eigsolve.jl
@@ -7,27 +7,22 @@ export eigenenergies, eigenstates, eigsolve
 export eigsolve_al
 
 @doc raw"""
-    struct EigsolveResult{T1<:Vector{<:Number}, T2<:AbstractMatrix{<:Number}, ObjType<:Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject},N}
-        values::T1
-        vectors::T2
-        type::ObjType
-        dims::SVector{N,Int}
-        iter::Int
-        numops::Int
-        converged::Bool
-    end
+    struct EigsolveResult
 
 A struct containing the eigenvalues, the eigenvectors, and some information from the solver
 
-# Fields
+# Fields (Attributes)
 - `values::AbstractVector`: the eigenvalues
 - `vectors::AbstractMatrix`: the transformation matrix (eigenvectors)
 - `type::Union{Nothing,QuantumObjectType}`: the type of [`QuantumObject`](@ref), or `nothing` means solving eigen equation for general matrix
-- `dims::SVector`: the `dims` of [`QuantumObject`](@ref)
+- `dimensions::Union{Nothing,AbstractDimensions}`: the `dimensions` of [`QuantumObject`](@ref), or `nothing` means solving eigen equation for general matrix
 - `iter::Int`: the number of iteration during the solving process
 - `numops::Int` : number of times the linear map was applied in krylov methods
 - `converged::Bool`: Whether the result is converged
 
+!!! note "`dims` property"
+    For a given `res::EigsolveResult`, `res.dims` or `getproperty(res, :dims)` returns its `dimensions` in the type of integer-vector.
+
 # Examples
 One can obtain the eigenvalues and the corresponding [`QuantumObject`](@ref)-type eigenvectors by:
 ```jldoctest
@@ -50,7 +45,7 @@ julia> λ
   1.0 + 0.0im
 
 julia> ψ
-2-element Vector{QuantumObject{Vector{ComplexF64}, KetQuantumObject, 1}}:
+2-element Vector{QuantumObject{Vector{ComplexF64}, KetQuantumObject, Dimensions{1, Tuple{Space}}}}:
 
 Quantum Object:   type=Ket   dims=[2]   size=(2,)
 2-element Vector{ComplexF64}:
@@ -72,29 +67,38 @@ struct EigsolveResult{
     T1<:Vector{<:Number},
     T2<:AbstractMatrix{<:Number},
     ObjType<:Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject},
-    N,
+    DimType<:Union{Nothing,AbstractDimensions},
 }
     values::T1
     vectors::T2
     type::ObjType
-    dims::SVector{N,Int}
+    dimensions::DimType
     iter::Int
     numops::Int
     converged::Bool
 end
 
+function Base.getproperty(res::EigsolveResult, key::Symbol)
+    # a comment here to avoid bad render by JuliaFormatter
+    if key === :dims
+        return dimensions_to_dims(getfield(res, :dimensions))
+    else
+        return getfield(res, key)
+    end
+end
+
 Base.iterate(res::EigsolveResult) = (res.values, Val(:vector_list))
 Base.iterate(res::EigsolveResult{T1,T2,Nothing}, ::Val{:vector_list}) where {T1,T2} =
     ([res.vectors[:, k] for k in 1:length(res.values)], Val(:vectors))
 Base.iterate(res::EigsolveResult{T1,T2,OperatorQuantumObject}, ::Val{:vector_list}) where {T1,T2} =
-    ([QuantumObject(res.vectors[:, k], Ket, res.dims) for k in 1:length(res.values)], Val(:vectors))
+    ([QuantumObject(res.vectors[:, k], Ket, res.dimensions) for k in 1:length(res.values)], Val(:vectors))
 Base.iterate(res::EigsolveResult{T1,T2,SuperOperatorQuantumObject}, ::Val{:vector_list}) where {T1,T2} =
-    ([QuantumObject(res.vectors[:, k], OperatorKet, res.dims) for k in 1:length(res.values)], Val(:vectors))
+    ([QuantumObject(res.vectors[:, k], OperatorKet, res.dimensions) for k in 1:length(res.values)], Val(:vectors))
 Base.iterate(res::EigsolveResult, ::Val{:vectors}) = (res.vectors, Val(:done))
 Base.iterate(res::EigsolveResult, ::Val{:done}) = nothing
 
 function Base.show(io::IO, res::EigsolveResult)
-    println(io, "EigsolveResult:   type=", res.type, "   dims=", res.dims)
+    println(io, "EigsolveResult:   type=", res.type, "   dims=", _get_dims_string(res.dimensions))
     println(io, "values:")
     show(io, MIME("text/plain"), res.values)
     print(io, "\n")
@@ -161,7 +165,7 @@ function _eigsolve(
     A,
     b::AbstractVector{T},
     type::ObjType,
-    dims::SVector,
+    dimensions::Union{Nothing,AbstractDimensions},
     k::Int = 1,
     m::Int = max(20, 2 * k + 1);
     tol::Real = 1e-8,
@@ -246,7 +250,7 @@ function _eigsolve(
     mul!(cache1, Vₘ, M(Uₘ * VR))
     vecs = cache1[:, 1:k]
 
-    return EigsolveResult(vals, vecs, type, dims, iter, numops, (iter < maxiter))
+    return EigsolveResult(vals, vecs, type, dimensions, iter, numops, (iter < maxiter))
 end
 
 @doc raw"""
@@ -283,7 +287,7 @@ function eigsolve(
         A.data;
         v0 = v0,
         type = A.type,
-        dims = A.dims,
+        dimensions = A.dimensions,
         sigma = sigma,
         k = k,
         krylovdim = krylovdim,
@@ -298,7 +302,7 @@ function eigsolve(
     A;
     v0::Union{Nothing,AbstractVector} = nothing,
     type::Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject} = nothing,
-    dims = SVector{0,Int}(),
+    dimensions = nothing,
     sigma::Union{Nothing,Real} = nothing,
     k::Int = 1,
     krylovdim::Int = max(20, 2 * k + 1),
@@ -311,10 +315,8 @@ function eigsolve(
     isH = ishermitian(A)
     v0 === nothing && (v0 = normalize!(rand(T, size(A, 1))))
 
-    dims = SVector(dims)
-
     if sigma === nothing
-        res = _eigsolve(A, v0, type, dims, k, krylovdim, tol = tol, maxiter = maxiter)
+        res = _eigsolve(A, v0, type, dimensions, k, krylovdim, tol = tol, maxiter = maxiter)
         vals = res.values
     else
         Aₛ = A - sigma * I
@@ -332,11 +334,11 @@ function eigsolve(
 
         Amap = EigsolveInverseMap(T, size(A), linsolve)
 
-        res = _eigsolve(Amap, v0, type, dims, k, krylovdim, tol = tol, maxiter = maxiter)
+        res = _eigsolve(Amap, v0, type, dimensions, k, krylovdim, tol = tol, maxiter = maxiter)
         vals = @. (1 + sigma * res.values) / res.values
     end
 
-    return EigsolveResult(vals, res.vectors, res.type, res.dims, res.iter, res.numops, res.converged)
+    return EigsolveResult(vals, res.vectors, res.type, res.dimensions, res.iter, res.numops, res.converged)
 end
 
 @doc raw"""
@@ -346,7 +348,7 @@ end
         c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         alg::OrdinaryDiffEqAlgorithm = Tsit5(),
         params::NamedTuple = NamedTuple(),
-        ρ0::AbstractMatrix = rand_dm(prod(H.dims)).data,
+        ρ0::AbstractMatrix = rand_dm(prod(H.dimensions)).data,
         k::Int = 1,
         krylovdim::Int = min(10, size(H, 1)),
         maxiter::Int = 200,
@@ -385,7 +387,7 @@ function eigsolve_al(
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
     params::NamedTuple = NamedTuple(),
-    ρ0::AbstractMatrix = rand_dm(prod(H.dims)).data,
+    ρ0::AbstractMatrix = rand_dm(prod(H.dimensions)).data,
     k::Int = 1,
     krylovdim::Int = min(10, size(H, 1)),
     maxiter::Int = 200,
@@ -396,7 +398,7 @@ function eigsolve_al(
     prob =
         mesolveProblem(
             L_evo,
-            QuantumObject(ρ0, type = Operator, dims = H.dims),
+            QuantumObject(ρ0, type = Operator, dims = H.dimensions),
             [zero(T), T];
             params = params,
             progress_bar = Val(false),
@@ -408,7 +410,7 @@ function eigsolve_al(
 
     Lmap = ArnoldiLindbladIntegratorMap(eltype(DT1), size(L_evo), integrator)
 
-    res = _eigsolve(Lmap, mat2vec(ρ0), L_evo.type, L_evo.dims, k, krylovdim, maxiter = maxiter, tol = eigstol)
+    res = _eigsolve(Lmap, mat2vec(ρ0), L_evo.type, L_evo.dimensions, k, krylovdim, maxiter = maxiter, tol = eigstol)
     # finish!(prog)
 
     vals = similar(res.values)
@@ -420,7 +422,7 @@ function eigsolve_al(
         @. vecs[:, i] = vec * exp(-1im * angle(vec[1]))
     end
 
-    return EigsolveResult(vals, vecs, res.type, res.dims, res.iter, res.numops, res.converged)
+    return EigsolveResult(vals, vecs, res.type, res.dimensions, res.iter, res.numops, res.converged)
 end
 
 @doc raw"""
@@ -464,7 +466,7 @@ function LinearAlgebra.eigen(
     E::mat2vec(sparse_to_dense(MT)) = F.values
     U::sparse_to_dense(MT) = F.vectors
 
-    return EigsolveResult(E, U, A.type, A.dims, 0, 0, true)
+    return EigsolveResult(E, U, A.type, A.dimensions, 0, 0, true)
 end
 
 @doc raw"""
diff --git a/src/qobj/functions.jl b/src/qobj/functions.jl
index 158514a4f..2b738a800 100644
--- a/src/qobj/functions.jl
+++ b/src/qobj/functions.jl
@@ -105,7 +105,7 @@ variance(O::QuantumObject{DT1,OperatorQuantumObject}, ψ::Vector{<:QuantumObject
 
 Converts a sparse QuantumObject to a dense QuantumObject.
 """
-sparse_to_dense(A::QuantumObject{<:AbstractVecOrMat}) = QuantumObject(sparse_to_dense(A.data), A.type, A.dims)
+sparse_to_dense(A::QuantumObject{<:AbstractVecOrMat}) = QuantumObject(sparse_to_dense(A.data), A.type, A.dimensions)
 sparse_to_dense(A::MT) where {MT<:AbstractSparseArray} = Array(A)
 for op in (:Transpose, :Adjoint)
     @eval sparse_to_dense(A::$op{T,<:AbstractSparseMatrix}) where {T<:BlasFloat} = Array(A)
@@ -132,7 +132,7 @@ sparse_to_dense(::Type{M}) where {M<:AbstractMatrix} = M
 Converts a dense QuantumObject to a sparse QuantumObject.
 """
 dense_to_sparse(A::QuantumObject{<:AbstractVecOrMat}, tol::Real = 1e-10) =
-    QuantumObject(dense_to_sparse(A.data, tol), A.type, A.dims)
+    QuantumObject(dense_to_sparse(A.data, tol), A.type, A.dimensions)
 function dense_to_sparse(A::MT, tol::Real = 1e-10) where {MT<:AbstractMatrix}
     idxs = findall(@. abs(A) > tol)
     row_indices = getindex.(idxs, 1)
@@ -180,12 +180,61 @@ julia> a.dims, O.dims
 ```
 """
 function LinearAlgebra.kron(
-    A::AbstractQuantumObject{DT1,OpType},
-    B::AbstractQuantumObject{DT2,OpType},
+    A::AbstractQuantumObject{DT1,OpType,<:Dimensions},
+    B::AbstractQuantumObject{DT2,OpType,<:Dimensions},
 ) where {DT1,DT2,OpType<:Union{KetQuantumObject,BraQuantumObject,OperatorQuantumObject}}
     QType = promote_op_type(A, B)
-    return QType(kron(A.data, B.data), A.type, vcat(A.dims, B.dims))
+    return QType(kron(A.data, B.data), A.type, Dimensions((A.dimensions.to..., B.dimensions.to...)))
 end
+
+# if A and B are both Operator but either one of them has GeneralDimensions
+for ADimType in (:Dimensions, :GeneralDimensions)
+    for BDimType in (:Dimensions, :GeneralDimensions)
+        if !(ADimType == BDimType == :Dimensions) # not for this case because it's already implemented
+            @eval begin
+                function LinearAlgebra.kron(
+                    A::AbstractQuantumObject{DT1,OperatorQuantumObject,<:$ADimType},
+                    B::AbstractQuantumObject{DT2,OperatorQuantumObject,<:$BDimType},
+                ) where {DT1,DT2}
+                    QType = promote_op_type(A, B)
+                    return QType(
+                        kron(A.data, B.data),
+                        Operator,
+                        GeneralDimensions(
+                            (get_dimensions_to(A)..., get_dimensions_to(B)...),
+                            (get_dimensions_from(A)..., get_dimensions_from(B)...),
+                        ),
+                    )
+                end
+            end
+        end
+    end
+end
+
+# if A and B are different type (must return Operator with GeneralDimensions)
+for AOpType in (:KetQuantumObject, :BraQuantumObject, :OperatorQuantumObject)
+    for BOpType in (:KetQuantumObject, :BraQuantumObject, :OperatorQuantumObject)
+        if (AOpType != BOpType)
+            @eval begin
+                function LinearAlgebra.kron(
+                    A::AbstractQuantumObject{DT1,$AOpType},
+                    B::AbstractQuantumObject{DT2,$BOpType},
+                ) where {DT1,DT2}
+                    QType = promote_op_type(A, B)
+                    return QType(
+                        kron(A.data, B.data),
+                        Operator,
+                        GeneralDimensions(
+                            (get_dimensions_to(A)..., get_dimensions_to(B)...),
+                            (get_dimensions_from(A)..., get_dimensions_from(B)...),
+                        ),
+                    )
+                end
+            end
+        end
+    end
+end
+
 LinearAlgebra.kron(A::AbstractQuantumObject) = A
 function LinearAlgebra.kron(A::Vector{<:AbstractQuantumObject})
     @warn "`tensor(A)` or `kron(A)` with `A` is a `Vector` can hurt performance. Try to use `tensor(A...)` or `kron(A...)` instead."
@@ -208,7 +257,7 @@ end
 Convert a quantum object from vector ([`OperatorKetQuantumObject`](@ref)-type) to matrix ([`OperatorQuantumObject`](@ref)-type)
 """
 vec2mat(A::QuantumObject{<:AbstractArray{T},OperatorKetQuantumObject}) where {T} =
-    QuantumObject(vec2mat(A.data), Operator, A.dims)
+    QuantumObject(vec2mat(A.data), Operator, A.dimensions)
 
 @doc raw"""
     mat2vec(A::QuantumObject)
@@ -216,7 +265,7 @@ vec2mat(A::QuantumObject{<:AbstractArray{T},OperatorKetQuantumObject}) where {T}
 Convert a quantum object from matrix ([`OperatorQuantumObject`](@ref)-type) to vector ([`OperatorKetQuantumObject`](@ref)-type)
 """
 mat2vec(A::QuantumObject{<:AbstractArray{T},OperatorQuantumObject}) where {T} =
-    QuantumObject(mat2vec(A.data), OperatorKet, A.dims)
+    QuantumObject(mat2vec(A.data), OperatorKet, A.dimensions)
 
 @doc raw"""
     mat2vec(A::AbstractMatrix)
diff --git a/src/qobj/operators.jl b/src/qobj/operators.jl
index d2a0d8fe5..9b0e1de5f 100644
--- a/src/qobj/operators.jl
+++ b/src/qobj/operators.jl
@@ -19,7 +19,7 @@ Returns a random unitary [`QuantumObject`](@ref).
 
 The `dimensions` can be either the following types:
 - `dimensions::Int`: Number of basis states in the Hilbert space.
-- `dimensions::Union{AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+- `dimensions::Union{Dimensions,AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
 
 The `distribution` specifies which of the method used to obtain the unitary matrix:
 - `:haar`: Haar random unitary matrix using the algorithm from reference 1
@@ -33,9 +33,9 @@ The `distribution` specifies which of the method used to obtain the unitary matr
 """
 rand_unitary(dimensions::Int, distribution::Union{Symbol,Val} = Val(:haar)) =
     rand_unitary(SVector(dimensions), makeVal(distribution))
-rand_unitary(dimensions::Union{AbstractVector{Int},Tuple}, distribution::Union{Symbol,Val} = Val(:haar)) =
+rand_unitary(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}, distribution::Union{Symbol,Val} = Val(:haar)) =
     rand_unitary(dimensions, makeVal(distribution))
-function rand_unitary(dimensions::Union{AbstractVector{Int},Tuple}, ::Val{:haar})
+function rand_unitary(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}, ::Val{:haar})
     N = prod(dimensions)
 
     # generate N x N matrix Z of complex standard normal random variates
@@ -50,7 +50,7 @@ function rand_unitary(dimensions::Union{AbstractVector{Int},Tuple}, ::Val{:haar}
     Λ ./= abs.(Λ) # rescaling the elements
     return QuantumObject(sparse_to_dense(Q * Diagonal(Λ)); type = Operator, dims = dimensions)
 end
-function rand_unitary(dimensions::Union{AbstractVector{Int},Tuple}, ::Val{:exp})
+function rand_unitary(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}, ::Val{:exp})
     N = prod(dimensions)
 
     # generate N x N matrix Z of complex standard normal random variates
@@ -61,7 +61,7 @@ function rand_unitary(dimensions::Union{AbstractVector{Int},Tuple}, ::Val{:exp})
 
     return sparse_to_dense(exp(-1.0im * H))
 end
-rand_unitary(dimensions::Union{AbstractVector{Int},Tuple}, ::Val{T}) where {T} =
+rand_unitary(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}, ::Val{T}) where {T} =
     throw(ArgumentError("Invalid distribution: $(T)"))
 
 @doc raw"""
@@ -432,12 +432,16 @@ Note that `type` can only be either [`Operator`](@ref) or [`SuperOperator`](@ref
 !!! note
     `qeye` is a synonym of `eye`.
 """
-eye(
+function eye(
     N::Int;
     type::ObjType = Operator,
     dims = nothing,
-) where {ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    QuantumObject(Diagonal(ones(ComplexF64, N)); type = type, dims = dims)
+) where {ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
+    if dims isa Nothing
+        dims = isa(type, OperatorQuantumObject) ? N : isqrt(N)
+    end
+    return QuantumObject(Diagonal(ones(ComplexF64, N)); type = type, dims = dims)
+end
 
 @doc raw"""
     fdestroy(N::Union{Int,Val}, j::Int)
@@ -498,7 +502,8 @@ end
 
 Generates the projection operator ``\hat{O} = |i \rangle\langle j|`` with Hilbert space dimension `N`.
 """
-projection(N::Int, i::Int, j::Int) = QuantumObject(sparse([i + 1], [j + 1], [1.0 + 0.0im], N, N), type = Operator)
+projection(N::Int, i::Int, j::Int) =
+    QuantumObject(sparse([i + 1], [j + 1], [1.0 + 0.0im], N, N), type = Operator, dims = N)
 
 @doc raw"""
     tunneling(N::Int, m::Int=1; sparse::Union{Bool,Val{<:Bool}}=Val(false))
@@ -534,7 +539,7 @@ Generates a discrete Fourier transform matrix ``\hat{F}_N`` for [Quantum Fourier
 
 The `dimensions` can be either the following types:
 - `dimensions::Int`: Number of basis states in the Hilbert space.
-- `dimensions::Union{AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+- `dimensions::Union{Dimensions,AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
 
 ``N`` represents the total dimension, and therefore the matrix is defined as
 
@@ -555,7 +560,7 @@ where ``\omega = \exp(\frac{2 \pi i}{N})``.
     It is highly recommended to use `qft(dimensions)` with `dimensions` as `Tuple` or `SVector` to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 qft(dimensions::Int) = QuantumObject(_qft_op(dimensions), Operator, dimensions)
-qft(dimensions::Union{AbstractVector{T},Tuple}) where {T} =
+qft(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}) =
     QuantumObject(_qft_op(prod(dimensions)), Operator, dimensions)
 function _qft_op(N::Int)
     ω = exp(2.0im * π / N)
diff --git a/src/qobj/quantum_object.jl b/src/qobj/quantum_object.jl
index 0b0269c9f..a2ba3ffdf 100644
--- a/src/qobj/quantum_object.jl
+++ b/src/qobj/quantum_object.jl
@@ -9,13 +9,16 @@ It also implements the fundamental functions in Julia standard library:
 export QuantumObject
 
 @doc raw"""
-    struct QuantumObject{MT<:AbstractArray,ObjType<:QuantumObjectType,N}
-        data::MT
+    struct QuantumObject{DataType<:AbstractArray,ObjType<:QuantumObjectType,DimType<:AbstractDimensions} <: AbstractQuantumObject{DataType,ObjType,DimType}
+        data::DataType
         type::ObjType
-        dims::SVector{N, Int}
+        dimensions::DimType
     end
 
-Julia struct representing any quantum objects.
+Julia structure representing any time-independent quantum objects. For time-dependent cases, see [`QuantumObjectEvolution`](@ref).
+
+!!! note "`dims` property"
+    For a given `H::QuantumObject`, `H.dims` or `getproperty(H, :dims)` returns its `dimensions` in the type of integer-vector.
 
 # Examples
 
@@ -32,28 +35,31 @@ Quantum Object:   type=Operator   dims=[20]   size=(20, 20)   ishermitian=false
 
 julia> a isa QuantumObject
 true
+
+julia> a.dims
+1-element SVector{1, Int64} with indices SOneTo(1):
+ 20
+
+julia> a.dimensions
+Dimensions{1, Tuple{Space}}((Space(20),))
 ```
 """
-struct QuantumObject{MT<:AbstractArray,ObjType<:QuantumObjectType,N} <: AbstractQuantumObject{MT,ObjType,N}
-    data::MT
+struct QuantumObject{DataType<:AbstractArray,ObjType<:QuantumObjectType,DimType<:AbstractDimensions} <:
+       AbstractQuantumObject{DataType,ObjType,DimType}
+    data::DataType
     type::ObjType
-    dims::SVector{N,Int}
+    dimensions::DimType
 
     function QuantumObject(data::MT, type::ObjType, dims) where {MT<:AbstractArray,ObjType<:QuantumObjectType}
-        _check_dims(dims)
+        dimensions = _gen_dimensions(dims)
 
         _size = _get_size(data)
-        _check_QuantumObject(type, dims, _size[1], _size[2])
-
-        N = length(dims)
+        _check_QuantumObject(type, dimensions, _size[1], _size[2])
 
-        return new{MT,ObjType,N}(data, type, SVector{N,Int}(dims))
+        return new{MT,ObjType,typeof(dimensions)}(data, type, dimensions)
     end
 end
 
-QuantumObject(A::AbstractArray, type::ObjType, dims::Integer) where {ObjType<:QuantumObjectType} =
-    QuantumObject(A, type, SVector{1,Int}(dims))
-
 @doc raw"""
     Qobj(A::AbstractArray; type = nothing, dims = nothing)
     QuantumObject(A::AbstractArray; type = nothing, dims = nothing)
@@ -81,10 +87,14 @@ function QuantumObject(
     end
 
     if dims isa Nothing
-        if type isa OperatorQuantumObject || type isa BraQuantumObject
-            dims = SVector{1,Int}(_size[2])
+        if type isa BraQuantumObject
+            dims = Dimensions(_size[2])
+        elseif type isa OperatorQuantumObject
+            dims =
+                (_size[1] == _size[2]) ? Dimensions(_size[1]) :
+                GeneralDimensions(SVector{2}(SVector{1}(_size[1]), SVector{1}(_size[2])))
         elseif type isa SuperOperatorQuantumObject || type isa OperatorBraQuantumObject
-            dims = SVector{1,Int}(isqrt(_size[2]))
+            dims = Dimensions(isqrt(_size[2]))
         end
     end
 
@@ -105,9 +115,9 @@ function QuantumObject(
     if dims isa Nothing
         _size = _get_size(A)
         if type isa KetQuantumObject
-            dims = SVector{1,Int}(_size[1])
+            dims = Dimensions(_size[1])
         elseif type isa OperatorKetQuantumObject
-            dims = SVector{1,Int}(isqrt(_size[1]))
+            dims = Dimensions(isqrt(_size[1]))
         end
     end
 
@@ -125,11 +135,12 @@ end
 function QuantumObject(
     A::QuantumObject{<:AbstractArray{T,N}};
     type::ObjType = A.type,
-    dims = A.dims,
+    dims = A.dimensions,
 ) where {T,N,ObjType<:QuantumObjectType}
     _size = N == 1 ? (length(A), 1) : size(A)
-    _check_QuantumObject(type, dims, _size[1], _size[2])
-    return QuantumObject(copy(A.data), type, dims)
+    dimensions = _gen_dimensions(dims)
+    _check_QuantumObject(type, dimensions, _size[1], _size[2])
+    return QuantumObject(copy(A.data), type, dimensions)
 end
 
 function Base.show(
@@ -146,7 +157,15 @@ function Base.show(
     },
 }
     op_data = QO.data
-    println(io, "\nQuantum Object:   type=", QO.type, "   dims=", QO.dims, "   size=", size(op_data))
+    println(
+        io,
+        "\nQuantum Object:   type=",
+        QO.type,
+        "   dims=",
+        _get_dims_string(QO.dimensions),
+        "   size=",
+        size(op_data),
+    )
     return show(io, MIME("text/plain"), op_data)
 end
 
@@ -157,7 +176,7 @@ function Base.show(io::IO, QO::QuantumObject)
         "\nQuantum Object:   type=",
         QO.type,
         "   dims=",
-        QO.dims,
+        _get_dims_string(QO.dimensions),
         "   size=",
         size(op_data),
         "   ishermitian=",
@@ -166,15 +185,16 @@ function Base.show(io::IO, QO::QuantumObject)
     return show(io, MIME("text/plain"), op_data)
 end
 
-Base.real(x::QuantumObject) = QuantumObject(real(x.data), x.type, x.dims)
-Base.imag(x::QuantumObject) = QuantumObject(imag(x.data), x.type, x.dims)
+Base.real(x::QuantumObject) = QuantumObject(real(x.data), x.type, x.dimensions)
+Base.imag(x::QuantumObject) = QuantumObject(imag(x.data), x.type, x.dimensions)
 
-SparseArrays.sparse(A::QuantumObject{<:AbstractArray{T}}) where {T} = QuantumObject(sparse(A.data), A.type, A.dims)
+SparseArrays.sparse(A::QuantumObject{<:AbstractArray{T}}) where {T} =
+    QuantumObject(sparse(A.data), A.type, A.dimensions)
 SparseArrays.nnz(A::QuantumObject{<:AbstractSparseArray}) = nnz(A.data)
 SparseArrays.nonzeros(A::QuantumObject{<:AbstractSparseArray}) = nonzeros(A.data)
 SparseArrays.rowvals(A::QuantumObject{<:AbstractSparseArray}) = rowvals(A.data)
 SparseArrays.droptol!(A::QuantumObject{<:AbstractSparseArray}, tol::Real) = (droptol!(A.data, tol); return A)
-SparseArrays.dropzeros(A::QuantumObject{<:AbstractSparseArray}) = QuantumObject(dropzeros(A.data), A.type, A.dims)
+SparseArrays.dropzeros(A::QuantumObject{<:AbstractSparseArray}) = QuantumObject(dropzeros(A.data), A.type, A.dimensions)
 SparseArrays.dropzeros!(A::QuantumObject{<:AbstractSparseArray}) = (dropzeros!(A.data); return A)
 
 @doc raw"""
@@ -191,7 +211,7 @@ SciMLOperators.cache_operator(
     L::AbstractQuantumObject{DT,OpType},
     u::AbstractVector,
 ) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    get_typename_wrapper(L)(cache_operator(L.data, sparse_to_dense(similar(u))), L.type, L.dims)
+    get_typename_wrapper(L)(cache_operator(L.data, sparse_to_dense(similar(u))), L.type, L.dimensions)
 
 function SciMLOperators.cache_operator(
     L::AbstractQuantumObject{DT1,OpType},
@@ -202,7 +222,7 @@ function SciMLOperators.cache_operator(
     OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     SType<:Union{KetQuantumObject,OperatorKetQuantumObject},
 }
-    check_dims(L, u)
+    check_dimensions(L, u)
 
     if isoper(L) && isoperket(u)
         throw(ArgumentError("The input state `u` must be a Ket if `L` is an Operator."))
@@ -213,14 +233,17 @@ function SciMLOperators.cache_operator(
 end
 
 # data type conversions
-Base.Vector(A::QuantumObject{<:AbstractVector}) = QuantumObject(Vector(A.data), A.type, A.dims)
-Base.Vector{T}(A::QuantumObject{<:AbstractVector}) where {T<:Number} = QuantumObject(Vector{T}(A.data), A.type, A.dims)
-Base.Matrix(A::QuantumObject{<:AbstractMatrix}) = QuantumObject(Matrix(A.data), A.type, A.dims)
-Base.Matrix{T}(A::QuantumObject{<:AbstractMatrix}) where {T<:Number} = QuantumObject(Matrix{T}(A.data), A.type, A.dims)
-SparseArrays.SparseVector(A::QuantumObject{<:AbstractVector}) = QuantumObject(SparseVector(A.data), A.type, A.dims)
+Base.Vector(A::QuantumObject{<:AbstractVector}) = QuantumObject(Vector(A.data), A.type, A.dimensions)
+Base.Vector{T}(A::QuantumObject{<:AbstractVector}) where {T<:Number} =
+    QuantumObject(Vector{T}(A.data), A.type, A.dimensions)
+Base.Matrix(A::QuantumObject{<:AbstractMatrix}) = QuantumObject(Matrix(A.data), A.type, A.dimensions)
+Base.Matrix{T}(A::QuantumObject{<:AbstractMatrix}) where {T<:Number} =
+    QuantumObject(Matrix{T}(A.data), A.type, A.dimensions)
+SparseArrays.SparseVector(A::QuantumObject{<:AbstractVector}) =
+    QuantumObject(SparseVector(A.data), A.type, A.dimensions)
 SparseArrays.SparseVector{T}(A::QuantumObject{<:SparseVector}) where {T<:Number} =
-    QuantumObject(SparseVector{T}(A.data), A.type, A.dims)
+    QuantumObject(SparseVector{T}(A.data), A.type, A.dimensions)
 SparseArrays.SparseMatrixCSC(A::QuantumObject{<:AbstractMatrix}) =
-    QuantumObject(SparseMatrixCSC(A.data), A.type, A.dims)
+    QuantumObject(SparseMatrixCSC(A.data), A.type, A.dimensions)
 SparseArrays.SparseMatrixCSC{T}(A::QuantumObject{<:SparseMatrixCSC}) where {T<:Number} =
-    QuantumObject(SparseMatrixCSC{T}(A.data), A.type, A.dims)
+    QuantumObject(SparseMatrixCSC{T}(A.data), A.type, A.dimensions)
diff --git a/src/qobj/quantum_object_base.jl b/src/qobj/quantum_object_base.jl
index df259307c..cec0992b3 100644
--- a/src/qobj/quantum_object_base.jl
+++ b/src/qobj/quantum_object_base.jl
@@ -14,7 +14,7 @@ export QuantumObjectType,
 export Bra, Ket, Operator, OperatorBra, OperatorKet, SuperOperator
 
 @doc raw"""
-    abstract type AbstractQuantumObject{DataType,ObjType,N}
+    abstract type AbstractQuantumObject{DataType,ObjType,DimType}
 
 Abstract type for all quantum objects like [`QuantumObject`](@ref) and [`QuantumObjectEvolution`](@ref).
 
@@ -24,7 +24,7 @@ julia> sigmax() isa AbstractQuantumObject
 true
 ```
 """
-abstract type AbstractQuantumObject{DataType,ObjType,N} end
+abstract type AbstractQuantumObject{DataType,ObjType,DimType} end
 
 abstract type QuantumObjectType end
 
@@ -152,68 +152,138 @@ Returns the length of the matrix or vector corresponding to the [`AbstractQuantu
 Base.length(A::AbstractQuantumObject) = length(A.data)
 
 Base.isequal(A::AbstractQuantumObject, B::AbstractQuantumObject) =
-    isequal(A.type, B.type) && isequal(A.dims, B.dims) && isequal(A.data, B.data)
+    isequal(A.type, B.type) && isequal(A.dimensions, B.dimensions) && isequal(A.data, B.data)
 Base.isapprox(A::AbstractQuantumObject, B::AbstractQuantumObject; kwargs...) =
-    isequal(A.type, B.type) && isequal(A.dims, B.dims) && isapprox(A.data, B.data; kwargs...)
+    isequal(A.type, B.type) && isequal(A.dimensions, B.dimensions) && isapprox(A.data, B.data; kwargs...)
 Base.:(==)(A::AbstractQuantumObject, B::AbstractQuantumObject) =
-    (A.type == B.type) && (A.dims == B.dims) && (A.data == B.data)
+    (A.type == B.type) && (A.dimensions == B.dimensions) && (A.data == B.data)
 
-function check_dims(A::AbstractQuantumObject, B::AbstractQuantumObject)
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects don't have the same Hilbert dimension."))
+function check_dimensions(dimensions_list::NTuple{N,AbstractDimensions}) where {N}
+    allequal(dimensions_list) ||
+        throw(DimensionMismatch("The quantum objects should have the same Hilbert `dimensions`."))
     return nothing
 end
-
-function _check_dims(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Integer,N}
-    _non_static_array_warning("dims", dims)
-    return (all(>(0), dims) && length(dims) > 0) ||
-           throw(DomainError(dims, "The argument dims must be of non-zero length and contain only positive integers."))
-end
-_check_dims(dims::Any) = throw(
-    ArgumentError(
-        "The argument dims must be a Tuple or a StaticVector of non-zero length and contain only positive integers.",
+check_dimensions(Qobj_tuple::NTuple{N,AbstractQuantumObject}) where {N} =
+    check_dimensions(getfield.(Qobj_tuple, :dimensions))
+check_dimensions(A::AbstractQuantumObject...) = check_dimensions(A)
+
+_check_QuantumObject(
+    type::ObjType,
+    dimensions::GeneralDimensions,
+    m::Int,
+    n::Int,
+) where {
+    ObjType<:Union{
+        KetQuantumObject,
+        BraQuantumObject,
+        SuperOperatorQuantumObject,
+        OperatorBraQuantumObject,
+        OperatorKetQuantumObject,
+    },
+} = throw(
+    DomainError(
+        _get_dims_string(dimensions),
+        "The given `dims` is not compatible with type = $type, should be a single list of integers.",
     ),
 )
 
-function _check_QuantumObject(type::KetQuantumObject, dims, m::Int, n::Int)
+function _check_QuantumObject(type::KetQuantumObject, dimensions::Dimensions, m::Int, n::Int)
     (n != 1) && throw(DomainError((m, n), "The size of the array is not compatible with Ket"))
-    (prod(dims) != m) && throw(DimensionMismatch("Ket with dims = $(dims) does not fit the array size = $((m, n))."))
+    (prod(dimensions) != m) && throw(
+        DimensionMismatch("Ket with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n))."),
+    )
     return nothing
 end
 
-function _check_QuantumObject(type::BraQuantumObject, dims, m::Int, n::Int)
+function _check_QuantumObject(type::BraQuantumObject, dimensions::Dimensions, m::Int, n::Int)
     (m != 1) && throw(DomainError((m, n), "The size of the array is not compatible with Bra"))
-    (prod(dims) != n) && throw(DimensionMismatch("Bra with dims = $(dims) does not fit the array size = $((m, n))."))
+    (prod(dimensions) != n) && throw(
+        DimensionMismatch("Bra with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n))."),
+    )
+    return nothing
+end
+
+function _check_QuantumObject(type::OperatorQuantumObject, dimensions::Dimensions, m::Int, n::Int)
+    L = prod(dimensions)
+    (L == m == n) || throw(
+        DimensionMismatch(
+            "Operator with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
+        ),
+    )
     return nothing
 end
 
-function _check_QuantumObject(type::OperatorQuantumObject, dims, m::Int, n::Int)
-    (m != n) && throw(DomainError((m, n), "The size of the array is not compatible with Operator"))
-    (prod(dims) != m) &&
-        throw(DimensionMismatch("Operator with dims = $(dims) does not fit the array size = $((m, n))."))
+function _check_QuantumObject(type::OperatorQuantumObject, dimensions::GeneralDimensions, m::Int, n::Int)
+    ((m == 1) || (n == 1)) && throw(DomainError((m, n), "The size of the array is not compatible with Operator"))
+    ((prod(dimensions.to) != m) || (prod(dimensions.from) != n)) && throw(
+        DimensionMismatch(
+            "Operator with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
+        ),
+    )
     return nothing
 end
 
-function _check_QuantumObject(type::SuperOperatorQuantumObject, dims, m::Int, n::Int)
+function _check_QuantumObject(type::SuperOperatorQuantumObject, dimensions::Dimensions, m::Int, n::Int)
     (m != n) && throw(DomainError((m, n), "The size of the array is not compatible with SuperOperator"))
-    (prod(dims) != sqrt(m)) &&
-        throw(DimensionMismatch("SuperOperator with dims = $(dims) does not fit the array size = $((m, n))."))
+    (prod(dimensions) != sqrt(m)) && throw(
+        DimensionMismatch(
+            "SuperOperator with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
+        ),
+    )
     return nothing
 end
 
-function _check_QuantumObject(type::OperatorKetQuantumObject, dims, m::Int, n::Int)
+function _check_QuantumObject(type::OperatorKetQuantumObject, dimensions::Dimensions, m::Int, n::Int)
     (n != 1) && throw(DomainError((m, n), "The size of the array is not compatible with OperatorKet"))
-    (prod(dims) != sqrt(m)) &&
-        throw(DimensionMismatch("OperatorKet with dims = $(dims) does not fit the array size = $((m, n))."))
+    (prod(dimensions) != sqrt(m)) && throw(
+        DimensionMismatch(
+            "OperatorKet with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
+        ),
+    )
     return nothing
 end
 
-function _check_QuantumObject(type::OperatorBraQuantumObject, dims, m::Int, n::Int)
+function _check_QuantumObject(type::OperatorBraQuantumObject, dimensions::Dimensions, m::Int, n::Int)
     (m != 1) && throw(DomainError((m, n), "The size of the array is not compatible with OperatorBra"))
-    (prod(dims) != sqrt(n)) &&
-        throw(DimensionMismatch("OperatorBra with dims = $(dims) does not fit the array size = $((m, n))."))
+    (prod(dimensions) != sqrt(n)) && throw(
+        DimensionMismatch(
+            "OperatorBra with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
+        ),
+    )
     return nothing
 end
 
+function Base.getproperty(A::AbstractQuantumObject, key::Symbol)
+    # a comment here to avoid bad render by JuliaFormatter
+    if key === :dims
+        return dimensions_to_dims(getfield(A, :dimensions))
+    else
+        return getfield(A, key)
+    end
+end
+
+# this returns `to` in GeneralDimensions representation
+get_dimensions_to(A::AbstractQuantumObject{DT,KetQuantumObject,<:Dimensions{N}}) where {DT,N} = A.dimensions.to
+get_dimensions_to(A::AbstractQuantumObject{DT,BraQuantumObject,<:Dimensions{N}}) where {DT,N} = space_one_list(N)
+get_dimensions_to(A::AbstractQuantumObject{DT,OperatorQuantumObject,<:Dimensions{N}}) where {DT,N} = A.dimensions.to
+get_dimensions_to(A::AbstractQuantumObject{DT,OperatorQuantumObject,<:GeneralDimensions{N}}) where {DT,N} =
+    A.dimensions.to
+get_dimensions_to(
+    A::AbstractQuantumObject{DT,ObjType,<:Dimensions{N}},
+) where {DT,ObjType<:Union{SuperOperatorQuantumObject,OperatorBraQuantumObject,OperatorKetQuantumObject},N} =
+    A.dimensions.to
+
+# this returns `from` in GeneralDimensions representation
+get_dimensions_from(A::AbstractQuantumObject{DT,KetQuantumObject,<:Dimensions{N}}) where {DT,N} = space_one_list(N)
+get_dimensions_from(A::AbstractQuantumObject{DT,BraQuantumObject,<:Dimensions{N}}) where {DT,N} = A.dimensions.to
+get_dimensions_from(A::AbstractQuantumObject{DT,OperatorQuantumObject,<:Dimensions{N}}) where {DT,N} = A.dimensions.to
+get_dimensions_from(A::AbstractQuantumObject{DT,OperatorQuantumObject,<:GeneralDimensions{N}}) where {DT,N} =
+    A.dimensions.from
+get_dimensions_from(
+    A::AbstractQuantumObject{DT,ObjType,<:Dimensions{N}},
+) where {DT,ObjType<:Union{SuperOperatorQuantumObject,OperatorBraQuantumObject,OperatorKetQuantumObject},N} =
+    A.dimensions.to
+
 # functions for getting Float or Complex element type
 _FType(A::AbstractQuantumObject) = _FType(eltype(A))
 _CType(A::AbstractQuantumObject) = _CType(eltype(A))
diff --git a/src/qobj/quantum_object_evo.jl b/src/qobj/quantum_object_evo.jl
index ef7f0be60..f055c61ae 100644
--- a/src/qobj/quantum_object_evo.jl
+++ b/src/qobj/quantum_object_evo.jl
@@ -5,10 +5,10 @@ This file defines the QuantumObjectEvolution (QobjEvo) structure.
 export QuantumObjectEvolution
 
 @doc raw"""
-    struct QuantumObjectEvolution{DT<:AbstractSciMLOperator,ObjType<:QuantumObjectType,N} <: AbstractQuantumObject
-        data::DT
+    struct QuantumObjectEvolution{DataType<:AbstractSciMLOperator,ObjType<:QuantumObjectType,DimType<:AbstractDimensions} <: AbstractQuantumObject{DataType,ObjType,DimType}
+        data::DataType
         type::ObjType
-        dims::SVector{N,Int}
+        dimensions::DimType
     end
 
 Julia struct representing any time-dependent quantum object. The `data` field is a `AbstractSciMLOperator` object that represents the time-dependent quantum object. It can be seen as
@@ -17,7 +17,12 @@ Julia struct representing any time-dependent quantum object. The `data` field is
 \hat{O}(t) = \sum_{i} c_i(p, t) \hat{O}_i
 ```
 
-where ``c_i(p, t)`` is a function that depends on the parameters `p` and time `t`, and ``\hat{O}_i`` are the operators that form the quantum object. The `type` field is the type of the quantum object, and the `dims` field is the dimensions of the quantum object. For more information about `type` and `dims`, see [`QuantumObject`](@ref). For more information about `AbstractSciMLOperator`, see the [SciML](https://docs.sciml.ai/SciMLOperators/stable/) documentation.
+where ``c_i(p, t)`` is a function that depends on the parameters `p` and time `t`, and ``\hat{O}_i`` are the operators that form the quantum object. 
+
+For time-independent cases, see [`QuantumObject`](@ref), and for more information about `AbstractSciMLOperator`, see the [SciML](https://docs.sciml.ai/SciMLOperators/stable/) documentation.
+
+!!! note "`dims` property"
+    For a given `H::QuantumObjectEvolution`, `H.dims` or `getproperty(H, :dims)` returns its `dimensions` in the type of integer-vector.
 
 # Examples
 This operator can be initialized in the same way as the QuTiP `QobjEvo` object. For example
@@ -104,13 +109,13 @@ Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=fal
 ```
 """
 struct QuantumObjectEvolution{
-    DT<:AbstractSciMLOperator,
+    DataType<:AbstractSciMLOperator,
     ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
-    N,
-} <: AbstractQuantumObject{DT,ObjType,N}
-    data::DT
+    DimType<:AbstractDimensions,
+} <: AbstractQuantumObject{DataType,ObjType,DimType}
+    data::DataType
     type::ObjType
-    dims::SVector{N,Int}
+    dimensions::DimType
 
     function QuantumObjectEvolution(
         data::DT,
@@ -120,14 +125,12 @@ struct QuantumObjectEvolution{
         (type == Operator || type == SuperOperator) ||
             throw(ArgumentError("The type $type is not supported for QuantumObjectEvolution."))
 
-        _check_dims(dims)
+        dimensions = _gen_dimensions(dims)
 
         _size = _get_size(data)
-        _check_QuantumObject(type, dims, _size[1], _size[2])
-
-        N = length(dims)
+        _check_QuantumObject(type, dimensions, _size[1], _size[2])
 
-        return new{DT,ObjType,N}(data, type, SVector{N,Int}(dims))
+        return new{DT,ObjType,typeof(dimensions)}(data, type, dimensions)
     end
 end
 
@@ -138,7 +141,7 @@ function Base.show(io::IO, QO::QuantumObjectEvolution)
         "\nQuantum Object Evo.:   type=",
         QO.type,
         "   dims=",
-        QO.dims,
+        _get_dims_string(QO.dimensions),
         "   size=",
         size(op_data),
         "   ishermitian=",
@@ -149,10 +152,6 @@ function Base.show(io::IO, QO::QuantumObjectEvolution)
     return show(io, MIME("text/plain"), op_data)
 end
 
-function QuantumObjectEvolution(data::AbstractSciMLOperator, type::QuantumObjectType, dims::Integer)
-    return QuantumObjectEvolution(data, type, SVector{1,Int}(dims))
-end
-
 @doc raw"""
     QobjEvo(data::AbstractSciMLOperator; type::QuantumObjectType = Operator, dims = nothing)
     QuantumObjectEvolution(data::AbstractSciMLOperator; type::QuantumObjectType = Operator, dims = nothing)
@@ -165,10 +164,12 @@ function QuantumObjectEvolution(data::AbstractSciMLOperator; type::QuantumObject
     _size = _get_size(data)
 
     if dims isa Nothing
-        if type == Operator
-            dims = SVector{1,Int}(_size[2])
-        elseif type == SuperOperator
-            dims = SVector{1,Int}(isqrt(_size[2]))
+        if type isa OperatorQuantumObject
+            dims =
+                (_size[1] == _size[2]) ? Dimensions(_size[1]) :
+                GeneralDimensions(SVector{2}(SVector{1}(_size[1]), SVector{1}(_size[2])))
+        elseif type isa SuperOperatorQuantumObject
+            dims = Dimensions(isqrt(_size[2]))
         end
     end
 
@@ -274,7 +275,7 @@ function QuantumObjectEvolution(
     type::Union{Nothing,QuantumObjectType} = nothing,
 )
     op, data = _QobjEvo_generate_data(op_func_list, α)
-    dims = op.dims
+    dims = op.dimensions
     if type isa Nothing
         type = op.type
     end
@@ -339,7 +340,7 @@ function QuantumObjectEvolution(
     if type isa Nothing
         type = op.type
     end
-    return QuantumObjectEvolution(_make_SciMLOperator(op, α), type, op.dims)
+    return QuantumObjectEvolution(_make_SciMLOperator(op, α), type, op.dimensions)
 end
 
 function QuantumObjectEvolution(
@@ -357,9 +358,9 @@ function QuantumObjectEvolution(
         )
     end
     if α isa Nothing
-        return QuantumObjectEvolution(op.data, type, op.dims)
+        return QuantumObjectEvolution(op.data, type, op.dimensions)
     end
-    return QuantumObjectEvolution(α * op.data, type, op.dims)
+    return QuantumObjectEvolution(α * op.data, type, op.dimensions)
 end
 
 #=
@@ -397,7 +398,7 @@ Parse the `op_func_list` and generate the data for the `QuantumObjectEvolution`
 
             op = :(op_func_list[$i][1])
             data_type = op_type.parameters[1]
-            dims_expr = (dims_expr..., :($op.dims))
+            dims_expr = (dims_expr..., :($op.dimensions))
             func_methods_expr = (func_methods_expr..., :(methods(op_func_list[$i][2], [Any, Real]))) # [Any, Real] means each func must accept 2 arguments
             if i == 1
                 first_op = :($op)
@@ -409,7 +410,7 @@ Parse the `op_func_list` and generate the data for the `QuantumObjectEvolution`
                 throw(ArgumentError("The element must be a Operator or SuperOperator."))
 
             data_type = op_type.parameters[1]
-            dims_expr = (dims_expr..., :(op_func_list[$i].dims))
+            dims_expr = (dims_expr..., :(op_func_list[$i].dimensions))
             if i == 1
                 first_op = :(op_func_list[$i])
             end
@@ -507,8 +508,7 @@ function (A::QuantumObjectEvolution)(
     p,
     t,
 ) where {DT1,DT2,QobjType<:Union{KetQuantumObject,OperatorKetQuantumObject}}
-    check_dims(ψout, ψin)
-    check_dims(ψout, A)
+    check_dimensions(A, ψout, ψin)
 
     if isoper(A) && isoperket(ψin)
         throw(ArgumentError("The input state must be a Ket if the QuantumObjectEvolution object is an Operator."))
@@ -535,7 +535,7 @@ function (A::QuantumObjectEvolution)(
     p,
     t,
 ) where {DT,QobjType<:Union{KetQuantumObject,OperatorKetQuantumObject}}
-    ψout = QuantumObject(similar(ψ.data), ψ.type, ψ.dims)
+    ψout = QuantumObject(similar(ψ.data), ψ.type, ψ.dimensions)
     return A(ψout, ψ, p, t)
 end
 
@@ -554,7 +554,7 @@ Calculate the time-dependent [`QuantumObjectEvolution`](@ref) object `A` at time
 function (A::QuantumObjectEvolution)(p, t)
     # We put 0 in the place of `u` because the time-dependence doesn't depend on the state
     update_coefficients!(A.data, 0, p, t)
-    return QuantumObject(concretize(A.data), A.type, A.dims)
+    return QuantumObject(concretize(A.data), A.type, A.dimensions)
 end
 
 (A::QuantumObjectEvolution)(t) = A(nothing, t)
diff --git a/src/qobj/space.jl b/src/qobj/space.jl
new file mode 100644
index 000000000..0b90ecf55
--- /dev/null
+++ b/src/qobj/space.jl
@@ -0,0 +1,41 @@
+#=
+This file defines the Hilbert space structure.
+=#
+
+export AbstractSpace, Space
+
+abstract type AbstractSpace end
+
+@doc raw"""
+    struct Space <: AbstractSpace
+        size::Int
+    end
+
+A structure that describes a single Hilbert space with size = `size`.
+"""
+struct Space <: AbstractSpace
+    size::Int
+
+    function Space(size::Int)
+        (size < 1) && throw(DomainError(size, "The size of Space must be positive integer (≥ 1)."))
+        return new(size)
+    end
+end
+
+dimensions_to_dims(s::Space) = SVector{1,Int}(s.size)
+
+# this creates a list of Space(1), it is used to generate `from` for Ket, and `to` for Bra)
+space_one_list(N::Int) = ntuple(i -> Space(1), Val(N))
+
+# TODO: introduce energy restricted space
+#=
+struct EnrSpace{N} <: AbstractSpace
+    size::Int
+    dims::SVector{N,Int}
+    n_excitations::Int
+    state2idx
+    idx2state
+end
+
+dimensions_to_dims(s::EnrSpace) = s.dims
+=#
diff --git a/src/qobj/states.jl b/src/qobj/states.jl
index 62c4b12c2..886982f7c 100644
--- a/src/qobj/states.jl
+++ b/src/qobj/states.jl
@@ -14,13 +14,13 @@ Returns a zero [`Ket`](@ref) vector with given argument `dimensions`.
 
 The `dimensions` can be either the following types:
 - `dimensions::Int`: Number of basis states in the Hilbert space.
-- `dimensions::Union{AbstractVector{Int}, Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+- `dimensions::Union{Dimensions,AbstractVector{Int}, Tuple}`: list of dimensions representing the each number of basis in the subsystems.
 
 !!! warning "Beware of type-stability!"
     It is highly recommended to use `zero_ket(dimensions)` with `dimensions` as `Tuple` or `SVector` to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 zero_ket(dimensions::Int) = QuantumObject(zeros(ComplexF64, dimensions), Ket, dimensions)
-zero_ket(dimensions::Union{AbstractVector{Int},Tuple}) =
+zero_ket(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}) =
     QuantumObject(zeros(ComplexF64, prod(dimensions)), Ket, dimensions)
 
 @doc raw"""
@@ -71,13 +71,13 @@ Generate a random normalized [`Ket`](@ref) vector with given argument `dimension
 
 The `dimensions` can be either the following types:
 - `dimensions::Int`: Number of basis states in the Hilbert space.
-- `dimensions::Union{AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+- `dimensions::Union{Dimensions,AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
 
 !!! warning "Beware of type-stability!"
     If you want to keep type stability, it is recommended to use `rand_ket(dimensions)` with `dimensions` as `Tuple` or `SVector` to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 rand_ket(dimensions::Int) = rand_ket(SVector(dimensions))
-function rand_ket(dimensions::Union{AbstractVector{Int},Tuple})
+function rand_ket(dimensions::Union{Dimensions,AbstractVector{Int},Tuple})
     N = prod(dimensions)
     ψ = rand(ComplexF64, N) .- (0.5 + 0.5im)
     return QuantumObject(normalize!(ψ); type = Ket, dims = dimensions)
@@ -144,13 +144,13 @@ Returns the maximally mixed density matrix with given argument `dimensions`.
 
 The `dimensions` can be either the following types:
 - `dimensions::Int`: Number of basis states in the Hilbert space.
-- `dimensions::Union{AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+- `dimensions::Union{Dimensions,AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
 
 !!! warning "Beware of type-stability!"
     If you want to keep type stability, it is recommended to use `maximally_mixed_dm(dimensions)` with `dimensions` as `Tuple` or `SVector` to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 maximally_mixed_dm(dimensions::Int) = QuantumObject(I(dimensions) / complex(dimensions), Operator, SVector(dimensions))
-function maximally_mixed_dm(dimensions::Union{AbstractVector{Int},Tuple})
+function maximally_mixed_dm(dimensions::Union{Dimensions,AbstractVector{Int},Tuple})
     N = prod(dimensions)
     return QuantumObject(I(N) / complex(N), Operator, dimensions)
 end
@@ -162,7 +162,7 @@ Generate a random density matrix from Ginibre ensemble with given argument `dime
 
 The `dimensions` can be either the following types:
 - `dimensions::Int`: Number of basis states in the Hilbert space.
-- `dimensions::Union{AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+- `dimensions::Union{Dimensions,AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
 
 The default keyword argument `rank = prod(dimensions)` (full rank).
 
@@ -174,7 +174,7 @@ The default keyword argument `rank = prod(dimensions)` (full rank).
 - [K. Życzkowski, et al., Generating random density matrices, Journal of Mathematical Physics 52, 062201 (2011)](http://dx.doi.org/10.1063/1.3595693)
 """
 rand_dm(dimensions::Int; rank::Int = prod(dimensions)) = rand_dm(SVector(dimensions), rank = rank)
-function rand_dm(dimensions::Union{AbstractVector{Int},Tuple}; rank::Int = prod(dimensions))
+function rand_dm(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}; rank::Int = prod(dimensions))
     N = prod(dimensions)
     (rank < 1) && throw(DomainError(rank, "The argument rank must be larger than 1."))
     (rank > N) && throw(DomainError(rank, "The argument rank cannot exceed dimensions."))
diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 010855122..90eb44ab8 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -77,7 +77,7 @@ The optional argument `Id_cache` can be used to pass a precomputed identity matr
 See also [`spost`](@ref) and [`sprepost`](@ref).
 """
 spre(A::AbstractQuantumObject{DT,OperatorQuantumObject}, Id_cache = I(size(A, 1))) where {DT} =
-    get_typename_wrapper(A)(_spre(A.data, Id_cache), SuperOperator, A.dims)
+    get_typename_wrapper(A)(_spre(A.data, Id_cache), SuperOperator, A.dimensions)
 
 @doc raw"""
     spost(B::AbstractQuantumObject, Id_cache=I(size(B,1)))
@@ -96,7 +96,7 @@ The optional argument `Id_cache` can be used to pass a precomputed identity matr
 See also [`spre`](@ref) and [`sprepost`](@ref).
 """
 spost(B::AbstractQuantumObject{DT,OperatorQuantumObject}, Id_cache = I(size(B, 1))) where {DT} =
-    get_typename_wrapper(B)(_spost(B.data, Id_cache), SuperOperator, B.dims)
+    get_typename_wrapper(B)(_spost(B.data, Id_cache), SuperOperator, B.dimensions)
 
 @doc raw"""
     sprepost(A::AbstractQuantumObject, B::AbstractQuantumObject)
@@ -116,8 +116,8 @@ function sprepost(
     A::AbstractQuantumObject{DT1,OperatorQuantumObject},
     B::AbstractQuantumObject{DT2,OperatorQuantumObject},
 ) where {DT1,DT2}
-    check_dims(A, B)
-    return promote_op_type(A, B)(_sprepost(A.data, B.data), SuperOperator, A.dims)
+    check_dimensions(A, B)
+    return promote_op_type(A, B)(_sprepost(A.data, B.data), SuperOperator, A.dimensions)
 end
 
 @doc raw"""
@@ -135,13 +135,13 @@ The optional argument `Id_cache` can be used to pass a precomputed identity matr
 See also [`spre`](@ref), [`spost`](@ref), and [`sprepost`](@ref).
 """
 lindblad_dissipator(O::AbstractQuantumObject{DT,OperatorQuantumObject}, Id_cache = I(size(O, 1))) where {DT} =
-    get_typename_wrapper(O)(_lindblad_dissipator(O.data, Id_cache), SuperOperator, O.dims)
+    get_typename_wrapper(O)(_lindblad_dissipator(O.data, Id_cache), SuperOperator, O.dimensions)
 
 # It is already a SuperOperator
 lindblad_dissipator(O::AbstractQuantumObject{DT,SuperOperatorQuantumObject}, Id_cache = nothing) where {DT} = O
 
 @doc raw"""
-    liouvillian(H::AbstractQuantumObject, c_ops::Union{Nothing,AbstractVector,Tuple}=nothing, Id_cache=I(prod(H.dims)))
+    liouvillian(H::AbstractQuantumObject, c_ops::Union{Nothing,AbstractVector,Tuple}=nothing, Id_cache=I(prod(H.dimensions)))
 
 Construct the Liouvillian [`SuperOperator`](@ref) for a system Hamiltonian ``\hat{H}`` and a set of collapse operators ``\{\hat{C}_n\}_n``:
 
@@ -162,7 +162,7 @@ See also [`spre`](@ref), [`spost`](@ref), and [`lindblad_dissipator`](@ref).
 function liouvillian(
     H::AbstractQuantumObject{DT,OpType},
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
-    Id_cache = I(prod(H.dims)),
+    Id_cache = I(prod(H.dimensions)),
 ) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
     L = liouvillian(H, Id_cache)
     if !(c_ops isa Nothing)
@@ -181,7 +181,7 @@ function liouvillian(
     return L
 end
 
-liouvillian(H::AbstractQuantumObject{DT,OperatorQuantumObject}, Id_cache::Diagonal = I(prod(H.dims))) where {DT} =
-    get_typename_wrapper(H)(_liouvillian(H.data, Id_cache), SuperOperator, H.dims)
+liouvillian(H::AbstractQuantumObject{DT,OperatorQuantumObject}, Id_cache::Diagonal = I(prod(H.dimensions))) where {DT} =
+    get_typename_wrapper(H)(_liouvillian(H.data, Id_cache), SuperOperator, H.dimensions)
 
 liouvillian(H::AbstractQuantumObject{DT,SuperOperatorQuantumObject}, Id_cache::Diagonal) where {DT} = H
diff --git a/src/spectrum.jl b/src/spectrum.jl
index a62a69ca6..af4d08b1a 100644
--- a/src/spectrum.jl
+++ b/src/spectrum.jl
@@ -84,8 +84,7 @@ function _spectrum(
     solver::ExponentialSeries;
     kwargs...,
 ) where {T1,T2,T3}
-    allequal((L.dims, A.dims, B.dims)) ||
-        throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
+    check_dimensions(L, A, B)
 
     rates, vecs, ρss = _spectrum_get_rates_vecs_ss(L, solver)
 
@@ -111,8 +110,7 @@ function _spectrum(
     solver::PseudoInverse;
     kwargs...,
 ) where {T1,T2,T3}
-    allequal((L.dims, A.dims, B.dims)) ||
-        throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
+    check_dimensions(L, A, B)
 
     ωList = convert(Vector{_FType(L)}, ωlist) # Convert it to support GPUs and avoid type instabilities
     Length = length(ωList)
@@ -123,7 +121,7 @@ function _spectrum(
     b = (spre(B) * ρss).data
 
     # multiply by operator A on the left (spre) and then perform trace operation
-    D = prod(L.dims)
+    D = prod(L.dimensions)
     _tr = SparseVector(D^2, [1 + n * (D + 1) for n in 0:(D-1)], ones(_CType(L), D)) # same as vec(system_identity_matrix)
     _tr_A = transpose(_tr) * spre(A).data
 
diff --git a/src/spin_lattice.jl b/src/spin_lattice.jl
index 9ce0094a3..99f356df3 100644
--- a/src/spin_lattice.jl
+++ b/src/spin_lattice.jl
@@ -51,7 +51,7 @@ function MultiSiteOperator(dims::Union{AbstractVector,Tuple}, pairs::Pair{<:Inte
 
     sites, ops = _get_unique_sites_ops(_sites, _ops)
 
-    collect(dims)[sites] == [op.dims[1] for op in ops] || throw(ArgumentError("The dimensions of the operators do not match the dimensions of the lattice."))
+    _dims[sites] == [get_dimensions_to(op)[1].size for op in ops] || throw(ArgumentError("The dimensions of the operators do not match the dimensions of the lattice."))
 
     data = kron(I(prod(_dims[1:sites[1]-1])), ops[1].data)
     for i in 2:length(sites)
diff --git a/src/steadystate.jl b/src/steadystate.jl
index 606b0907e..43b23a2c3 100644
--- a/src/steadystate.jl
+++ b/src/steadystate.jl
@@ -96,7 +96,7 @@ function _steadystate(
     kwargs...,
 ) where {T}
     L_tmp = L.data
-    N = prod(L.dims)
+    N = prod(L.dimensions)
     weight = norm(L_tmp, 1) / length(L_tmp)
 
     v0 = _get_dense_similar(L_tmp, N^2)
@@ -126,7 +126,7 @@ function _steadystate(
 
     ρss = reshape(ρss_vec, N, N)
     ρss = (ρss + ρss') / 2 # Hermitianize
-    return QuantumObject(ρss, Operator, L.dims)
+    return QuantumObject(ρss, Operator, L.dimensions)
 end
 
 function _steadystate(
@@ -134,7 +134,7 @@ function _steadystate(
     solver::SteadyStateEigenSolver;
     kwargs...,
 ) where {T}
-    N = prod(L.dims)
+    N = prod(L.dimensions)
 
     kwargs = merge((sigma = 1e-8, k = 1), (; kwargs...))
 
@@ -142,7 +142,7 @@ function _steadystate(
     ρss = reshape(ρss_vec, N, N)
     ρss /= tr(ρss)
     ρss = (ρss + ρss') / 2 # Hermitianize
-    return QuantumObject(ρss, Operator, L.dims)
+    return QuantumObject(ρss, Operator, L.dimensions)
 end
 
 function _steadystate(
@@ -150,7 +150,7 @@ function _steadystate(
     solver::SteadyStateDirectSolver,
 ) where {T}
     L_tmp = L.data
-    N = prod(L.dims)
+    N = prod(L.dimensions)
     weight = norm(L_tmp, 1) / length(L_tmp)
 
     v0 = _get_dense_similar(L_tmp, N^2)
@@ -170,7 +170,7 @@ function _steadystate(
     ρss_vec = L_tmp \ v0 # This is still not supported on GPU, yet
     ρss = reshape(ρss_vec, N, N)
     ρss = (ρss + ρss') / 2 # Hermitianize
-    return QuantumObject(ρss, Operator, L.dims)
+    return QuantumObject(ρss, Operator, L.dimensions)
 end
 
 _steadystate(
@@ -411,13 +411,13 @@ function _steadystate_floquet(
     ρ0 = reshape(ρtot[offset1+1:offset2], Ns, Ns)
     ρ0_tr = tr(ρ0)
     ρ0 = ρ0 / ρ0_tr
-    ρ0 = QuantumObject((ρ0 + ρ0') / 2, type = Operator, dims = L_0.dims)
+    ρ0 = QuantumObject((ρ0 + ρ0') / 2, type = Operator, dims = L_0.dimensions)
     ρtot = ρtot / ρ0_tr
 
     ρ_list = [ρ0]
     for i in 0:n_max-1
         ρi_m = reshape(ρtot[offset1-(i+1)*N+1:offset1-i*N], Ns, Ns)
-        ρi_m = QuantumObject(ρi_m, type = Operator, dims = L_0.dims)
+        ρi_m = QuantumObject(ρi_m, type = Operator, dims = L_0.dimensions)
         push!(ρ_list, ρi_m)
     end
 
@@ -434,8 +434,7 @@ function _steadystate_floquet(
     tol::R = 1e-8,
     kwargs...,
 ) where {R<:Real}
-    ((L_0.dims == L_p.dims) && (L_0.dims == L_m.dims)) ||
-        throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
+    check_dimensions(L_0, L_p, L_m)
 
     L_eff = liouvillian_floquet(L_0, L_p, L_m, ωd; n_max = n_max, tol = tol)
 
diff --git a/src/time_evolution/lr_mesolve.jl b/src/time_evolution/lr_mesolve.jl
index 8b71e94d9..f8a62cb8f 100644
--- a/src/time_evolution/lr_mesolve.jl
+++ b/src/time_evolution/lr_mesolve.jl
@@ -401,7 +401,7 @@ function lr_mesolveProblem(
     opt::NamedTuple = lr_mesolve_options_default,
     kwargs...,
 ) where {DT1,T2}
-    Hdims = H.dims
+    Hdims = H.dimensions
 
     # Formulation of problem
     H -= 0.5im * mapreduce(op -> op' * op, +, c_ops)
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index db8ab0319..99b5d3fad 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -299,7 +299,7 @@ function mcsolveEnsembleProblem(
     ensemble_prob = TimeEvolutionProblem(
         EnsembleProblem(prob_mc.prob, prob_func = _prob_func, output_func = _output_func[1], safetycopy = false),
         prob_mc.times,
-        prob_mc.dims,
+        prob_mc.dimensions,
         (progr = _output_func[2], channel = _output_func[3]),
     )
 
@@ -472,7 +472,7 @@ function mcsolve(
 )
     sol = _mcsolve_solve_ens(ens_prob_mc, alg, ensemble_method, ntraj)
 
-    dims = ens_prob_mc.dims
+    dims = ens_prob_mc.dimensions
     _sol_1 = sol[:, 1]
     _expvals_sol_1 = _mcsolve_get_expvals(_sol_1)
 
diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
index 2d9f42efc..d119785b3 100644
--- a/src/time_evolution/mesolve.jl
+++ b/src/time_evolution/mesolve.jl
@@ -73,7 +73,7 @@ function mesolveProblem(
     tlist = convert(Vector{_FType(ψ0)}, tlist) # Convert it to support GPUs and avoid type instabilities for OrdinaryDiffEq.jl
 
     L_evo = _mesolve_make_L_QobjEvo(H, c_ops)
-    check_dims(L_evo, ψ0)
+    check_dimensions(L_evo, ψ0)
 
     T = Base.promote_eltype(L_evo, ψ0)
     ρ0 = sparse_to_dense(_CType(T), mat2vec(ket2dm(ψ0).data)) # Convert it to dense vector with complex element type
@@ -89,7 +89,7 @@ function mesolveProblem(
     tspan = (tlist[1], tlist[end])
     prob = ODEProblem{getVal(inplace),FullSpecialize}(L, ρ0, tspan, params; kwargs3...)
 
-    return TimeEvolutionProblem(prob, tlist, L_evo.dims)
+    return TimeEvolutionProblem(prob, tlist, L_evo.dimensions)
 end
 
 @doc raw"""
@@ -179,7 +179,7 @@ end
 function mesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit5())
     sol = solve(prob.prob, alg)
 
-    ρt = map(ϕ -> QuantumObject(vec2mat(ϕ), type = Operator, dims = prob.dims), sol.u)
+    ρt = map(ϕ -> QuantumObject(vec2mat(ϕ), type = Operator, dims = prob.dimensions), sol.u)
 
     return TimeEvolutionSol(
         prob.times,
diff --git a/src/time_evolution/sesolve.jl b/src/time_evolution/sesolve.jl
index 28562e280..341a2cc55 100644
--- a/src/time_evolution/sesolve.jl
+++ b/src/time_evolution/sesolve.jl
@@ -1,8 +1,8 @@
 export sesolveProblem, sesolve
 
 _sesolve_make_U_QobjEvo(H::QuantumObjectEvolution{<:MatrixOperator}) =
-    QobjEvo(MatrixOperator(-1im * H.data.A), dims = H.dims, type = Operator)
-_sesolve_make_U_QobjEvo(H::QuantumObject) = QobjEvo(MatrixOperator(-1im * H.data), dims = H.dims, type = Operator)
+    QobjEvo(MatrixOperator(-1im * H.data.A), dims = H.dimensions, type = Operator)
+_sesolve_make_U_QobjEvo(H::QuantumObject) = QobjEvo(MatrixOperator(-1im * H.data), dims = H.dimensions, type = Operator)
 _sesolve_make_U_QobjEvo(H::Union{QuantumObjectEvolution,Tuple}) = QobjEvo(H, -1im)
 
 @doc raw"""
@@ -62,7 +62,7 @@ function sesolveProblem(
 
     H_evo = _sesolve_make_U_QobjEvo(H) # Multiply by -i
     isoper(H_evo) || throw(ArgumentError("The Hamiltonian must be an Operator."))
-    check_dims(H_evo, ψ0)
+    check_dimensions(H_evo, ψ0)
 
     T = Base.promote_eltype(H_evo, ψ0)
     ψ0 = sparse_to_dense(_CType(T), get_data(ψ0)) # Convert it to dense vector with complex element type
@@ -78,7 +78,7 @@ function sesolveProblem(
     tspan = (tlist[1], tlist[end])
     prob = ODEProblem{getVal(inplace),FullSpecialize}(U, ψ0, tspan, params; kwargs3...)
 
-    return TimeEvolutionProblem(prob, tlist, H_evo.dims)
+    return TimeEvolutionProblem(prob, tlist, H_evo.dimensions)
 end
 
 @doc raw"""
@@ -152,7 +152,7 @@ end
 function sesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit5())
     sol = solve(prob.prob, alg)
 
-    ψt = map(ϕ -> QuantumObject(ϕ, type = Ket, dims = prob.dims), sol.u)
+    ψt = map(ϕ -> QuantumObject(ϕ, type = Ket, dims = prob.dimensions), sol.u)
 
     return TimeEvolutionSol(
         prob.times,
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index 2dc7fac16..cdd2de9bb 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -166,8 +166,8 @@ function ssesolveProblem(
 
     H_eff_evo = _mcsolve_make_Heff_QobjEvo(H, sc_ops)
     isoper(H_eff_evo) || throw(ArgumentError("The Hamiltonian must be an Operator."))
-    check_dims(H_eff_evo, ψ0)
-    dims = H_eff_evo.dims
+    check_dimensions(H_eff_evo, ψ0)
+    dims = H_eff_evo.dimensions
 
     ψ0 = sparse_to_dense(_CType(ψ0), get_data(ψ0))
 
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index b0f701f3f..ec429b1d1 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -15,16 +15,28 @@ A Julia constructor for handling the `ODEProblem` of the time evolution of quant
 
 - `prob::AbstractSciMLProblem`: The `ODEProblem` of the time evolution.
 - `times::Abstractvector`: The time list of the evolution.
-- `dims::Abstractvector`: The dimensions of the Hilbert space.
+- `dimensions::AbstractDimensions`: The dimensions of the Hilbert space.
 - `kwargs::KWT`: Generic keyword arguments.
+
+!!! note "`dims` property"
+    For a given `prob::TimeEvolutionProblem`, `prob.dims` or `getproperty(prob, :dims)` returns its `dimensions` in the type of integer-vector.
 """
-struct TimeEvolutionProblem{PT<:AbstractSciMLProblem,TT<:AbstractVector,DT<:AbstractVector,KWT}
+struct TimeEvolutionProblem{PT<:AbstractSciMLProblem,TT<:AbstractVector,DT<:AbstractDimensions,KWT}
     prob::PT
     times::TT
-    dims::DT
+    dimensions::DT
     kwargs::KWT
 end
 
+function Base.getproperty(prob::TimeEvolutionProblem, key::Symbol)
+    # a comment here to avoid bad render by JuliaFormatter
+    if key === :dims
+        return dimensions_to_dims(getfield(prob, :dimensions))
+    else
+        return getfield(prob, key)
+    end
+end
+
 TimeEvolutionProblem(prob, times, dims) = TimeEvolutionProblem(prob, times, dims, nothing)
 
 @doc raw"""
@@ -210,9 +222,7 @@ function liouvillian_floquet(
     n_max::Int = 3,
     tol::Real = 1e-15,
 ) where {T1,T2,T3}
-    ((L₀.dims == Lₚ.dims) && (L₀.dims == Lₘ.dims)) ||
-        throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
-
+    check_dimensions(L₀, Lₚ, Lₘ)
     return _liouvillian_floquet(L₀, Lₚ, Lₘ, ω, n_max, tol)
 end
 
@@ -253,11 +263,11 @@ function liouvillian_generalized(
 ) where {MT<:AbstractMatrix}
     (length(fields) == length(T_list)) || throw(DimensionMismatch("The number of fields, ωs and Ts must be the same."))
 
-    dims = (N_trunc isa Nothing) ? H.dims : SVector(N_trunc)
+    dims = (N_trunc isa Nothing) ? H.dimensions : SVector(N_trunc)
     final_size = prod(dims)
     result = eigen(H)
     E = real.(result.values[1:final_size])
-    U = QuantumObject(result.vectors, result.type, result.dims)
+    U = QuantumObject(result.vectors, result.type, result.dimensions)
 
     H_d = QuantumObject(Diagonal(complex(E)), type = Operator, dims = dims)
 
@@ -328,6 +338,6 @@ function _liouvillian_floquet(
         T = -(L_0 + 1im * n_i * ω * I + L_p * T) \ L_m_dense
     end
 
-    tol == 0 && return QuantumObject(L_0 + L_m * S + L_p * T, SuperOperator, L₀.dims)
-    return QuantumObject(dense_to_sparse(L_0 + L_m * S + L_p * T, tol), SuperOperator, L₀.dims)
+    tol == 0 && return QuantumObject(L_0 + L_m * S + L_p * T, SuperOperator, L₀.dimensions)
+    return QuantumObject(dense_to_sparse(L_0 + L_m * S + L_p * T, tol), SuperOperator, L₀.dimensions)
 end
diff --git a/src/time_evolution/time_evolution_dynamical.jl b/src/time_evolution/time_evolution_dynamical.jl
index 650f406e0..81cf3260d 100644
--- a/src/time_evolution/time_evolution_dynamical.jl
+++ b/src/time_evolution/time_evolution_dynamical.jl
@@ -149,7 +149,7 @@ function dfd_mesolveProblem(
     tol_list::Vector{<:Number} = fill(1e-8, length(maxdims)),
     kwargs...,
 ) where {DT1,T2<:Integer,StateOpType<:Union{KetQuantumObject,OperatorQuantumObject}}
-    length(ψ0.dims) != length(maxdims) &&
+    length(ψ0.dimensions) != length(maxdims) &&
         throw(DimensionMismatch("'dim_list' and 'maxdims' do not have the same dimension."))
 
     dim_list = MVector(ψ0.dims)
@@ -362,7 +362,7 @@ function dsf_mesolveProblem(
     op_l_vec = map(op -> mat2vec(get_data(op)'), op_list)
     # Create the Krylov subspace with kron(H₀.data, H₀.data) just for initialize
     expv_cache = arnoldi(kron(H₀.data, H₀.data), mat2vec(ket2dm(ψ0).data), krylov_dim)
-    dsf_identity = I(prod(H₀.dims))
+    dsf_identity = I(prod(H₀.dimensions))
     dsf_displace_cache_left = sum(op -> ScalarOperator(one(T)) * MatrixOperator(kron(op.data, dsf_identity)), op_list)
     dsf_displace_cache_left_dag =
         sum(op -> ScalarOperator(one(T)) * MatrixOperator(kron(sparse(op.data'), dsf_identity)), op_list)
diff --git a/test/core-test/low_rank_dynamics.jl b/test/core-test/low_rank_dynamics.jl
index ef767f304..fc1d0769c 100644
--- a/test/core-test/low_rank_dynamics.jl
+++ b/test/core-test/low_rank_dynamics.jl
@@ -9,7 +9,7 @@
     M = latt.N + 1       # Number of states in the LR basis
 
     # Define initial state
-    ϕ = Vector{QuantumObject{Vector{ComplexF64},KetQuantumObject,M - 1}}(undef, M)
+    ϕ = Vector{QuantumObject{Vector{ComplexF64},KetQuantumObject,Dimensions{M - 1,NTuple{M - 1,Space}}}}(undef, M)
     ϕ[1] = kron(fill(basis(2, 1), N_modes)...)
 
     i = 1
diff --git a/test/core-test/negativity_and_partial_transpose.jl b/test/core-test/negativity_and_partial_transpose.jl
index ba10db898..ea1c1580e 100644
--- a/test/core-test/negativity_and_partial_transpose.jl
+++ b/test/core-test/negativity_and_partial_transpose.jl
@@ -35,6 +35,7 @@
             end
         end
         @test_throws ArgumentError partial_transpose(A_dense, [true])
+        @test_throws ArgumentError partial_transpose(Qobj(zeros(ComplexF64, 3, 2)), [true]) # invalid GeneralDimensions
 
         @testset "Type Inference (partial_transpose)" begin
             @inferred partial_transpose(A_dense, [true, false, true])
diff --git a/test/core-test/quantum_objects.jl b/test/core-test/quantum_objects.jl
index 9ad54a760..91738bdc9 100644
--- a/test/core-test/quantum_objects.jl
+++ b/test/core-test/quantum_objects.jl
@@ -13,16 +13,27 @@
     # DomainError: incompatible between size of array and type
     @testset "DomainError" begin
         a = rand(ComplexF64, 3, 2)
-        for t in [nothing, Operator, SuperOperator, Bra, OperatorBra]
+        for t in [SuperOperator, Bra, OperatorBra]
             @test_throws DomainError Qobj(a, type = t)
         end
+
         a = rand(ComplexF64, 2, 2, 2)
         for t in [nothing, Ket, Bra, Operator, SuperOperator, OperatorBra, OperatorKet]
             @test_throws DomainError Qobj(a, type = t)
         end
+
         a = rand(ComplexF64, 1, 2)
         @test_throws DomainError Qobj(a, type = Operator)
         @test_throws DomainError Qobj(a, type = SuperOperator)
+
+        @test_throws DomainError Qobj(rand(ComplexF64, 2, 1), type = Operator) # should be type = Bra
+
+        # check that Ket, Bra, SuperOperator, OperatorKet, and OperatorBra don't support GeneralDimensions
+        @test_throws DomainError Qobj(rand(ComplexF64, 2), type = Ket, dims = ((2,), (1,)))
+        @test_throws DomainError Qobj(rand(ComplexF64, 1, 2), type = Bra, dims = ((1,), (2,)))
+        @test_throws DomainError Qobj(rand(ComplexF64, 4, 4), type = SuperOperator, dims = ((2,), (2,)))
+        @test_throws DomainError Qobj(rand(ComplexF64, 4), type = OperatorKet, dims = ((2,), (1,)))
+        @test_throws DomainError Qobj(rand(ComplexF64, 1, 4), type = OperatorBra, dims = ((1,), (2,)))
     end
 
     # unsupported type of dims
@@ -74,6 +85,7 @@
         a = sprand(ComplexF64, 100, 100, 0.1)
         a2 = Qobj(a)
         a3 = Qobj(a, type = SuperOperator)
+        a4 = Qobj(sprand(ComplexF64, 100, 10, 0.1)) # GeneralDimensions
         @test isket(a2) == false
         @test isbra(a2) == false
         @test isoper(a2) == true
@@ -83,6 +95,7 @@
         @test iscached(a2) == true
         @test isconstant(a2) == true
         @test isunitary(a2) == false
+        @test a2.dims == [100]
         @test isket(a3) == false
         @test isbra(a3) == false
         @test isoper(a3) == false
@@ -92,7 +105,20 @@
         @test iscached(a3) == true
         @test isconstant(a3) == true
         @test isunitary(a3) == false
+        @test a3.dims == [10]
+        @test isket(a4) == false
+        @test isbra(a4) == false
+        @test isoper(a4) == true
+        @test issuper(a4) == false
+        @test isoperket(a4) == false
+        @test isoperbra(a4) == false
+        @test iscached(a4) == true
+        @test isconstant(a4) == true
+        @test isunitary(a4) == false
+        @test a4.dims == [[100], [10]]
         @test_throws DimensionMismatch Qobj(a, dims = 2)
+        @test_throws DimensionMismatch Qobj(a4.data, dims = 2)
+        @test_throws DimensionMismatch Qobj(a4.data, dims = ((100,), (2,)))
     end
 
     @testset "OperatorKet and OperatorBra" begin
@@ -235,6 +261,16 @@
         @test opstring ==
               "\nQuantum Object:   type=Operator   dims=$a_dims   size=$a_size   ishermitian=$a_isherm\n$datastring"
 
+        # GeneralDimensions
+        Gop = tensor(a, ψ)
+        opstring = sprint((t, s) -> show(t, "text/plain", s), Gop)
+        datastring = sprint((t, s) -> show(t, "text/plain", s), Gop.data)
+        Gop_dims = [[N, N], [N, 1]]
+        Gop_size = size(Gop)
+        Gop_isherm = isherm(Gop)
+        @test opstring ==
+              "\nQuantum Object:   type=Operator   dims=$Gop_dims   size=$Gop_size   ishermitian=$Gop_isherm\n$datastring"
+
         a = spre(a)
         opstring = sprint((t, s) -> show(t, "text/plain", s), a)
         datastring = sprint((t, s) -> show(t, "text/plain", s), a.data)
@@ -310,24 +346,29 @@
         for T in [ComplexF32, ComplexF64]
             N = 4
             a = rand(T, N)
-            @inferred QuantumObject{typeof(a),KetQuantumObject} Qobj(a)
+            @inferred QuantumObject{typeof(a),KetQuantumObject,Dimensions{1}} Qobj(a)
             for type in [Ket, OperatorKet]
                 @inferred Qobj(a, type = type)
             end
 
-            UnionType =
-                Union{QuantumObject{Matrix{T},BraQuantumObject,1},QuantumObject{Matrix{T},OperatorQuantumObject,1}}
+            UnionType = Union{
+                QuantumObject{Matrix{T},BraQuantumObject,Dimensions{1,Tuple{Space}}},
+                QuantumObject{Matrix{T},OperatorQuantumObject,Dimensions{1,Tuple{Space}}},
+            }
             a = rand(T, 1, N)
             @inferred UnionType Qobj(a)
             for type in [Bra, OperatorBra]
                 @inferred Qobj(a, type = type)
             end
 
+            UnionType2 = Union{
+                QuantumObject{Matrix{T},OperatorQuantumObject,GeneralDimensions{1,Tuple{Space},Tuple{Space}}},
+                QuantumObject{Matrix{T},OperatorQuantumObject,Dimensions{1,Tuple{Space}}},
+            }
             a = rand(T, N, N)
             @inferred UnionType Qobj(a)
-            for type in [Operator, SuperOperator]
-                @inferred Qobj(a, type = type)
-            end
+            @inferred UnionType2 Qobj(a, type = Operator)
+            @inferred Qobj(a, type = SuperOperator)
         end
 
         @testset "Math Operation" begin
@@ -629,8 +670,26 @@
         ρ = kron(ρ1, ρ2)
         ρ1_ptr = ptrace(ρ, 1)
         ρ2_ptr = ptrace(ρ, 2)
-        @test ρ1.data ≈ ρ1_ptr.data atol = 1e-10
-        @test ρ2.data ≈ ρ2_ptr.data atol = 1e-10
+
+        # use GeneralDimensions to do partial trace
+        ρ1_compound = Qobj(zeros(ComplexF64, 2, 2), dims = ((2, 1), (2, 1)))
+        II = qeye(2)
+        basis_list = [basis(2, i) for i in 0:1]
+        for b in basis_list
+            ρ1_compound += tensor(II, b') * ρ * tensor(II, b)
+        end
+        ρ2_compound = Qobj(zeros(ComplexF64, 2, 2), dims = ((1, 2), (1, 2)))
+        for b in basis_list
+            ρ2_compound += tensor(b', II) * ρ * tensor(b, II)
+        end
+        @test ρ1.data ≈ ρ1_ptr.data ≈ ρ1_compound.data
+        @test ρ2.data ≈ ρ2_ptr.data ≈ ρ2_compound.data
+        @test ρ1.dims != ρ1_compound.dims
+        @test ρ2.dims != ρ2_compound.dims
+        ρ1_compound = ptrace(ρ1_compound, 1)
+        ρ2_compound = ptrace(ρ2_compound, 2)
+        @test ρ1.dims == ρ1_compound.dims
+        @test ρ2.dims == ρ2_compound.dims
 
         ψlist = [rand_ket(2), rand_ket(3), rand_ket(4), rand_ket(5)]
         ρlist = [rand_dm(2), rand_dm(3), rand_dm(4), rand_dm(5)]
@@ -683,6 +742,7 @@
         @test_throws ArgumentError ptrace(ρtotal, (0, 2))
         @test_throws ArgumentError ptrace(ρtotal, (2, 5))
         @test_throws ArgumentError ptrace(ρtotal, (2, 2, 3))
+        @test_throws ArgumentError ptrace(Qobj(zeros(ComplexF64, 3, 2)), 1) # invalid GeneralDimensions
 
         @testset "Type Inference (ptrace)" begin
             @inferred ptrace(ρ, 1)
@@ -695,6 +755,7 @@
     end
 
     @testset "permute" begin
+        # standard Dimensions
         ket_a = Qobj(rand(ComplexF64, 2))
         ket_b = Qobj(rand(ComplexF64, 3))
         ket_c = Qobj(rand(ComplexF64, 4))
@@ -729,10 +790,22 @@
         @test_throws ArgumentError permute(op_bdca, wrong_order1)
         @test_throws ArgumentError permute(op_bdca, wrong_order2)
 
+        # GeneralDimensions
+        Gop_d = Qobj(rand(ComplexF64, 5, 6))
+        compound_bdca = permute(tensor(ket_a, op_b, bra_c, Gop_d), (2, 4, 3, 1))
+        compound_dacb = permute(tensor(ket_a, op_b, bra_c, Gop_d), (4, 1, 3, 2))
+        @test compound_bdca ≈ tensor(op_b, Gop_d, bra_c, ket_a)
+        @test compound_dacb ≈ tensor(Gop_d, ket_a, bra_c, op_b)
+        @test compound_bdca.dims == [[3, 5, 1, 2], [3, 6, 4, 1]]
+        @test compound_dacb.dims == [[5, 2, 1, 3], [6, 1, 4, 3]]
+        @test isoper(compound_bdca)
+        @test isoper(compound_dacb)
+
         @testset "Type Inference (permute)" begin
             @inferred permute(ket_bdca, (2, 4, 3, 1))
             @inferred permute(bra_bdca, (2, 4, 3, 1))
             @inferred permute(op_bdca, (2, 4, 3, 1))
+            @inferred permute(compound_bdca, (2, 4, 3, 1))
         end
     end
 end
diff --git a/test/core-test/quantum_objects_evo.jl b/test/core-test/quantum_objects_evo.jl
index f5b39c557..c64004e5b 100644
--- a/test/core-test/quantum_objects_evo.jl
+++ b/test/core-test/quantum_objects_evo.jl
@@ -2,9 +2,10 @@
     # DomainError: incompatible between size of array and type
     @testset "Thrown Errors" begin
         a = MatrixOperator(rand(ComplexF64, 3, 2))
-        for t in [Operator, SuperOperator]
-            @test_throws DomainError QobjEvo(a, type = t)
-        end
+        @test_throws DomainError QobjEvo(a, type = SuperOperator)
+
+        a = MatrixOperator(rand(ComplexF64, 4, 4))
+        @test_throws DomainError QobjEvo(a, type = SuperOperator, dims = ((2,), (2,)))
 
         a = MatrixOperator(rand(ComplexF64, 3, 2))
         for t in (Ket, Bra, OperatorKet, OperatorBra)
@@ -131,10 +132,13 @@
         N = 4
         for T in [ComplexF32, ComplexF64]
             a = MatrixOperator(rand(T, N, N))
-            @inferred QobjEvo(a)
-            for type in [Operator, SuperOperator]
-                @inferred QobjEvo(a, type = type)
-            end
+            UnionType = Union{
+                QuantumObjectEvolution{typeof(a),OperatorQuantumObject,GeneralDimensions{1,Tuple{Space},Tuple{Space}}},
+                QuantumObjectEvolution{typeof(a),OperatorQuantumObject,Dimensions{1,Tuple{Space}}},
+            }
+            @inferred UnionType QobjEvo(a)
+            @inferred UnionType QobjEvo(a, type = Operator)
+            @inferred QobjEvo(a, type = SuperOperator)
         end
 
         a = destroy(N)
diff --git a/test/core-test/steady_state.jl b/test/core-test/steady_state.jl
index 507f1227a..630529559 100644
--- a/test/core-test/steady_state.jl
+++ b/test/core-test/steady_state.jl
@@ -46,10 +46,6 @@
         @inferred steadystate(H, c_ops, solver = solver)
         @inferred steadystate(L, solver = solver)
 
-        solver = SteadyStateLinearSolver()
-        @inferred steadystate(H, c_ops, solver = solver)
-        @inferred steadystate(L, solver = solver)
-
         solver = SteadyStateEigenSolver()
         @inferred steadystate(H, c_ops, solver = solver)
         @inferred steadystate(L, solver = solver)