diff --git a/CHANGELOG.md b/CHANGELOG.md
index 24cd4648f..63edc6772 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
+- Add `spre` and `spost` methods for `ComposedOperator` and cache propagator in every time evolution solver. ([#596])
+
 ## [v0.39.0]
 Release date: 2025-11-17
 
@@ -374,3 +376,4 @@ Release date: 2024-11-13
 [#588]: https://github.com/qutip/QuantumToolbox.jl/issues/588
 [#589]: https://github.com/qutip/QuantumToolbox.jl/issues/589
 [#591]: https://github.com/qutip/QuantumToolbox.jl/issues/591
+[#596]: https://github.com/qutip/QuantumToolbox.jl/issues/596
diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index 44b7c10e6..b54cf2e81 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -50,6 +50,7 @@ import SciMLOperators:
     ScalarOperator,
     ScaledOperator,
     AddedOperator,
+    ComposedOperator,
     IdentityOperator,
     update_coefficients!,
     concretize
diff --git a/src/qobj/arithmetic_and_attributes.jl b/src/qobj/arithmetic_and_attributes.jl
index 14e06aae4..6902e053f 100644
--- a/src/qobj/arithmetic_and_attributes.jl
+++ b/src/qobj/arithmetic_and_attributes.jl
@@ -38,6 +38,8 @@ for op in (:(+), :(-), :(*))
     @eval begin
         function Base.$op(A::AbstractQuantumObject, B::AbstractQuantumObject)
             check_dimensions(A, B)
+            A.type != B.type &&
+                throw(DomainError((A.type, B.type), "The quantum objects should have the same type to do $op."))
             QType = promote_op_type(A, B)
             return QType($(op)(A.data, B.data), A.type, A.dimensions)
         end
@@ -50,54 +52,24 @@ for op in (:(+), :(-), :(*))
     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.",
-        ),
+for type in (:Operator, :SuperOperator)
+    @eval function Base.:(*)(
+        A::AbstractQuantumObject{$type,ProductDimensions},
+        B::AbstractQuantumObject{$type,ProductDimensions},
     )
-    return nothing
-end
-
-for ADimType in (:Dimensions, :GeneralDimensions)
-    for BDimType in (:Dimensions, :GeneralDimensions)
-        if ADimType == BDimType == :Dimensions
-            @eval begin
-                function Base.:(*)(
-                    A::AbstractQuantumObject{Operator,<:$ADimType},
-                    B::AbstractQuantumObject{Operator,<:$BDimType},
-                )
-                    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 Base.:(*)(
-                    A::AbstractQuantumObject{Operator,<:$ADimType},
-                    B::AbstractQuantumObject{Operator,<:$BDimType},
-                )
-                    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
+        check_mul_dimensions(get_dimensions_from(A), get_dimensions_to(B))
+        QType = promote_op_type(A, B)
+        return QType(A.data * B.data, $type(), ProductDimensions(get_dimensions_to(A), get_dimensions_from(B)))
     end
 end
 
-function Base.:(*)(A::AbstractQuantumObject{Operator}, B::QuantumObject{Ket,<:Dimensions})
+function Base.:(*)(A::AbstractQuantumObject{Operator}, B::QuantumObject{Ket})
     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}, 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)
@@ -124,6 +96,15 @@ function Base.:(*)(A::QuantumObject{OperatorBra}, B::AbstractQuantumObject{Super
     return QuantumObject(A.data * B.data, OperatorBra(), A.dimensions)
 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
+
 Base.:(^)(A::QuantumObject, n::T) where {T<:Number} = QuantumObject(^(A.data, n), A.type, A.dimensions)
 Base.:(/)(A::AbstractQuantumObject, n::T) where {T<:Number} = get_typename_wrapper(A)(A.data / n, A.type, A.dimensions)
 
@@ -672,7 +653,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 +661,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)
 
@@ -748,6 +729,5 @@ _dims_and_perm(::Operator, dims::SVector{N,Int}, order::AbstractVector{Int}, L::
 _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], isnothing(dimensions.from) ? nothing : dimensions.from[order])
diff --git a/src/qobj/dimensions.jl b/src/qobj/dimensions.jl
index 692da0f05..0b7ec9f48 100644
--- a/src/qobj/dimensions.jl
+++ b/src/qobj/dimensions.jl
@@ -2,103 +2,100 @@
 This file defines the Dimensions structures, which can describe composite Hilbert spaces.
 =#
 
-export AbstractDimensions, Dimensions, GeneralDimensions
+export AbstractDimensions, ProductDimensions
 
 abstract type AbstractDimensions{M,N} end
 
 @doc raw"""
-    struct Dimensions{N,T<:Tuple} <: AbstractDimensions{N, N}
-        to::T
-    end
-
-A structure that describes the Hilbert [`Space`](@ref) of each subsystems.
-"""
-struct Dimensions{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)
-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(Tuple(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}
+    struct ProductDimensions{M,N,T1<:Tuple,T2<:Union{<:Tuple, Nothing}} <: AbstractDimensions{M,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 describes the left-hand side (`to`) and right-hand side (`from`) dimensions of a quantum object.
+
+The `from` field can be different from `nothing` only for non-square [`Operator`](@ref) and [`SuperOperator`](@ref) quantum objects.
 """
-struct GeneralDimensions{M,N,T1<:Tuple,T2<:Tuple} <: AbstractDimensions{M,N}
+struct ProductDimensions{M,N,T1<:Tuple,T2<:Union{<:Tuple,Nothing}} <: 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} =
-        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}
-    (length(dims) != 2) && throw(ArgumentError("Invalid dims = $dims"))
+    function ProductDimensions(to::Union{AbstractVector,NTuple}, from::Union{NTuple,Nothing})
+        M = length(to)
+        N = isnothing(from) ? M : length(from)
 
-    _non_static_array_warning("dims[1]", dims[1])
-    _non_static_array_warning("dims[2]", dims[2])
+        _non_static_array_warning("dims", to)
+        isnothing(from) || _non_static_array_warning("dims", from)
 
-    L1 = length(dims[1])
-    L2 = length(dims[2])
-    (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."))
+        to_space = _dims_tuple_of_space(to)
+        from_space = _dims_tuple_of_space(from)
 
-    return GeneralDimensions(Tuple(Space.(dims[1])), Tuple(Space.(dims[2])))
+        new{M,N,typeof(to_space),typeof(from_space)}(to_space, from_space)
+    end
 end
+function ProductDimensions(dims::Union{AbstractVector,Tuple})
+    (length(dims) != 2) && throw(ArgumentError("Invalid dims = $dims"))
 
-_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)
+    return ProductDimensions(dims[1], dims[2])
+end
+
+ProductDimensions(dims::Union{Int,AbstractSpace}) = ProductDimensions((dims,), nothing)
+ProductDimensions(dims::Union{AbstractVector{<:Integer},NTuple{N,Integer}}) where {N} = ProductDimensions(dims, nothing)
+ProductDimensions(dims::Union{AbstractVector{<:AbstractSpace},NTuple{N,AbstractSpace}}) where {N} = ProductDimensions(dims, nothing)
+ProductDimensions(dims::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) =
-    SVector{2}(dimensions_to_dims(dimensions.to), dimensions_to_dims(dimensions.from))
-
+function dimensions_to_dims(dimensions::ProductDimensions)
+    dims_to = dimensions_to_dims(dimensions.to)
+    isnothing(dimensions.from) && return dims_to
+    dims_from = dimensions_to_dims(dimensions.from)
+    return SVector{2}(dims_to, dims_from)
+end
 dimensions_to_dims(::Nothing) = nothing # for EigsolveResult.dimensions = nothing
 
-Base.length(::AbstractDimensions{N}) where {N} = N
+hilbert_dimensions_to_size(dimensions::ProductDimensions) =
+    (hilbert_dimensions_to_size(dimensions.to), hilbert_dimensions_to_size(dimensions.from))
+hilbert_dimensions_to_size(dimensions::NTuple{N,AbstractSpace}) where {N} = prod(hilbert_dimensions_to_size, dimensions)
+hilbert_dimensions_to_size(dim::Int) = dim
+hilbert_dimensions_to_size(dimensions::Union{AbstractVector, NTuple{N,Integer}}) where {N} = prod(dimensions)
+hilbert_dimensions_to_size(::Nothing) = nothing
+
+liouvillian_dimensions_to_size(dimensions::ProductDimensions) =
+    (liouvillian_dimensions_to_size(dimensions.to), liouvillian_dimensions_to_size(dimensions.from))
+liouvillian_dimensions_to_size(dimensions::NTuple{N,AbstractSpace}) where {N} =
+    prod(liouvillian_dimensions_to_size, dimensions)
+liouvillian_dimensions_to_size(::Nothing) = nothing
 
-# 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)
+Base.length(::AbstractDimensions{N}) where {N} = N
 
-Base.transpose(dimensions::Dimensions) = dimensions
-Base.transpose(dimensions::GeneralDimensions) = GeneralDimensions(dimensions.from, dimensions.to) # switch `to` and `from`
+Base.transpose(dimensions::ProductDimensions) = ProductDimensions(dimensions.from, dimensions.to) # switch `to` and `from`
+Base.transpose(dimensions::ProductDimensions{M,N,T1,Nothing}) where {M,N,T1<:Tuple} = dimensions
 Base.adjoint(dimensions::AbstractDimensions) = transpose(dimensions)
 
+Base.:(==)(dim1::ProductDimensions, dim2::ProductDimensions) = (dim1.to == dim2.to) && (dim1.from == dim2.from)
+
+_dims_tuple_of_space(dims::NTuple{N,AbstractSpace}) where {N} = dims
+function _dims_tuple_of_space(dims::Union{AbstractVector{<:Integer},SVector{M, <:Integer}, NTuple{M, Integer}}) where {M}
+    _non_static_array_warning("dims", dims)
+
+    N = length(dims)
+    N > 0 || throw(DomainError(N, "The length of `dims` must be larger or equal to 1."))
+
+    return ntuple(dim -> HilbertSpace(dims[dim]), Val(N))
+end
+_dims_tuple_of_space(::Nothing) = nothing
+
+_gen_dimensions(dims::AbstractDimensions) = dims
+_gen_dimensions(dims) = ProductDimensions(dims)
+
 # 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)
+function _get_dims_string(dimensions::ProductDimensions)
     dims = dimensions_to_dims(dimensions)
+    isnothing(dimensions.from) && return string(dims)
     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
+space_one_list(dimensions::NTuple{N,AbstractSpace}) where {N} =
+    ntuple(i -> one(dimensions[i]), Val(sum(length, dimensions)))
diff --git a/src/qobj/functions.jl b/src/qobj/functions.jl
index 1d2cf5689..c632c78c0 100644
--- a/src/qobj/functions.jl
+++ b/src/qobj/functions.jl
@@ -187,40 +187,20 @@ julia> a.dims, O.dims
 ```
 """
 function Base.kron(
-    A::AbstractQuantumObject{OpType,<:Dimensions},
-    B::AbstractQuantumObject{OpType,<:Dimensions},
-) where {OpType<:Union{Ket,Bra,Operator}}
+    A::AbstractQuantumObject{ObjType,<:ProductDimensions},
+    B::AbstractQuantumObject{ObjType,<:ProductDimensions},
+) where {ObjType<: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...)))
-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 Base.kron(
-                    A::AbstractQuantumObject{Operator,<:$ADimType},
-                    B::AbstractQuantumObject{Operator,<:$BDimType},
-                )
-                    QType = promote_op_type(A, B)
-                    _lazy_tensor_warning(A.data, B.data)
-                    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
+    A_dims_from = get_dimensions_from(A)
+    B_dims_from = get_dimensions_from(B)
+    AB_dims_to = (get_dimensions_to(A)..., get_dimensions_to(B)...)
+    AB_dims_from = (A.dimensions.from === B.dimensions.from === nothing) ? nothing : (A_dims_from..., B_dims_from...)
+
+    return QType(kron(A.data, B.data), A.type, ProductDimensions(AB_dims_to, AB_dims_from))
 end
 
-# if A and B are different type (must return Operator with GeneralDimensions)
 for AOpType in (:Ket, :Bra, :Operator)
     for BOpType in (:Ket, :Bra, :Operator)
         if (AOpType != BOpType)
@@ -228,14 +208,13 @@ for AOpType in (:Ket, :Bra, :Operator)
                 function Base.kron(A::AbstractQuantumObject{$AOpType}, B::AbstractQuantumObject{$BOpType})
                     QType = promote_op_type(A, B)
                     _lazy_tensor_warning(A.data, B.data)
-                    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)...),
-                        ),
-                    )
+
+                    A_dims_from = get_dimensions_from(A)
+                    B_dims_from = get_dimensions_from(B)
+                    AB_dims_to = (get_dimensions_to(A)..., get_dimensions_to(B)...)
+                    AB_dims_from = (A.dimensions.from === B.dimensions.from === nothing) ? nothing : (A_dims_from..., B_dims_from...)
+
+                    return QType(kron(A.data, B.data), Operator(), ProductDimensions(AB_dims_to, AB_dims_from))
                 end
             end
         end
