diff --git a/docs/src/getting_started/type_stability.md b/docs/src/getting_started/type_stability.md
index bd22d3481..6b3a13e8a 100644
--- a/docs/src/getting_started/type_stability.md
+++ b/docs/src/getting_started/type_stability.md
@@ -183,7 +183,7 @@ The `dims` field contains the dimensions of the subsystems (in this case, three
 
 ```@example type-stability
 function reshape_operator_data(dims)
-    op = Qobj(randn(prod(dims), prod(dims)), type=Operator(), dims=dims)
+    op = Qobj(randn(hilbert_dimensions_to_size(dims)...), type=Operator(), dims=dims)
     op_dims = op.dims
     op_data = op.data
     return reshape(op_data, vcat(op_dims, op_dims)...)
diff --git a/docs/src/resources/api.md b/docs/src/resources/api.md
index 1bdcb9448..272f66490 100644
--- a/docs/src/resources/api.md
+++ b/docs/src/resources/api.md
@@ -17,10 +17,10 @@ end
 ## [Quantum object (Qobj) and type](@id doc-API:Quantum-object-and-type)
 
 ```@docs
-Space
+HilbertSpace
 EnrSpace
-Dimensions
-GeneralDimensions
+ProductDimensions
+GeneralProductDimensions
 AbstractQuantumObject
 Bra
 Ket
@@ -326,6 +326,8 @@ convert_unit
 row_major_reshape
 meshgrid
 enr_state_dictionaries
+hilbert_dimensions_to_size
+liouville_dimensions_to_size
 ```
 
 ## [Visualization](@id doc-API:Visualization)
diff --git a/ext/QuantumToolboxMakieExt.jl b/ext/QuantumToolboxMakieExt.jl
index 015acae18..2c7f0c716 100644
--- a/ext/QuantumToolboxMakieExt.jl
+++ b/ext/QuantumToolboxMakieExt.jl
@@ -226,7 +226,7 @@ function _plot_fock_distribution(
     kwargs...,
 ) where {SType<:Union{Bra,Ket,Operator}}
     ρ = ket2dm(ρ)
-    D = prod(ρ.dims)
+    D = hilbert_dimensions_to_size(ρ.dims)[1]
     isapprox(tr(ρ), 1, atol = 1e-4) || (@warn "The input ρ should be normalized.")
 
     xvec = 0:(D-1)
diff --git a/src/correlations.jl b/src/correlations.jl
index a8ac92f27..9dcd94c11 100644
--- a/src/correlations.jl
+++ b/src/correlations.jl
@@ -113,7 +113,7 @@ function correlation_2op_2t(
     reverse::Bool = false,
     kwargs...,
 ) where {HOpType<:Union{Operator,SuperOperator},StateOpType<:Union{Ket,Operator}}
-    C = eye(prod(H.dimensions), dims = H.dimensions)
+    C = eye(hilbert_dimensions_to_size(H.dimensions)[1], dims = H.dimensions)
     if reverse
         corr = correlation_3op_2t(H, ψ0, tlist, τlist, c_ops, A, B, C; kwargs...)
     else
diff --git a/src/entropy.jl b/src/entropy.jl
index b4e4bf91c..54d2754d8 100644
--- a/src/entropy.jl
+++ b/src/entropy.jl
@@ -221,7 +221,7 @@ Calculate the [concurrence](https://en.wikipedia.org/wiki/Concurrence_(quantum_c
 - [Hill-Wootters1997](@citet)
 """
 function concurrence(ρ::QuantumObject{OpType}) where {OpType<:Union{Ket,Operator}}
-    (ρ.dimensions == Dimensions((Space(2), Space(2)))) || throw(
+    (ρ.dimensions == ProductDimensions((HilbertSpace(2), HilbertSpace(2)))) || throw(
         ArgumentError(
             "The `concurrence` only works for a two-qubit state, invalid dims = $(_get_dims_string(ρ.dimensions)).",
         ),
diff --git a/src/negativity.jl b/src/negativity.jl
index 286ce37a0..23b481ce0 100644
--- a/src/negativity.jl
+++ b/src/negativity.jl
@@ -76,7 +76,7 @@ function partial_transpose(ρ::QuantumObject{Operator}, mask::Vector{Bool})
     (length(mask) != length(ρ.dimensions)) &&
         throw(ArgumentError("The length of \`mask\` should be equal to the length of \`ρ.dims\`."))
 
-    isa(ρ.dimensions, GeneralDimensions) &&
+    isa(ρ.dimensions, GeneralProductDimensions) &&
         (get_dimensions_to(ρ) != get_dimensions_from(ρ)) &&
         throw(ArgumentError("Invalid partial transpose for dims = $(_get_dims_string(ρ.dimensions))"))
 
@@ -100,7 +100,7 @@ function _partial_transpose(ρ::QuantumObject{Operator}, mask::Vector{Bool})
     return QuantumObject(
         reshape(permutedims(reshape(ρ.data, (dims_rev..., dims_rev...)), pt_idx), size(ρ)),
         Operator(),
-        Dimensions(ρ.dimensions.to),
+        ProductDimensions(ρ.dimensions.to),
     )
 end
 
diff --git a/src/qobj/arithmetic_and_attributes.jl b/src/qobj/arithmetic_and_attributes.jl
index 14e06aae4..ac97d74e2 100644
--- a/src/qobj/arithmetic_and_attributes.jl
+++ b/src/qobj/arithmetic_and_attributes.jl
@@ -59,9 +59,9 @@ function check_mul_dimensions(from::NTuple{NA,AbstractSpace}, to::NTuple{NB,Abst
     return nothing
 end
 
-for ADimType in (:Dimensions, :GeneralDimensions)
-    for BDimType in (:Dimensions, :GeneralDimensions)
-        if ADimType == BDimType == :Dimensions
+for ADimType in (:ProductDimensions, :GeneralProductDimensions)
+    for BDimType in (:ProductDimensions, :GeneralProductDimensions)
+        if ADimType == BDimType == :ProductDimensions
             @eval begin
                 function Base.:(*)(
                     A::AbstractQuantumObject{Operator,<:$ADimType},
@@ -83,7 +83,7 @@ for ADimType in (:Dimensions, :GeneralDimensions)
                     return QType(
                         A.data * B.data,
                         Operator(),
-                        GeneralDimensions(get_dimensions_to(A), get_dimensions_from(B)),
+                        GeneralProductDimensions(get_dimensions_to(A), get_dimensions_from(B)),
                     )
                 end
             end
@@ -91,13 +91,13 @@ for ADimType in (:Dimensions, :GeneralDimensions)
     end
 end
 
-function Base.:(*)(A::AbstractQuantumObject{Operator}, B::QuantumObject{Ket,<:Dimensions})
+function Base.:(*)(A::AbstractQuantumObject{Operator}, B::QuantumObject{Ket,<:ProductDimensions})
     check_mul_dimensions(get_dimensions_from(A), get_dimensions_to(B))
-    return QuantumObject(A.data * B.data, Ket(), Dimensions(get_dimensions_to(A)))
+    return QuantumObject(A.data * B.data, Ket(), ProductDimensions(get_dimensions_to(A)))
 end
-function Base.:(*)(A::QuantumObject{Bra,<:Dimensions}, B::AbstractQuantumObject{Operator})
+function Base.:(*)(A::QuantumObject{Bra,<:ProductDimensions}, B::AbstractQuantumObject{Operator})
     check_mul_dimensions(get_dimensions_from(A), get_dimensions_to(B))
-    return QuantumObject(A.data * B.data, Bra(), Dimensions(get_dimensions_from(B)))
+    return QuantumObject(A.data * B.data, Bra(), ProductDimensions(get_dimensions_from(B)))
 end
 function Base.:(*)(A::QuantumObject{Ket}, B::QuantumObject{Bra})
     check_dimensions(A, B)
@@ -523,7 +523,7 @@ function ptrace(QO::QuantumObject{Ket}, sel::Union{AbstractVector{Int},Tuple})
 
     _sort_sel = sort(SVector{length(sel),Int}(sel))
     ρtr, dkeep = _ptrace_ket(QO.data, QO.dims, _sort_sel)
-    return QuantumObject(ρtr, type = Operator(), dims = Dimensions(dkeep))
+    return QuantumObject(ρtr, type = Operator(), dims = ProductDimensions(dkeep))
 end
 
 ptrace(QO::QuantumObject{Bra}, sel::Union{AbstractVector{Int},Tuple}) = ptrace(QO', sel)
@@ -532,7 +532,7 @@ function ptrace(QO::QuantumObject{Operator}, sel::Union{AbstractVector{Int},Tupl
     any(s -> s isa EnrSpace, QO.dimensions.to) && throw(ArgumentError("ptrace does not support EnrSpace"))
 
     # TODO: support for special cases when some of the subsystems have same `to` and `from` space
-    isa(QO.dimensions, GeneralDimensions) &&
+    isa(QO.dimensions, GeneralProductDimensions) &&
         (get_dimensions_to(QO) != get_dimensions_from(QO)) &&
         throw(ArgumentError("Invalid partial trace for dims = $(_get_dims_string(QO.dimensions))"))
 
@@ -553,7 +553,7 @@ function ptrace(QO::QuantumObject{Operator}, sel::Union{AbstractVector{Int},Tupl
     dims = dimensions_to_dims(get_dimensions_to(QO))
     _sort_sel = sort(SVector{length(sel),Int}(sel))
     ρtr, dkeep = _ptrace_oper(QO.data, dims, _sort_sel)
-    return QuantumObject(ρtr, type = Operator(), dims = Dimensions(dkeep))
+    return QuantumObject(ρtr, type = Operator(), dims = ProductDimensions(dkeep))
 end
 ptrace(QO::QuantumObject, sel::Int) = ptrace(QO, SVector(sel))
 
@@ -672,7 +672,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{Ket})
-    a = destroy(prod(ψ.dimensions))
+    a = destroy(hilbert_dimensions_to_size(ψ.dimensions)[1])
     α = expect(a, ψ)
     D = exp(α * a' - conj(α) * a)
 
@@ -680,7 +680,7 @@ function get_coherence(ψ::QuantumObject{Ket})
 end
 
 function get_coherence(ρ::QuantumObject{Operator})
-    a = destroy(prod(ρ.dimensions))
+    a = destroy(hilbert_dimensions_to_size(ρ.dimensions)[1])
     α = expect(a, ρ)
     D = exp(α * a' - conj(α) * a)
 
@@ -740,14 +740,14 @@ end
 _dims_and_perm(::ObjType, dims::SVector{N,Int}, order::AbstractVector{Int}, L::Int) where {ObjType<:Union{Ket,Bra},N} =
     reverse(dims), reverse((L + 1) .- order)
 
-# if dims originates from Dimensions
+# if dims originates from ProductDimensions
 _dims_and_perm(::Operator, 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
+# if dims originates from GeneralProductDimensions
 _dims_and_perm(::Operator, 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])
+_order_dimensions(dimensions::ProductDimensions, order::AbstractVector{Int}) = ProductDimensions(dimensions.to[order])
+_order_dimensions(dimensions::GeneralProductDimensions, order::AbstractVector{Int}) =
+    GeneralProductDimensions(dimensions.to[order], dimensions.from[order])
diff --git a/src/qobj/dimensions.jl b/src/qobj/dimensions.jl
index 692da0f05..c1cf19093 100644
--- a/src/qobj/dimensions.jl
+++ b/src/qobj/dimensions.jl
@@ -1,56 +1,57 @@
 #=
-This file defines the Dimensions structures, which can describe composite Hilbert spaces.
+This file defines the ProductDimensions structures, which can describe composite Hilbert spaces.
 =#
 
-export AbstractDimensions, Dimensions, GeneralDimensions
+export AbstractDimensions, ProductDimensions, GeneralProductDimensions
+export hilbert_dimensions_to_size, liouville_dimensions_to_size
 
 abstract type AbstractDimensions{M,N} end
 
 @doc raw"""
-    struct Dimensions{N,T<:Tuple} <: AbstractDimensions{N, N}
+    struct ProductDimensions{N,T<:Tuple} <: AbstractDimensions{N, N}
         to::T
     end
 
-A structure that describes the Hilbert [`Space`](@ref) of each subsystems.
+A structure that embodies the [`AbstractSpace`](@ref) of each subsystem in a composite Hilbert space.
 """
-struct Dimensions{N,T<:Tuple} <: AbstractDimensions{N,N}
+struct ProductDimensions{N,T<:Tuple} <: AbstractDimensions{N,N}
     to::T
 
     # make sure the elements in the tuple are all AbstractSpace
-    Dimensions(to::NTuple{N,AbstractSpace}) where {N} = new{N,typeof(to)}(to)
+    ProductDimensions(to::NTuple{N,AbstractSpace}) where {N} = new{N,typeof(to)}(to)
 end
-function Dimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Integer,N}
+function ProductDimensions(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(Tuple(Space.(dims)))
+    return ProductDimensions(Tuple(HilbertSpace.(dims)))
 end
-Dimensions(dims::Int) = Dimensions(Space(dims))
-Dimensions(dims::DimType) where {DimType<:AbstractSpace} = Dimensions((dims,))
-Dimensions(dims::Any) = throw(
+ProductDimensions(dims::Int) = ProductDimensions(HilbertSpace(dims))
+ProductDimensions(dims::DimType) where {DimType<:AbstractSpace} = ProductDimensions((dims,))
+ProductDimensions(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}
+    struct GeneralProductDimensions{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).
+A structure that embodies the left-hand side (`to`) and right-hand side (`from`) [`AbstractSpace`](@ref) of a quantum object.
 """
-struct GeneralDimensions{M,N,T1<:Tuple,T2<:Tuple} <: AbstractDimensions{M,N}
+struct GeneralProductDimensions{M,N,T1<:Tuple,T2<:Tuple} <: AbstractDimensions{M,N}
     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{M,AbstractSpace}, from::NTuple{N,AbstractSpace}) where {M,N} =
+    GeneralProductDimensions(to::NTuple{M,AbstractSpace}, from::NTuple{N,AbstractSpace}) where {M,N} =
         new{M,N,typeof(to),typeof(from)}(to, from)
 end
-function GeneralDimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Union{AbstractVector,NTuple},N}
+function GeneralProductDimensions(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])
@@ -61,44 +62,86 @@ function GeneralDimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<
     (L1 > 0) || throw(DomainError(L1, "The length of `dims[1]` must be larger or equal to 1."))
     (L2 > 0) || throw(DomainError(L2, "The length of `dims[2]` must be larger or equal to 1."))
 
-    return GeneralDimensions(Tuple(Space.(dims[1])), Tuple(Space.(dims[2])))
+    return GeneralProductDimensions(Tuple(HilbertSpace.(dims[1])), Tuple(HilbertSpace.(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{<:Integer},NTuple{N,Integer}}) where {N} = ProductDimensions(dims)
 _gen_dimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Union{AbstractVector,NTuple},N} =
-    GeneralDimensions(dims)
-_gen_dimensions(dims::Any) = Dimensions(dims)
+    GeneralProductDimensions(dims)
+_gen_dimensions(dims::Any) = ProductDimensions(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) =
+dimensions_to_dims(dimensions::ProductDimensions) = dimensions_to_dims(dimensions.to)
+dimensions_to_dims(dimensions::GeneralProductDimensions) =
     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)
+"""
+    hilbert_dimensions_to_size(dimensions)
+
+Returns the matrix dimensions `(m, n)` of an [`Operator`](@ref) with the given `dimensions`.
+
+For [`ProductDimensions`](@ref), returns `(m, m)` where `m` is the product of all subsystem Hilbert space dimensions.
+For [`GeneralProductDimensions`](@ref), returns `(m, n)` where `m` is the product of the `to` dimensions
+and `n` is the product of the `from` dimensions.
+
+If `dimensions` is an `Integer` or a vector/tuple of `Integer`s, it is automatically treated as [`ProductDimensions`](@ref).
+"""
+function hilbert_dimensions_to_size(dimensions::ProductDimensions)
+    m = prod(hilbert_dimensions_to_size, dimensions.to)
+    return (m, m)
+end
+function hilbert_dimensions_to_size(dimensions::GeneralProductDimensions)
+    m = prod(hilbert_dimensions_to_size, dimensions.to)
+    n = prod(hilbert_dimensions_to_size, dimensions.from)
+    return (m, n)
+end
+hilbert_dimensions_to_size(dimensions::Union{<:Integer,AbstractVector{<:Integer},NTuple{N,Integer}}) where {N} =
+    hilbert_dimensions_to_size(ProductDimensions(dimensions))
+
+"""
+    liouville_dimensions_to_size(dimensions)
+
+Returns the matrix dimensions `(m, n)` of a [`SuperOperator`](@ref) with the given `dimensions`.
+
+For [`ProductDimensions`](@ref), returns `(m, m)` where `m` is the product of all subsystem Liouville space dimensions
+(each Hilbert dimension `d` contributes `d²` to the product).
+For [`GeneralProductDimensions`](@ref), returns `(m, n)` where `m` is the product of the `to` dimensions
+and `n` is the product of the `from` dimensions.
+
+If `dimensions` is an `Integer` or a vector/tuple of `Integer`s, it is automatically treated as [`ProductDimensions`](@ref).
+"""
+function liouville_dimensions_to_size(dimensions::ProductDimensions)
+    m = prod(liouville_dimensions_to_size, dimensions.to)
+    return (m, m)
+end
+function liouville_dimensions_to_size(dimensions::GeneralProductDimensions)
+    m = prod(liouville_dimensions_to_size, dimensions.to)
+    n = prod(liouville_dimensions_to_size, dimensions.from)
+    return (m, n)
+end
+liouville_dimensions_to_size(dimensions::Union{<:Integer,AbstractVector{<:Integer},NTuple{N,Integer}}) where {N} =
+    liouville_dimensions_to_size(ProductDimensions(dimensions))
 
-Base.transpose(dimensions::Dimensions) = dimensions
-Base.transpose(dimensions::GeneralDimensions) = GeneralDimensions(dimensions.from, dimensions.to) # switch `to` and `from`
+Base.transpose(dimensions::ProductDimensions) = dimensions
+Base.transpose(dimensions::GeneralProductDimensions) = GeneralProductDimensions(dimensions.from, dimensions.to) # switch `to` and `from`
 Base.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)
+_get_dims_string(dimensions::ProductDimensions) = string(dimensions_to_dims(dimensions))
+function _get_dims_string(dimensions::GeneralProductDimensions)
     dims = dimensions_to_dims(dimensions)
     return "[$(string(dims[1])), $(string(dims[2]))]"
 end
 _get_dims_string(::Nothing) = "nothing" # for EigsolveResult.dimensions = nothing
 
-Base.:(==)(dim1::Dimensions, dim2::Dimensions) = dim1.to == dim2.to
-Base.:(==)(dim1::GeneralDimensions, dim2::GeneralDimensions) = (dim1.to == dim2.to) && (dim1.from == dim2.from)
-Base.:(==)(dim1::Dimensions, dim2::GeneralDimensions) = false
-Base.:(==)(dim1::GeneralDimensions, dim2::Dimensions) = false
+Base.:(==)(dim1::ProductDimensions, dim2::ProductDimensions) = dim1.to == dim2.to
+Base.:(==)(dim1::GeneralProductDimensions, dim2::GeneralProductDimensions) =
+    (dim1.to == dim2.to) && (dim1.from == dim2.from)
+Base.:(==)(dim1::ProductDimensions, dim2::GeneralProductDimensions) = false
+Base.:(==)(dim1::GeneralProductDimensions, dim2::ProductDimensions) = false
diff --git a/src/qobj/eigsolve.jl b/src/qobj/eigsolve.jl
index dabd91a56..c86993f14 100644
--- a/src/qobj/eigsolve.jl
+++ b/src/qobj/eigsolve.jl
@@ -45,7 +45,7 @@ julia> λ
   1.0 + 0.0im
 
 julia> ψ
-2-element Vector{QuantumObject{Ket, Dimensions{1, Tuple{Space}}, Vector{ComplexF64}}}:
+2-element Vector{QuantumObject{Ket, ProductDimensions{1, Tuple{HilbertSpace}}, Vector{ComplexF64}}}:
 
 Quantum Object:   type=Ket()   dims=[2]   size=(2,)
 2-element Vector{ComplexF64}:
@@ -424,7 +424,7 @@ end
         c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         alg::AbstractODEAlgorithm = DP5(),
         params::NamedTuple = NamedTuple(),
-        ρ0::AbstractMatrix = rand_dm(prod(H.dimensions)).data,
+        ρ0::AbstractMatrix = rand_dm(hilbert_dimensions_to_size(H.dimensions)[1]).data,
         eigvals::Int = 1,
         krylovdim::Int = min(10, size(H, 1)),
         maxiter::Int = 200,
@@ -467,7 +467,7 @@ function eigsolve_al(
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::AbstractODEAlgorithm = DP5(),
     params::NamedTuple = NamedTuple(),
-    ρ0::AbstractMatrix = rand_dm(prod(H.dimensions)).data,
+    ρ0::AbstractMatrix = rand_dm(hilbert_dimensions_to_size(H.dimensions)[1]).data,
     eigvals::Int = 1,
     krylovdim::Int = min(10, size(H, 1)),
     maxiter::Int = 200,
diff --git a/src/qobj/energy_restricted.jl b/src/qobj/energy_restricted.jl
index f65ea84ca..0e343b746 100644
--- a/src/qobj/energy_restricted.jl
+++ b/src/qobj/energy_restricted.jl
@@ -70,10 +70,15 @@ function Base.show(io::IO, s::EnrSpace)
     return nothing
 end
 
+Base.length(::EnrSpace{N}) where {N} = N
+
 Base.:(==)(s_enr1::EnrSpace, s_enr2::EnrSpace) = (s_enr1.size == s_enr2.size) && (s_enr1.dims == s_enr2.dims)
 
 dimensions_to_dims(s_enr::EnrSpace) = s_enr.dims
 
+hilbert_dimensions_to_size(s_enr::EnrSpace) = s_enr.size
+liouville_dimensions_to_size(s_enr::EnrSpace) = s_enr.size^2
+
 @doc raw"""
     enr_state_dictionaries(dims, n_excitations)
 
diff --git a/src/qobj/functions.jl b/src/qobj/functions.jl
index 1d2cf5689..644ef1d59 100644
--- a/src/qobj/functions.jl
+++ b/src/qobj/functions.jl
@@ -187,18 +187,18 @@ julia> a.dims, O.dims
 ```
 """
 function Base.kron(
-    A::AbstractQuantumObject{OpType,<:Dimensions},
-    B::AbstractQuantumObject{OpType,<:Dimensions},
+    A::AbstractQuantumObject{OpType,<:ProductDimensions},
+    B::AbstractQuantumObject{OpType,<:ProductDimensions},
 ) where {OpType<:Union{Ket,Bra,Operator}}
     QType = promote_op_type(A, B)
     _lazy_tensor_warning(A.data, B.data)
-    return QType(kron(A.data, B.data), A.type, Dimensions((A.dimensions.to..., B.dimensions.to...)))
+    return QType(kron(A.data, B.data), A.type, ProductDimensions((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
+# if A and B are both Operator but either one of them has GeneralProductDimensions
+for ADimType in (:ProductDimensions, :GeneralProductDimensions)
+    for BDimType in (:ProductDimensions, :GeneralProductDimensions)
+        if !(ADimType == BDimType == :ProductDimensions) # not for this case because it's already implemented
             @eval begin
                 function Base.kron(
                     A::AbstractQuantumObject{Operator,<:$ADimType},
@@ -209,7 +209,7 @@ for ADimType in (:Dimensions, :GeneralDimensions)
                     return QType(
                         kron(A.data, B.data),
                         Operator(),
-                        GeneralDimensions(
+                        GeneralProductDimensions(
                             (get_dimensions_to(A)..., get_dimensions_to(B)...),
                             (get_dimensions_from(A)..., get_dimensions_from(B)...),
                         ),
@@ -220,7 +220,7 @@ for ADimType in (:Dimensions, :GeneralDimensions)
     end
 end
 
-# if A and B are different type (must return Operator with GeneralDimensions)
+# if A and B are different type (must return Operator with GeneralProductDimensions)
 for AOpType in (:Ket, :Bra, :Operator)
     for BOpType in (:Ket, :Bra, :Operator)
         if (AOpType != BOpType)
@@ -231,7 +231,7 @@ for AOpType in (:Ket, :Bra, :Operator)
                     return QType(
                         kron(A.data, B.data),
                         Operator(),
-                        GeneralDimensions(
+                        GeneralProductDimensions(
                             (get_dimensions_to(A)..., get_dimensions_to(B)...),
                             (get_dimensions_from(A)..., get_dimensions_from(B)...),
                         ),
diff --git a/src/qobj/operators.jl b/src/qobj/operators.jl
index fb0c2041f..389990a46 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{Dimensions,AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+- `dimensions::Union{ProductDimensions,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,10 +33,12 @@ 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{Dimensions,AbstractVector{Int},Tuple}, distribution::Union{Symbol,Val} = Val(:haar)) =
-    rand_unitary(dimensions, makeVal(distribution))
-function rand_unitary(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}, ::Val{:haar})
-    N = prod(dimensions)
+rand_unitary(
+    dimensions::Union{ProductDimensions,AbstractVector{Int},Tuple},
+    distribution::Union{Symbol,Val} = Val(:haar),
+) = rand_unitary(dimensions, makeVal(distribution))
+function rand_unitary(dimensions::Union{ProductDimensions,AbstractVector{Int},Tuple}, ::Val{:haar})
+    N = hilbert_dimensions_to_size(dimensions)[1]
 
     # generate N x N matrix Z of complex standard normal random variates
     Z = randn(ComplexF64, N, N)
@@ -50,8 +52,8 @@ function rand_unitary(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}, :
     Λ ./= abs.(Λ) # rescaling the elements
     return QuantumObject(to_dense(Q * Diagonal(Λ)); type = Operator(), dims = dimensions)
 end
-function rand_unitary(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}, ::Val{:exp})
-    N = prod(dimensions)
+function rand_unitary(dimensions::Union{ProductDimensions,AbstractVector{Int},Tuple}, ::Val{:exp})
+    N = hilbert_dimensions_to_size(dimensions)[1]
 
     # generate N x N matrix Z of complex standard normal random variates
     Z = randn(ComplexF64, N, N)
@@ -61,7 +63,7 @@ function rand_unitary(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}, :
 
     return to_dense(exp(-1.0im * H))
 end
-rand_unitary(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}, ::Val{T}) where {T} =
+rand_unitary(dimensions::Union{ProductDimensions,AbstractVector{Int},Tuple}, ::Val{T}) where {T} =
     throw(ArgumentError("Invalid distribution: $(T)"))
 
 @doc raw"""
@@ -530,7 +532,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{Dimensions,AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+- `dimensions::Union{ProductDimensions,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
 
@@ -551,8 +553,8 @@ where ``\omega = \exp(\frac{2 \pi i}{N})``.
     It is highly recommended to use `qft(dimensions)` with `dimensions` as `Tuple` or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) 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{Dimensions,AbstractVector{Int},Tuple}) =
-    QuantumObject(_qft_op(prod(dimensions)), Operator(), dimensions)
+qft(dimensions::Union{ProductDimensions,AbstractVector{Int},Tuple}) =
+    QuantumObject(_qft_op(hilbert_dimensions_to_size(dimensions)[1]), Operator(), dimensions)
 function _qft_op(N::Int)
     ω = exp(2.0im * π / N)
     arr = 0:(N-1)
diff --git a/src/qobj/quantum_object.jl b/src/qobj/quantum_object.jl
index 053b9b050..0fb46d225 100644
--- a/src/qobj/quantum_object.jl
+++ b/src/qobj/quantum_object.jl
@@ -41,7 +41,7 @@ julia> a.dims
  20
 
 julia> a.dimensions
-Dimensions{1, Tuple{Space}}((Space(20),))
+ProductDimensions{1, Tuple{HilbertSpace}}((HilbertSpace(20),))
 ```
 """
 struct QuantumObject{ObjType<:QuantumObjectType,DimType<:AbstractDimensions,DataType<:AbstractArray} <:
@@ -55,8 +55,7 @@ struct QuantumObject{ObjType<:QuantumObjectType,DimType<:AbstractDimensions,Data
 
         ObjType = _check_type(type)
 
-        _size = _get_size(data)
-        _check_QuantumObject(type, dimensions, _size[1], _size[2])
+        _check_QuantumObject(type, dimensions, size(data, 1), size(data, 2))
 
         return new{ObjType,typeof(dimensions),DT}(data, type, dimensions)
     end
@@ -72,12 +71,10 @@ Generate [`QuantumObject`](@ref) with a given `A::AbstractArray` and specified `
     `Qobj` is a synonym of `QuantumObject`.
 """
 function QuantumObject(A::AbstractMatrix{T}; type = nothing, dims = nothing) where {T}
-    _size = _get_size(A)
-
     _check_type(type)
 
     if type isa Nothing
-        type = (_size[1] == 1 && _size[2] > 1) ? Bra() : Operator() # default type
+        type = (size(A, 1) == 1 && size(A, 2) > 1) ? Bra() : Operator() # default type
     elseif !(type isa Operator) && !(type isa SuperOperator) && !(type isa Bra) && !(type isa OperatorBra)
         throw(
             ArgumentError(
@@ -88,13 +85,13 @@ function QuantumObject(A::AbstractMatrix{T}; type = nothing, dims = nothing) whe
 
     if dims isa Nothing
         if type isa Bra
-            dims = Dimensions(_size[2])
+            dims = ProductDimensions(size(A, 2))
         elseif type isa Operator
             dims =
-                (_size[1] == _size[2]) ? Dimensions(_size[1]) :
-                GeneralDimensions(SVector{2}(SVector{1}(_size[1]), SVector{1}(_size[2])))
+                (size(A, 1) == size(A, 2)) ? ProductDimensions(size(A, 1)) :
+                GeneralProductDimensions(SVector{2}(SVector{1}(size(A, 1)), SVector{1}(size(A, 2))))
         elseif type isa SuperOperator || type isa OperatorBra
-            dims = Dimensions(isqrt(_size[2]))
+            dims = ProductDimensions(isqrt(size(A, 2)))
         end
     end
 
@@ -110,11 +107,10 @@ function QuantumObject(A::AbstractVector{T}; type = nothing, dims = nothing) whe
     end
 
     if dims isa Nothing
-        _size = _get_size(A)
         if type isa Ket
-            dims = Dimensions(_size[1])
+            dims = ProductDimensions(size(A, 1))
         elseif type isa OperatorKet
-            dims = Dimensions(isqrt(_size[1]))
+            dims = ProductDimensions(isqrt(size(A, 1)))
         end
     end
 
@@ -126,10 +122,9 @@ function QuantumObject(A::AbstractArray{T,N}; type = nothing, dims = nothing) wh
 end
 
 function QuantumObject(A::QuantumObject; type = A.type, dims = A.dimensions)
-    _size = _get_size(A.data)
     dimensions = _gen_dimensions(dims)
     _check_type(type)
-    _check_QuantumObject(type, dimensions, _size[1], _size[2])
+    _check_QuantumObject(type, dimensions, size(A.data, 1), size(A.data, 2))
     return QuantumObject(copy(A.data), type, dimensions)
 end
 
diff --git a/src/qobj/quantum_object_base.jl b/src/qobj/quantum_object_base.jl
index 6ae1024d3..3470a7790 100644
--- a/src/qobj/quantum_object_base.jl
+++ b/src/qobj/quantum_object_base.jl
@@ -130,7 +130,7 @@ check_dimensions(A::AbstractQuantumObject...) = check_dimensions(A)
 
 _check_QuantumObject(
     type::ObjType,
-    dimensions::GeneralDimensions,
+    dimensions::GeneralProductDimensions,
     m::Int,
     n::Int,
 ) where {ObjType<:Union{Ket,Bra,SuperOperator,OperatorBra,OperatorKet}} = throw(
@@ -140,25 +140,25 @@ _check_QuantumObject(
     ),
 )
 
-function _check_QuantumObject(type::Ket, dimensions::Dimensions, m::Int, n::Int)
+function _check_QuantumObject(type::Ket, dimensions::ProductDimensions, m::Int, n::Int)
     (n != 1) && throw(DomainError((m, n), "The size of the array is not compatible with Ket"))
-    (prod(dimensions) != m) && throw(
+    (hilbert_dimensions_to_size(dimensions)[1] != 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::Bra, dimensions::Dimensions, m::Int, n::Int)
+function _check_QuantumObject(type::Bra, dimensions::ProductDimensions, m::Int, n::Int)
     (m != 1) && throw(DomainError((m, n), "The size of the array is not compatible with Bra"))
-    (prod(dimensions) != n) && throw(
+    (hilbert_dimensions_to_size(dimensions)[1] != 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::Operator, dimensions::Dimensions, m::Int, n::Int)
-    L = prod(dimensions)
-    (L == m == n) || throw(
+function _check_QuantumObject(type::Operator, dimensions::ProductDimensions, m::Int, n::Int)
+    obj_size = hilbert_dimensions_to_size(dimensions)[1]
+    (obj_size == m == n) || throw(
         DimensionMismatch(
             "Operator with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
         ),
@@ -166,9 +166,10 @@ function _check_QuantumObject(type::Operator, dimensions::Dimensions, m::Int, n:
     return nothing
 end
 
-function _check_QuantumObject(type::Operator, dimensions::GeneralDimensions, m::Int, n::Int)
+function _check_QuantumObject(type::Operator, dimensions::GeneralProductDimensions, 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(
+    obj_size = hilbert_dimensions_to_size(dimensions)
+    ((obj_size[1] != m) || (obj_size[2] != n)) && throw(
         DimensionMismatch(
             "Operator with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
         ),
@@ -176,9 +177,9 @@ function _check_QuantumObject(type::Operator, dimensions::GeneralDimensions, m::
     return nothing
 end
 
-function _check_QuantumObject(type::SuperOperator, dimensions::Dimensions, m::Int, n::Int)
+function _check_QuantumObject(type::SuperOperator, dimensions::ProductDimensions, m::Int, n::Int)
     (m != n) && throw(DomainError((m, n), "The size of the array is not compatible with SuperOperator"))
-    (prod(dimensions) != sqrt(m)) && throw(
+    (liouville_dimensions_to_size(dimensions)[1] != m) && throw(
         DimensionMismatch(
             "SuperOperator with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
         ),
@@ -186,9 +187,9 @@ function _check_QuantumObject(type::SuperOperator, dimensions::Dimensions, m::In
     return nothing
 end
 
-function _check_QuantumObject(type::OperatorKet, dimensions::Dimensions, m::Int, n::Int)
+function _check_QuantumObject(type::OperatorKet, dimensions::ProductDimensions, m::Int, n::Int)
     (n != 1) && throw(DomainError((m, n), "The size of the array is not compatible with OperatorKet"))
-    (prod(dimensions) != sqrt(m)) && throw(
+    (liouville_dimensions_to_size(dimensions)[1] != m) && throw(
         DimensionMismatch(
             "OperatorKet with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
         ),
@@ -196,9 +197,9 @@ function _check_QuantumObject(type::OperatorKet, dimensions::Dimensions, m::Int,
     return nothing
 end
 
-function _check_QuantumObject(type::OperatorBra, dimensions::Dimensions, m::Int, n::Int)
+function _check_QuantumObject(type::OperatorBra, dimensions::ProductDimensions, m::Int, n::Int)
     (m != 1) && throw(DomainError((m, n), "The size of the array is not compatible with OperatorBra"))
-    (prod(dimensions) != sqrt(n)) && throw(
+    (liouville_dimensions_to_size(dimensions)[1] != n) && throw(
         DimensionMismatch(
             "OperatorBra with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
         ),
@@ -220,30 +221,24 @@ function Base.getproperty(A::AbstractQuantumObject, key::Symbol)
     end
 end
 
-# this returns `to` in GeneralDimensions representation
-get_dimensions_to(A::AbstractQuantumObject{Ket,<:Dimensions}) = A.dimensions.to
-get_dimensions_to(A::AbstractQuantumObject{Bra,<:Dimensions}) = space_one_list(A.dimensions.to)
-get_dimensions_to(A::AbstractQuantumObject{Operator,<:Dimensions}) = A.dimensions.to
-get_dimensions_to(A::AbstractQuantumObject{Operator,<:GeneralDimensions}) = A.dimensions.to
+# this returns `to` in GeneralProductDimensions representation
+get_dimensions_to(A::AbstractQuantumObject{Ket,<:ProductDimensions}) = A.dimensions.to
+get_dimensions_to(A::AbstractQuantumObject{Bra,<:ProductDimensions}) = hilbertspace_one_list(A.dimensions.to)
+get_dimensions_to(A::AbstractQuantumObject{Operator,<:ProductDimensions}) = A.dimensions.to
+get_dimensions_to(A::AbstractQuantumObject{Operator,<:GeneralProductDimensions}) = A.dimensions.to
 get_dimensions_to(
-    A::AbstractQuantumObject{ObjType,<:Dimensions},
+    A::AbstractQuantumObject{ObjType,<:ProductDimensions},
 ) where {ObjType<:Union{SuperOperator,OperatorBra,OperatorKet}} = A.dimensions.to
 
-# this returns `from` in GeneralDimensions representation
-get_dimensions_from(A::AbstractQuantumObject{Ket,<:Dimensions}) = space_one_list(A.dimensions.to)
-get_dimensions_from(A::AbstractQuantumObject{Bra,<:Dimensions}) = A.dimensions.to
-get_dimensions_from(A::AbstractQuantumObject{Operator,<:Dimensions}) = A.dimensions.to
-get_dimensions_from(A::AbstractQuantumObject{Operator,<:GeneralDimensions}) = A.dimensions.from
+# this returns `from` in GeneralProductDimensions representation
+get_dimensions_from(A::AbstractQuantumObject{Ket,<:ProductDimensions}) = hilbertspace_one_list(A.dimensions.to)
+get_dimensions_from(A::AbstractQuantumObject{Bra,<:ProductDimensions}) = A.dimensions.to
+get_dimensions_from(A::AbstractQuantumObject{Operator,<:ProductDimensions}) = A.dimensions.to
+get_dimensions_from(A::AbstractQuantumObject{Operator,<:GeneralProductDimensions}) = A.dimensions.from
 get_dimensions_from(
-    A::AbstractQuantumObject{ObjType,<:Dimensions},
+    A::AbstractQuantumObject{ObjType,<:ProductDimensions},
 ) where {ObjType<:Union{SuperOperator,OperatorBra,OperatorKet}} = A.dimensions.to
 
-# this creates a list of Space(1), it is used to generate `from` for Ket, and `to` for Bra
-space_one_list(dimensions::NTuple{N,AbstractSpace}) where {N} =
-    ntuple(i -> Space(1), Val(sum(_get_dims_length, dimensions)))
-_get_dims_length(::Space) = 1
-_get_dims_length(::EnrSpace{N}) where {N} = N
-
 # functions for getting Float or Complex element type
 _float_type(A::AbstractQuantumObject) = _float_type(eltype(A))
 _complex_float_type(A::AbstractQuantumObject) = _complex_float_type(eltype(A))
diff --git a/src/qobj/quantum_object_evo.jl b/src/qobj/quantum_object_evo.jl
index 8e4de483b..a2d1016e8 100644
--- a/src/qobj/quantum_object_evo.jl
+++ b/src/qobj/quantum_object_evo.jl
@@ -124,8 +124,7 @@ struct QuantumObjectEvolution{
 
         dimensions = _gen_dimensions(dims)
 
-        _size = _get_size(data)
-        _check_QuantumObject(type, dimensions, _size[1], _size[2])
+        _check_QuantumObject(type, dimensions, size(data, 1), size(data, 2))
 
         return new{ObjType,typeof(dimensions),DT}(data, type, dimensions)
     end
@@ -158,16 +157,15 @@ Generate a [`QuantumObjectEvolution`](@ref) object from a [`SciMLOperator`](http
 Note that `QobjEvo` is a synonym of `QuantumObjectEvolution`
 """
 function QuantumObjectEvolution(data::AbstractSciMLOperator; type = Operator(), dims = nothing)
-    _size = _get_size(data)
     _check_type(type)
 
     if dims isa Nothing
         if type isa Operator
             dims =
-                (_size[1] == _size[2]) ? Dimensions(_size[1]) :
-                GeneralDimensions(SVector{2}(SVector{1}(_size[1]), SVector{1}(_size[2])))
+                (size(data, 1) == size(data, 2)) ? ProductDimensions(size(data, 1)) :
+                GeneralProductDimensions(SVector{2}(SVector{1}(size(data, 1)), SVector{1}(size(data, 2))))
         elseif type isa SuperOperator
-            dims = Dimensions(isqrt(_size[2]))
+            dims = ProductDimensions(isqrt(size(data, 2)))
         end
     end
 
diff --git a/src/qobj/space.jl b/src/qobj/space.jl
index 90b199759..9cc051bbd 100644
--- a/src/qobj/space.jl
+++ b/src/qobj/space.jl
@@ -2,24 +2,33 @@
 This file defines the Hilbert space structure.
 =#
 
-export AbstractSpace, Space
+export AbstractSpace, HilbertSpace
 
 abstract type AbstractSpace end
 
 @doc raw"""
-    struct Space <: AbstractSpace
+    struct HilbertSpace <: AbstractSpace
         size::Int
     end
 
 A structure that describes a single Hilbert space with size = `size`.
 """
-struct Space <: AbstractSpace
+struct HilbertSpace <: AbstractSpace
     size::Int
 
-    function Space(size::Int)
-        (size < 1) && throw(DomainError(size, "The size of Space must be positive integer (≥ 1)."))
+    function HilbertSpace(size::Int)
+        (size < 1) && throw(DomainError(size, "The size of `HilbertSpace` must be positive integer (≥ 1)."))
         return new(size)
     end
 end
 
-dimensions_to_dims(s::Space) = SVector{1,Int}(s.size)
+Base.length(s::HilbertSpace) = 1
+
+dimensions_to_dims(s::HilbertSpace) = SVector{1,Int}(s.size)
+
+hilbert_dimensions_to_size(s::HilbertSpace) = s.size
+liouville_dimensions_to_size(s::HilbertSpace) = s.size^2
+
+# this creates a list of HilbertSpace(1), it is used to generate `from` for Ket, and `to` for Bra
+hilbertspace_one_list(dimensions::NTuple{N,AbstractSpace}) where {N} =
+    ntuple(i -> HilbertSpace(1), Val(sum(length, dimensions)))
diff --git a/src/qobj/states.jl b/src/qobj/states.jl
index dd79bee8d..1981ced98 100644
--- a/src/qobj/states.jl
+++ b/src/qobj/states.jl
@@ -14,14 +14,14 @@ 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{Dimensions,AbstractVector{Int}, Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+- `dimensions::Union{ProductDimensions,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` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) 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{Dimensions,AbstractVector{Int},Tuple}) =
-    QuantumObject(zeros(ComplexF64, prod(dimensions)), Ket(), dimensions)
+zero_ket(dimensions::Union{ProductDimensions,AbstractVector{Int},Tuple}) =
+    QuantumObject(zeros(ComplexF64, hilbert_dimensions_to_size(dimensions)[1]), Ket(), dimensions)
 
 @doc raw"""
     fock(N::Int, j::Int=0; dims::Union{Int,AbstractVector{Int},Tuple}=N, sparse::Union{Bool,Val}=Val(false))
@@ -70,14 +70,14 @@ 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{Dimensions,AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+- `dimensions::Union{ProductDimensions,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` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) 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{Dimensions,AbstractVector{Int},Tuple})
-    N = prod(dimensions)
+function rand_ket(dimensions::Union{ProductDimensions,AbstractVector{Int},Tuple})
+    N = hilbert_dimensions_to_size(dimensions)[1]
     ψ = rand(ComplexF64, N) .- (0.5 + 0.5im)
     return QuantumObject(normalize!(ψ); type = Ket(), dims = dimensions)
 end
@@ -142,28 +142,28 @@ 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{Dimensions,AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+- `dimensions::Union{ProductDimensions,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` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 maximally_mixed_dm(dimensions::Int) =
     QuantumObject(diagm(0 => fill(ComplexF64(1 / dimensions), dimensions)), Operator(), SVector(dimensions))
-function maximally_mixed_dm(dimensions::Union{Dimensions,AbstractVector{Int},Tuple})
-    N = prod(dimensions)
+function maximally_mixed_dm(dimensions::Union{ProductDimensions,AbstractVector{Int},Tuple})
+    N = hilbert_dimensions_to_size(dimensions)[1]
     return QuantumObject(diagm(0 => fill(ComplexF64(1 / N), N)), Operator(), dimensions)
 end
 
 @doc raw"""
-    rand_dm(dimensions; rank::Int=prod(dimensions))
+    rand_dm(dimensions; rank::Int=hilbert_dimensions_to_size(dimensions)[1])
 
 Generate a random density matrix from Ginibre ensemble with given argument `dimensions` and `rank`, ensuring that it is positive semi-definite and trace equals to `1`.
 
 The `dimensions` can be either the following types:
 - `dimensions::Int`: Number of basis states in the Hilbert space.
-- `dimensions::Union{Dimensions,AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+- `dimensions::Union{ProductDimensions,AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
 
-The default keyword argument `rank = prod(dimensions)` (full rank).
+The default keyword argument `rank = hilbert_dimensions_to_size(dimensions)[1]` (full rank).
 
 !!! warning "Beware of type-stability!"
     If you want to keep type stability, it is recommended to use `rand_dm(dimensions; rank=rank)` with `dimensions` as `Tuple` or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) instead of `Vector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
@@ -172,9 +172,13 @@ The default keyword argument `rank = prod(dimensions)` (full rank).
 - [J. Ginibre, Statistical ensembles of complex, quaternion, and real matrices, Journal of Mathematical Physics 6.3 (1965): 440-449](https://doi.org/10.1063/1.1704292)
 - [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{Dimensions,AbstractVector{Int},Tuple}; rank::Int = prod(dimensions))
-    N = prod(dimensions)
+rand_dm(dimensions::Int; rank::Int = hilbert_dimensions_to_size(dimensions)[1]) =
+    rand_dm(SVector(dimensions), rank = rank)
+function rand_dm(
+    dimensions::Union{ProductDimensions,AbstractVector{Int},Tuple};
+    rank::Int = hilbert_dimensions_to_size(dimensions)[1],
+)
+    N = hilbert_dimensions_to_size(dimensions)[1]
     (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/spectrum.jl b/src/spectrum.jl
index a4ce15d67..772693b9f 100644
--- a/src/spectrum.jl
+++ b/src/spectrum.jl
@@ -136,7 +136,7 @@ function _spectrum(
     b = (spre(B) * ρss).data
 
     # multiply by operator A on the left (spre) and then perform trace operation
-    D = prod(L.dimensions)
+    D = hilbert_dimensions_to_size(L.dimensions)[1]
     _tr = SparseVector(D^2, [1 + n * (D + 1) for n in 0:(D-1)], ones(_complex_float_type(L), D)) # same as vec(system_identity_matrix)
     _tr_A = transpose(_tr) * spre(A).data
 
@@ -182,7 +182,7 @@ function _spectrum(
     vT = typeof(vₖ)
 
     # Calculate <w₁| = <I|A
-    D = prod(L.dimensions)
+    D = hilbert_dimensions_to_size(L.dimensions)[1]
     Ivec = SparseVector(D^2, [1 + n * (D + 1) for n in 0:(D-1)], ones(cT, D)) # same as vec(system_identity_matrix)
     wₖ = spre(A).data' * vT(Ivec)
 
diff --git a/src/steadystate.jl b/src/steadystate.jl
index d45916da7..2c988cb99 100644
--- a/src/steadystate.jl
+++ b/src/steadystate.jl
@@ -138,7 +138,7 @@ end
 
 function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateLinearSolver; kwargs...)
     L_tmp = L.data
-    N = prod(L.dimensions)
+    N = hilbert_dimensions_to_size(L.dimensions)[1]
     weight = norm(L_tmp, 1) / length(L_tmp)
 
     v0 = _dense_similar(L_tmp, N^2)
@@ -166,7 +166,7 @@ function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateLinear
 end
 
 function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateEigenSolver; kwargs...)
-    N = prod(L.dimensions)
+    N = hilbert_dimensions_to_size(L.dimensions)[1]
 
     kwargs = merge((sigma = 1e-8, eigvals = 1), (; kwargs...))
 
@@ -179,7 +179,7 @@ end
 
 function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateDirectSolver)
     L_tmp = L.data
-    N = prod(L.dimensions)
+    N = hilbert_dimensions_to_size(L.dimensions)[1]
     weight = norm(L_tmp, 1) / length(L_tmp)
 
     v0 = _dense_similar(L_tmp, N^2)
diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl
index 1560de45f..abe9a4f7e 100644
--- a/src/time_evolution/smesolve.jl
+++ b/src/time_evolution/smesolve.jl
@@ -111,7 +111,7 @@ function smesolveProblem(
 
     K = cache_operator(get_data(L_evo), ρ0)
 
-    Id_op = IdentityOperator(prod(dims)^2)
+    Id_op = IdentityOperator(liouville_dimensions_to_size(dims)[1])
     D_l = map(sc_ops_evo_data) do op
         # TODO: # Currently, we are assuming a time-independent MatrixOperator
         # Also, the u state may become non-hermitian, so Tr[Sn * ρ + ρ * Sn'] != real(Tr[Sn * ρ]) / 2
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index df3372858..ab2a36eee 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -107,13 +107,18 @@ function ssesolveProblem(
 
     # Here the coefficients depend on the state, so this is a non-linear operator, which should be implemented with FunctionOperator instead. However, the nonlinearity is only on the coefficients, and it should be safe.
     K_l = sum(
-        op -> _ScalarOperator_e(op, +) * op + _ScalarOperator_e2_2(op, -) * IdentityOperator(prod(dims)),
+        op ->
+            _ScalarOperator_e(op, +) * op +
+            _ScalarOperator_e2_2(op, -) * IdentityOperator(hilbert_dimensions_to_size(dims)[1]),
         sc_ops_evo_data,
     )
 
     K = cache_operator(get_data(-1im * QuantumObjectEvolution(H_eff_evo)) + K_l, ψ0)
 
-    D_l = map(op -> op + _ScalarOperator_e(op, -) * IdentityOperator(prod(dims)), sc_ops_evo_data)
+    D_l = map(
+        op -> op + _ScalarOperator_e(op, -) * IdentityOperator(hilbert_dimensions_to_size(dims)[1]),
+        sc_ops_evo_data,
+    )
     D = DiffusionOperator(D_l)
 
     kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_SDE_SOLVER_OPTIONS; kwargs...)
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index e7c3791d3..fdd3f548b 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -564,7 +564,7 @@ function liouvillian_generalized(
     (length(fields) == length(T_list)) || throw(DimensionMismatch("The number of fields, ωs and Ts must be the same."))
 
     dims = (N_trunc isa Nothing) ? H.dimensions : SVector(N_trunc)
-    final_size = prod(dims)
+    final_size = hilbert_dimensions_to_size(dims)[1]
     result = eigen(H)
     E = real.(result.values[1:final_size])
     U = QuantumObject(result.vectors, result.type, result.dimensions)
diff --git a/src/time_evolution/time_evolution_dynamical.jl b/src/time_evolution/time_evolution_dynamical.jl
index bbe70a9f2..47fbd544b 100644
--- a/src/time_evolution/time_evolution_dynamical.jl
+++ b/src/time_evolution/time_evolution_dynamical.jl
@@ -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 = Eye(prod(H₀.dimensions))
+    dsf_identity = Eye(hilbert_dimensions_to_size(H₀.dimensions)[1])
     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/src/utilities.jl b/src/utilities.jl
index 6d647c08b..ba7d79a98 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -143,10 +143,6 @@ makeVal(x) = Val(x)
 getVal(x::Val{T}) where {T} = T
 getVal(x) = x # getVal for any other type
 
-_get_size(A::AbstractMatrix) = size(A)
-_get_size(A::AbstractVector) = (length(A), 1)
-_get_size(A::AbstractSciMLOperator) = size(A)
-
 _non_static_array_warning(argname, arg::Tuple{}) =
     throw(ArgumentError("The argument $argname must be a Tuple or a StaticVector of non-zero length."))
 _non_static_array_warning(argname, arg::Union{SVector{N,T},MVector{N,T},NTuple{N,T}}) where {N,T} = nothing
diff --git a/test/core-test/enr_state_operator.jl b/test/core-test/enr_state_operator.jl
index 50bed88c8..fd74af62c 100644
--- a/test/core-test/enr_state_operator.jl
+++ b/test/core-test/enr_state_operator.jl
@@ -22,14 +22,14 @@
     end
 
     @testset "kron" begin
-        # normal Space
+        # normal HilbertSpace
         D1 = 4
         D2 = 5
         dims_s = (D1, D2)
         ρ_s = rand_dm(dims_s)
         I_s = qeye(D1) ⊗ qeye(D2)
         size_s = prod(dims_s)
-        space_s = (Space(D1), Space(D2))
+        space_s = (HilbertSpace(D1), HilbertSpace(D2))
 
         # EnrSpace
         dims_enr = (2, 2, 3)
@@ -49,15 +49,20 @@
         @test opstring ==
               "\nQuantum Object:   type=Operator()   dims=$ρ_tot_dims   size=$((ρ_tot_size, ρ_tot_size))   ishermitian=$ρ_tot_isherm\n$datastring"
 
-        # use GeneralDimensions to do partial trace
-        new_dims1 = GeneralDimensions((Space(1), Space(1), space_enr), (Space(1), Space(1), space_enr))
+        # use GeneralProductDimensions to do partial trace
+        new_dims1 = GeneralProductDimensions(
+            (HilbertSpace(1), HilbertSpace(1), space_enr),
+            (HilbertSpace(1), HilbertSpace(1), space_enr),
+        )
         ρ_enr_compound = Qobj(zeros(ComplexF64, size_enr, size_enr), dims = new_dims1)
         basis_list = [tensor(basis(D1, i), basis(D2, j)) for i in 0:(D1-1) for j in 0:(D2-1)]
         for b in basis_list
             ρ_enr_compound += tensor(b', I_enr) * ρ_tot * tensor(b, I_enr)
         end
-        new_dims2 =
-            GeneralDimensions((space_s..., Space(1), Space(1), Space(1)), (space_s..., Space(1), Space(1), Space(1)))
+        new_dims2 = GeneralProductDimensions(
+            (space_s..., HilbertSpace(1), HilbertSpace(1), HilbertSpace(1)),
+            (space_s..., HilbertSpace(1), HilbertSpace(1), HilbertSpace(1)),
+        )
         ρ_s_compound = Qobj(zeros(ComplexF64, size_s, size_s), dims = new_dims2)
         basis_list = [enr_fock(space_enr, space_enr.idx2state[idx]) for idx in 1:space_enr.size]
         for b in basis_list
diff --git a/test/core-test/low_rank_dynamics.jl b/test/core-test/low_rank_dynamics.jl
index 099b8b590..bf1aa000a 100644
--- a/test/core-test/low_rank_dynamics.jl
+++ b/test/core-test/low_rank_dynamics.jl
@@ -11,7 +11,7 @@
     M = latt.N + 1       # Number of states in the LR basis
 
     # Define initial state
-    ϕ = Vector{QuantumObject{Ket,Dimensions{M - 1,NTuple{M - 1,Space}},Vector{ComplexF64}}}(undef, M)
+    ϕ = Vector{QuantumObject{Ket,ProductDimensions{M - 1,NTuple{M - 1,HilbertSpace}},Vector{ComplexF64}}}(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 a554d6065..2a3fc8a12 100644
--- a/test/core-test/negativity_and_partial_transpose.jl
+++ b/test/core-test/negativity_and_partial_transpose.jl
@@ -52,7 +52,7 @@
             end
         end
         @test_throws ArgumentError partial_transpose(A_dense, [true])
-        @test_throws ArgumentError partial_transpose(Qobj(zeros(ComplexF64, 3, 2)), [true]) # invalid GeneralDimensions
+        @test_throws ArgumentError partial_transpose(Qobj(zeros(ComplexF64, 3, 2)), [true]) # invalid GeneralProductDimensions
 
         @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 7551824d0..a78943d9c 100644
--- a/test/core-test/quantum_objects.jl
+++ b/test/core-test/quantum_objects.jl
@@ -32,7 +32,7 @@
 
         @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
+        # check that Ket, Bra, SuperOperator, OperatorKet, and OperatorBra don't support GeneralProductDimensions
         @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,)))
@@ -89,7 +89,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
+        a4 = Qobj(sprand(ComplexF64, 100, 10, 0.1)) # GeneralProductDimensions
         a5 = QuantumObject(rand(ComplexF64, 2*3*4, 5), dims = ((2, 3, 4), (5,)))
         @test isket(a2) == false
         @test isbra(a2) == false
@@ -271,7 +271,7 @@
         @test opstring ==
               "\nQuantum Object:   type=Operator()   dims=$a_dims   size=$a_size   ishermitian=$a_isherm\n$datastring"
 
-        # GeneralDimensions
+        # GeneralProductDimensions
         Gop = tensor(a, ψ)
         opstring = sprint((t, s) -> show(t, "text/plain", s), Gop)
         datastring = sprint((t, s) -> show(t, "text/plain", s), Gop.data)
@@ -367,8 +367,8 @@
             end
 
             UnionType = Union{
-                QuantumObject{Bra,Dimensions{1,Tuple{Space}},Matrix{T}},
-                QuantumObject{Operator,Dimensions{1,Tuple{Space}},Matrix{T}},
+                QuantumObject{Bra,ProductDimensions{1,Tuple{HilbertSpace}},Matrix{T}},
+                QuantumObject{Operator,ProductDimensions{1,Tuple{HilbertSpace}},Matrix{T}},
             }
             a = rand(T, 1, N)
             @inferred UnionType Qobj(a)
@@ -377,8 +377,8 @@
             end
 
             UnionType2 = Union{
-                QuantumObject{Operator,GeneralDimensions{1,1,Tuple{Space},Tuple{Space}},Matrix{T}},
-                QuantumObject{Operator,Dimensions{1,Tuple{Space}},Matrix{T}},
+                QuantumObject{Operator,GeneralProductDimensions{1,1,Tuple{HilbertSpace},Tuple{HilbertSpace}},Matrix{T}},
+                QuantumObject{Operator,ProductDimensions{1,Tuple{HilbertSpace}},Matrix{T}},
             }
             a = rand(T, N, N)
             @inferred UnionType Qobj(a)
@@ -670,7 +670,7 @@
         ρ1_ptr = ptrace(ρ, 1)
         ρ2_ptr = ptrace(ρ, 2)
 
-        # use GeneralDimensions to do partial trace
+        # use GeneralProductDimensions 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]
@@ -741,7 +741,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
+        @test_throws ArgumentError ptrace(Qobj(zeros(ComplexF64, 3, 2)), 1) # invalid GeneralProductDimensions
 
         @testset "Type Inference (ptrace)" begin
             @inferred ptrace(ρ, 1)
@@ -754,7 +754,7 @@
     end
 
     @testset "permute" begin
-        # standard Dimensions
+        # standard ProductDimensions
         ket_a = Qobj(rand(ComplexF64, 2))
         ket_b = Qobj(rand(ComplexF64, 3))
         ket_c = Qobj(rand(ComplexF64, 4))
@@ -789,7 +789,7 @@
         @test_throws ArgumentError permute(op_bdca, wrong_order1)
         @test_throws ArgumentError permute(op_bdca, wrong_order2)
 
-        # GeneralDimensions
+        # GeneralProductDimensions
         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))
diff --git a/test/core-test/quantum_objects_evo.jl b/test/core-test/quantum_objects_evo.jl
index 653d5a1ab..60aea7a61 100644
--- a/test/core-test/quantum_objects_evo.jl
+++ b/test/core-test/quantum_objects_evo.jl
@@ -139,8 +139,12 @@
         for T in [ComplexF32, ComplexF64]
             a = MatrixOperator(rand(T, N, N))
             UnionType = Union{
-                QuantumObjectEvolution{Operator,GeneralDimensions{1,Tuple{Space},Tuple{Space}},typeof(a)},
-                QuantumObjectEvolution{Operator,Dimensions{1,Tuple{Space}},typeof(a)},
+                QuantumObjectEvolution{
+                    Operator,
+                    GeneralProductDimensions{1,Tuple{HilbertSpace},Tuple{HilbertSpace}},
+                    typeof(a),
+                },
+                QuantumObjectEvolution{Operator,ProductDimensions{1,Tuple{HilbertSpace}},typeof(a)},
             }
             @inferred UnionType QobjEvo(a)
             @inferred UnionType QobjEvo(a, type = Operator())