diff --git a/src/qobj/operators.jl b/src/qobj/operators.jl
index fb0c2041f..a4aa0d893 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{AbstractDimensions,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{AbstractDimensions,AbstractVector{Int},Tuple},
+    distribution::Union{Symbol,Val} = Val(:haar),
+) = rand_unitary(dimensions, makeVal(distribution))
+function rand_unitary(dimensions::Union{AbstractDimensions,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{AbstractDimensions,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{AbstractDimensions,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{AbstractDimensions,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{AbstractDimensions,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..6fd43d31c 100644
--- a/src/qobj/quantum_object.jl
+++ b/src/qobj/quantum_object.jl
@@ -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
+    if isnothing(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(
@@ -86,15 +83,15 @@ function QuantumObject(A::AbstractMatrix{T}; type = nothing, dims = nothing) whe
         )
     end
 
-    if dims isa Nothing
+    if isnothing(dims)
         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])))
+            dims_to = (size(A, 1), )
+            dims_from = (size(A, 1) == size(A, 2)) ? nothing : (size(A, 2), )
+            dims = ProductDimensions(dims_to, dims_from)
         elseif type isa SuperOperator || type isa OperatorBra
-            dims = Dimensions(isqrt(_size[2]))
+            dims = ProductDimensions(isqrt(size(A, 2)))
         end
     end
 
@@ -103,18 +100,18 @@ end
 
 function QuantumObject(A::AbstractVector{T}; type = nothing, dims = nothing) where {T}
     _check_type(type)
-    if type isa Nothing
+
+    if isnothing(type)
         type = Ket() # default type
     elseif !(type isa Ket) && !(type isa OperatorKet)
         throw(ArgumentError("The argument type must be Ket() or OperatorKet() if the input array is a vector."))
     end
 
-    if dims isa Nothing
-        _size = _get_size(A)
+    if isnothing(dims)
         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 +123,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, 1), size(A, 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..0426309e2 100644
--- a/src/qobj/quantum_object_base.jl
+++ b/src/qobj/quantum_object_base.jl
@@ -128,67 +128,54 @@ 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{Ket,Bra,SuperOperator,OperatorBra,OperatorKet}} = 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::Ket, dimensions::Dimensions, m::Int, n::Int)
+function _check_QuantumObject(::Ket, dimensions::AbstractDimensions, m::Int, n::Int)
     (n != 1) && throw(DomainError((m, n), "The size of the array is not compatible with Ket"))
-    (prod(dimensions) != m) && throw(
+    dims_size = hilbert_dimensions_to_size(dimensions)
+    (dims_size[1] != m || !isnothing(dims_size[2])) && 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(::Bra, dimensions::AbstractDimensions, m::Int, n::Int)
     (m != 1) && throw(DomainError((m, n), "The size of the array is not compatible with Bra"))
-    (prod(dimensions) != n) && throw(
+    dims_size = hilbert_dimensions_to_size(dimensions)
+    (dims_size[1] != n || !isnothing(dims_size[2])) && 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(
-        DimensionMismatch(
-            "Operator with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
-        ),
-    )
-    return nothing
-end
-
-function _check_QuantumObject(type::Operator, dimensions::GeneralDimensions, m::Int, n::Int)
+function _check_QuantumObject(::Operator, dimensions::AbstractDimensions, 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)).",
-        ),
-    )
+    dims_size = hilbert_dimensions_to_size(dimensions)
+    if dims_size[1] != m || (!isnothing(dims_size[2]) && dims_size[2] != n) || (isnothing(dims_size[2]) && m != n)
+        throw(
+            DimensionMismatch(
+                "Operator with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
+            ),
+        )
+    end
     return nothing
 end
 
-function _check_QuantumObject(type::SuperOperator, dimensions::Dimensions, 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(
-        DimensionMismatch(
-            "SuperOperator with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
-        ),
-    )
+function _check_QuantumObject(::SuperOperator, dimensions::AbstractDimensions, m::Int, n::Int)
+    ((m == 1) || (n == 1)) && throw(DomainError((m, n), "The size of the array is not compatible with SuperOperator"))
+    dims_size = liouvillian_dimensions_to_size(dimensions)
+    if dims_size[1] != m || (!isnothing(dims_size[2]) && dims_size[2] != n) || (isnothing(dims_size[2]) && m != n)
+        throw(
+            DimensionMismatch(
+                "SuperOperator with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
+            ),
+        )
+    end
     return nothing
 end
 
-function _check_QuantumObject(type::OperatorKet, dimensions::Dimensions, m::Int, n::Int)
+function _check_QuantumObject(::OperatorKet, dimensions::AbstractDimensions, 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(
+    dims_size = liouvillian_dimensions_to_size(dimensions)
+    (dims_size[1] != m || !isnothing(dims_size[2])) && throw(
         DimensionMismatch(
             "OperatorKet with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
         ),
@@ -196,9 +183,10 @@ 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(::OperatorBra, dimensions::AbstractDimensions, 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(
+    dims_size = liouvillian_dimensions_to_size(dimensions)
+    (dims_size[1] != n || !isnothing(dims_size[2])) && throw(
         DimensionMismatch(
             "OperatorBra with dims = $(_get_dims_string(dimensions)) does not fit the array size = $((m, n)).",
         ),
@@ -220,29 +208,10 @@ 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
-get_dimensions_to(
-    A::AbstractQuantumObject{ObjType,<:Dimensions},
-) 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
-get_dimensions_from(
-    A::AbstractQuantumObject{ObjType,<:Dimensions},
-) 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
+get_dimensions_to(A::AbstractQuantumObject) = A.dimensions.to
+get_dimensions_from(A::AbstractQuantumObject) = isnothing(A.dimensions.from) ? A.dimensions.to : A.dimensions.from
+get_dimensions_from(A::AbstractQuantumObject{Ket}) = space_one_list(A.dimensions.to)
+get_dimensions_to(A::AbstractQuantumObject{Bra}) = space_one_list(A.dimensions.to)
 
 # functions for getting Float or Complex element type
 _float_type(A::AbstractQuantumObject) = _float_type(eltype(A))
diff --git a/src/qobj/quantum_object_evo.jl b/src/qobj/quantum_object_evo.jl
index 8e4de483b..247a2c8f0 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 isnothing(dims)
         if type isa Operator
-            dims =
-                (_size[1] == _size[2]) ? Dimensions(_size[1]) :
-                GeneralDimensions(SVector{2}(SVector{1}(_size[1]), SVector{1}(_size[2])))
+            dims_to = (size(data, 1), )
+            dims_from = (size(data, 1) == size(data, 2)) ? nothing : (size(data, 2), )
+            dims = ProductDimensions(dims_to, dims_from)
         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..53a535723 100644
--- a/src/qobj/space.jl
+++ b/src/qobj/space.jl
@@ -2,24 +2,29 @@
 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)
+    function HilbertSpace(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)
+Base.one(::HilbertSpace) = HilbertSpace(1)
+Base.length(::HilbertSpace) = 1
+
+dimensions_to_dims(s::HilbertSpace) = SVector{1,Int}(s.size)
+hilbert_dimensions_to_size(s::HilbertSpace) = s.size
+liouvillian_dimensions_to_size(s::HilbertSpace) = s.size^2
diff --git a/src/qobj/states.jl b/src/qobj/states.jl
index dd79bee8d..dd73b6351 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{Dimensions,AbstractVector{Int}, Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+- `dimensions::Union{AbstractDimensions,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}) =
+zero_ket(dimensions::Union{AbstractDimensions,AbstractVector{Int},Tuple}) =
     QuantumObject(zeros(ComplexF64, prod(dimensions)), Ket(), dimensions)
 
 @doc raw"""
@@ -70,13 +70,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{Dimensions,AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
+- `dimensions::Union{AbstractDimensions,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})
+function rand_ket(dimensions::Union{AbstractDimensions,AbstractVector{Int},Tuple})
     N = prod(dimensions)
     ψ = rand(ComplexF64, N) .- (0.5 + 0.5im)
     return QuantumObject(normalize!(ψ); type = Ket(), dims = dimensions)
@@ -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{AbstractDimensions,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})
+function maximally_mixed_dm(dimensions::Union{AbstractDimensions,AbstractVector{Int},Tuple})
     N = prod(dimensions)
     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{AbstractDimensions,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{AbstractDimensions,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/qobj/superoperators.jl b/src/qobj/superoperators.jl
index ce2bb5901..0c4b63069 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -22,6 +22,7 @@ _sprepost(A, B) = _spre(A) * _spost(B) # for any other input types
 _spre(A::MatrixOperator) = MatrixOperator(_spre(A.A))
 _spre(A::ScaledOperator) = ScaledOperator(A.λ, _spre(A.L))
 _spre(A::AddedOperator) = AddedOperator(map(op -> _spre(op), A.ops))
+_spre(A::ComposedOperator) = ComposedOperator(map(_spre, A.ops), nothing)
 function _spre(A::AbstractSciMLOperator)
     Id = Eye(size(A, 1))
     _lazy_tensor_warning(Id, A)
@@ -31,6 +32,7 @@ end
 _spost(B::MatrixOperator) = MatrixOperator(_spost(B.A))
 _spost(B::ScaledOperator) = ScaledOperator(B.λ, _spost(B.L))
 _spost(B::AddedOperator) = AddedOperator(map(op -> _spost(op), B.ops))
+_spost(B::ComposedOperator) = ComposedOperator(map(_spost, reverse(B.ops)), nothing)
 function _spost(B::AbstractSciMLOperator)
     B_T = transpose(B)
     Id = Eye(size(B, 1))
diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
index 29f84bc0b..b9aee6595 100644
--- a/src/time_evolution/mesolve.jl
+++ b/src/time_evolution/mesolve.jl
@@ -112,7 +112,7 @@ function mesolveProblem(
     else
         to_dense(_complex_float_type(T), mat2vec(ket2dm(ψ0).data))
     end
-    L = L_evo.data
+    L = cache_operator(L_evo.data, ρ0)
 
     kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...)
     kwargs3 = _merge_tstops(kwargs2, isconstant(L), tlist)
diff --git a/src/time_evolution/sesolve.jl b/src/time_evolution/sesolve.jl
index 93e1c82b1..df0b2cd93 100644
--- a/src/time_evolution/sesolve.jl
+++ b/src/time_evolution/sesolve.jl
@@ -80,7 +80,7 @@ function sesolveProblem(
 
     T = Base.promote_eltype(H_evo, ψ0)
     ψ0 = to_dense(_complex_float_type(T), get_data(ψ0)) # Convert it to dense vector with complex element type
-    U = H_evo.data
+    U = cache_operator(H_evo.data, ψ0)
 
     kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...)
     kwargs3 = _merge_tstops(kwargs2, isconstant(U), tlist)
diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl
index cdcaa5542..1560de45f 100644
--- a/src/time_evolution/smesolve.jl
+++ b/src/time_evolution/smesolve.jl
@@ -109,7 +109,7 @@ function smesolveProblem(
 
     sc_ops_evo_data = Tuple(map(get_data ∘ QobjEvo, sc_ops_list))
 
-    K = get_data(L_evo)
+    K = cache_operator(get_data(L_evo), ρ0)
 
     Id_op = IdentityOperator(prod(dims)^2)
     D_l = map(sc_ops_evo_data) do op
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index a47b3763a..df3372858 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -111,7 +111,7 @@ function ssesolveProblem(
         sc_ops_evo_data,
     )
 
-    K = get_data(-1im * QuantumObjectEvolution(H_eff_evo)) + K_l
+    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 = DiffusionOperator(D_l)
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/quantum_objects.jl b/test/core-test/quantum_objects.jl
index 7551824d0..a599099e3 100644
--- a/test/core-test/quantum_objects.jl
+++ b/test/core-test/quantum_objects.jl
@@ -17,7 +17,7 @@
     # DomainError: incompatible between size of array and type
     @testset "DomainError" begin
         a = rand(ComplexF64, 3, 2)
-        for t in [SuperOperator(), Bra(), OperatorBra()]
+        for t in (Bra(), OperatorBra())
             @test_throws DomainError Qobj(a, type = t)
         end
 
@@ -33,17 +33,14 @@
         @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,)))
+        @test_throws DimensionMismatch Qobj(rand(ComplexF64, 2), type = Ket(), dims = ((2,), (1,)))
+        @test_throws DimensionMismatch Qobj(rand(ComplexF64, 1, 2), type = Bra(), dims = ((1,), (2,)))
     end
 
     # unsupported type of dims
     @testset "unsupported dims" begin
-        @test_throws ArgumentError Qobj(rand(2, 2), dims = 2.0)
-        @test_throws ArgumentError Qobj(rand(2, 2), dims = 2.0 + 0.0im)
+        @test_throws MethodError Qobj(rand(2, 2), dims = 2.0)
+        @test_throws MethodError Qobj(rand(2, 2), dims = 2.0 + 0.0im)
         @test_throws DomainError Qobj(rand(2, 2), dims = 0)
         @test_throws DomainError Qobj(rand(2, 2), dims = (2, -2))
         @test_logs (
@@ -367,8 +364,8 @@
             end
 
             UnionType = Union{
-                QuantumObject{Bra,Dimensions{1,Tuple{Space}},Matrix{T}},
-                QuantumObject{Operator,Dimensions{1,Tuple{Space}},Matrix{T}},
+                QuantumObject{Bra,ProductDimensions{1,1,Tuple{HilbertSpace},Nothing},Matrix{T}},
+                QuantumObject{Operator,ProductDimensions{1,1,Tuple{HilbertSpace},Nothing},Matrix{T}},
             }
             a = rand(T, 1, N)
             @inferred UnionType Qobj(a)
@@ -377,8 +374,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,ProductDimensions{1,1,Tuple{HilbertSpace},Tuple{HilbertSpace}},Matrix{T}},
+                QuantumObject{Operator,ProductDimensions{1,1,Tuple{HilbertSpace},Nothing},Matrix{T}},
             }
             a = rand(T, N, N)
             @inferred UnionType Qobj(a)
diff --git a/test/core-test/quantum_objects_evo.jl b/test/core-test/quantum_objects_evo.jl
index 1b596b325..653d5a1ab 100644
--- a/test/core-test/quantum_objects_evo.jl
+++ b/test/core-test/quantum_objects_evo.jl
@@ -216,7 +216,7 @@
             coef2(p, t) * conj(coef3(p, t)) * (spre(a) * spost(X') - 0.5 * spre(X' * a) - 0.5 * spost(X' * a)) + # cross terms
             conj(coef2(p, t)) * coef3(p, t) * (spre(X) * spost(a') - 0.5 * spre(a' * X) - 0.5 * spost(a' * X))   # cross terms
         L_ti = liouvillian(H_ti) + D1_ti + D2_ti
-        L_td = @test_logs (:warn,) (:warn,) liouvillian(H_td, c_ops) # warnings from lazy tensor in `lindblad_dissipator(c_op2)`
+        L_td = @test_logs liouvillian(H_td, c_ops) # Test there are no warnings from lazy tensor
         ρvec = mat2vec(rand_dm(N))
         @test L_td(p, t) ≈ L_ti
         @test iscached(L_td) == false
@@ -237,7 +237,7 @@
 
         coef_wrong1(t) = nothing
         coef_wrong2(p, t::ComplexF64) = nothing
-        @test_logs (:warn,) (:warn,) liouvillian(H_td * H_td) # warnings from lazy tensor
+        @test_logs liouvillian(H_td * H_td) # Check there are no warnings from lazy tensor
         @test_throws ArgumentError QobjEvo(a, coef_wrong1)
         @test_throws ArgumentError QobjEvo(a, coef_wrong2)
         @test_throws MethodError QobjEvo([[a, coef1], a' * a, [a', coef2]])
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index 0418f45ed..d918b52d1 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -1256,3 +1256,26 @@ end
     @test all(isapprox.(Z_analytic, sol_se.expect[2, :]; atol = 1e-6))
     @test all(isapprox.(Z_analytic, sol_me.expect[2, :]; atol = 1e-6))
 end
+
+@testitem "ComposedOperator support in time evolution" begin
+    N = 10  # number of basis states
+    a = destroy(N)
+    H = a' * a
+
+    ψ0 = basis(N, 9)  # initial state
+
+    γ0 = 0.5
+    γ1(p, t) = sqrt(p[1] * exp(-t))
+    c_ops_1 = [QobjEvo(2a, γ1)] # This generates a ScaledOperator internally
+    c_ops_2 = [QobjEvo(a, γ1) + QobjEvo(a, γ1)] # This generates a ComposedOperator internally
+
+    e_ops = [a' * a]
+
+    tlist = range(0, 10, 100)
+    params = [γ0]
+    sol_1 = mesolve(H, ψ0, tlist, c_ops_1, e_ops = e_ops; progress_bar = Val(false), params = params, saveat = tlist)
+    sol_2 = mesolve(H, ψ0, tlist, c_ops_2, e_ops = e_ops; progress_bar = Val(false), params = params, saveat = tlist)
+
+    @test sol_1.expect[1, :] ≈ sol_2.expect[1, :] atol=1e-10
+    @test all(x -> isapprox(x[1].data, x[2].data; atol = 1e-10), zip(sol_1.states, sol_2.states))
+end