From 842e668ed2e7e3f6ffdb34cb3b50c3fc9fdd4f0c Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Sun, 15 Sep 2024 20:58:53 +0800
Subject: [PATCH 01/10] update README

---
 README.md | 2 --
 1 file changed, 2 deletions(-)

diff --git a/README.md b/README.md
index 369f9cec..cd4219aa 100644
--- a/README.md
+++ b/README.md
@@ -144,8 +144,6 @@ We can extract the expectation value of the number operator $\hat{a}^\dagger \ha
 
 We can easily pass the computation to the GPU, by simply passing all the `Qobj`s to the GPU:
 
-> **_NOTE:_** The described feature requires `Julia 1.9+`.
-
 ```julia
 using QuantumToolbox
 using CUDA

From 871a29aecc67c41ef50a818939679636af7866ef Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Sun, 15 Sep 2024 20:59:37 +0800
Subject: [PATCH 02/10] extended methods for `expect` and `variance`

---
 src/qobj/functions.jl      | 19 +++++++++++++++++--
 src/qobj/superoperators.jl |  3 +++
 test/quantum_objects.jl    |  7 +++++++
 3 files changed, 27 insertions(+), 2 deletions(-)

diff --git a/src/qobj/functions.jl b/src/qobj/functions.jl
index 667f7930..58bfd4f7 100644
--- a/src/qobj/functions.jl
+++ b/src/qobj/functions.jl
@@ -17,7 +17,7 @@ ket2dm(ψ::QuantumObject{<:AbstractArray{T},KetQuantumObject}) where {T} = ψ *
 ket2dm(ρ::QuantumObject{<:AbstractArray{T},OperatorQuantumObject}) where {T} = ρ
 
 @doc raw"""
-    expect(O::QuantumObject, ψ::QuantumObject)
+    expect(O::QuantumObject, ψ::Union{QuantumObject,Vector{QuantumObject}})
 
 Expectation value of the [`Operator`](@ref) `O` with the state `ψ`. The state can be a [`Ket`](@ref), [`Bra`](@ref) or [`Operator`](@ref).
 
@@ -27,6 +27,8 @@ If `ψ` is a density matrix ([`Operator`](@ref)), the function calculates ``\tex
 
 The function returns a real number if `O` is of `Hermitian` type or `Symmetric` type, and returns a complex number otherwise. You can make an operator `O` hermitian by using `Hermitian(O)`.
 
+Note that `ψ` can also be given as a list of [`QuantumObject`](@ref), it returns a list of expectation values.
+
 # Examples
 
 ```
@@ -77,20 +79,33 @@ function expect(
 ) where {TF<:Number,TR<:Real,T2}
     return real(tr(O * ρ))
 end
+function expect(
+    O::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
+    ρ::Vector{<:QuantumObject},
+) where {T1}
+    _expect = _ρ -> expect(O, _ρ)
+    return _expect.(ρ)
+end
 
 @doc raw"""
-    variance(O::QuantumObject, ψ::QuantumObject)
+    variance(O::QuantumObject, ψ::Union{QuantumObject,Vector{QuantumObject}})
 
 Variance of the [`Operator`](@ref) `O`: ``\langle\hat{O}^2\rangle - \langle\hat{O}\rangle^2``,
 
 where ``\langle\hat{O}\rangle`` is the expectation value of `O` with the state `ψ` (see also [`expect`](@ref)), and the state `ψ` can be a [`Ket`](@ref), [`Bra`](@ref) or [`Operator`](@ref).
 
 The function returns a real number if `O` is hermitian, and returns a complex number otherwise.
+
+Note that `ψ` can also be given as a list of [`QuantumObject`](@ref), it returns a list of expectation values.
 """
 variance(
     O::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
     ψ::QuantumObject{<:AbstractArray{T2}},
 ) where {T1,T2} = expect(O^2, ψ) - expect(O, ψ)^2
+variance(
+    O::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
+    ψ::Vector{<:QuantumObject},
+) where {T1} = expect(O^2, ψ) .- expect(O, ψ).^2
 
 @doc raw"""
     sparse_to_dense(A::QuantumObject)
diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl
index 61945013..62d761ee 100644
--- a/src/qobj/superoperators.jl
+++ b/src/qobj/superoperators.jl
@@ -25,6 +25,7 @@ Since the density matrix is vectorized in [`OperatorKet`](@ref) form: ``|\hat{\r
 ```math
 \mathcal{O} \left(\hat{A}\right) \left[ \hat{\rho} \right] = \hat{\mathbb{1}} \otimes \hat{A} ~ |\hat{\rho}\rangle\rangle
 ```
+(see the section in documentation: [Superoperators and Vectorized Operators](@ref doc:Superoperators-and-Vectorized-Operators) for more details)
 
 The optional argument `Id_cache` can be used to pass a precomputed identity matrix. This can be useful when
 the same function is applied multiple times with a known Hilbert space dimension.
@@ -42,6 +43,7 @@ Since the density matrix is vectorized in [`OperatorKet`](@ref) form: ``|\hat{\r
 ```math
 \mathcal{O} \left(\hat{B}\right) \left[ \hat{\rho} \right] = \hat{B}^T \otimes \hat{\mathbb{1}} ~ |\hat{\rho}\rangle\rangle
 ```
+(see the section in documentation: [Superoperators and Vectorized Operators](@ref doc:Superoperators-and-Vectorized-Operators) for more details)
 
 The optional argument `Id_cache` can be used to pass a precomputed identity matrix. This can be useful when
 the same function is applied multiple times with a known Hilbert space dimension.
@@ -59,6 +61,7 @@ Since the density matrix is vectorized in [`OperatorKet`](@ref) form: ``|\hat{\r
 ```math
 \mathcal{O} \left(\hat{A}, \hat{B}\right) \left[ \hat{\rho} \right] = \hat{B}^T \otimes \hat{A} ~ |\hat{\rho}\rangle\rangle = \textrm{spre}(A) * \textrm{spost}(B) ~ |\hat{\rho}\rangle\rangle
 ```
+(see the section in documentation: [Superoperators and Vectorized Operators](@ref doc:Superoperators-and-Vectorized-Operators) for more details)
 
 See also [`spre`](@ref) and [`spost`](@ref).
 """
diff --git a/test/quantum_objects.jl b/test/quantum_objects.jl
index e61ff996..cbf54667 100644
--- a/test/quantum_objects.jl
+++ b/test/quantum_objects.jl
@@ -384,6 +384,11 @@
         @test expect(a, ρ) ≈ tr(a * ρ)
         @test variance(a, ρ) ≈ tr(a^2 * ρ) - tr(a * ρ)^2
 
+        # when input is a vector of states
+        xlist = [1.0, 1.0im, -1.0, -1.0im]
+        ψlist = [normalize!(basis(N, 4) + x * basis(N, 3)) for x in xlist]
+        @test all(expect(a', ψlist) .≈ xlist)
+
         @testset "Type Inference (expect)" begin
             @inferred expect(a, ψ)
             @inferred expect(a, ψ')
@@ -391,6 +396,8 @@
             @inferred variance(a, ψ')
             @inferred expect(a, ρ)
             @inferred variance(a, ρ)
+            @inferred expect(a, ψlist)
+            @inferred variance(a, ψlist)
         end
     end
 

From f58de01caa4186f5c975f9519ed24ad81546e822 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Sun, 15 Sep 2024 21:00:21 +0800
Subject: [PATCH 03/10] add new section to  docs

---
 docs/make.jl                                  |   2 +-
 docs/src/index.md                             |   3 -
 docs/src/{users_guide => }/type_stability.md  |   0
 .../QuantumObject/QuantumObject.md            |   5 +-
 docs/src/users_guide/extensions/cuda.md       |   3 -
 docs/src/users_guide/states_and_operators.md  | 406 +++++++++++++++++-
 docs/src/users_guide/tensor.md                |   2 +-
 7 files changed, 404 insertions(+), 17 deletions(-)
 rename docs/src/{users_guide => }/type_stability.md (100%)

diff --git a/docs/make.jl b/docs/make.jl
index fccd116d..9255bca5 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -21,6 +21,7 @@ const PAGES = [
     "Getting Started" => [
         "Introduction" => "index.md",
         "Key differences from QuTiP" => "qutip_differences.md",
+        "The Importance of Type-Stability" => "type_stability.md",
         # "Cite QuantumToolbox.jl" => "cite.md",
     ],
     "Users Guide" => [
@@ -28,7 +29,6 @@ const PAGES = [
             "users_guide/QuantumObject/QuantumObject.md",
             "users_guide/QuantumObject/QuantumObject_functions.md",
         ],
-        "The Importance of Type-Stability" => "users_guide/type_stability.md",
         "Manipulating States and Operators" => "users_guide/states_and_operators.md",
         "Tensor Products and Partial Traces" => "users_guide/tensor.md",
         "Time Evolution and Dynamics" => [
diff --git a/docs/src/index.md b/docs/src/index.md
index 149df299..4d1642c6 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -97,9 +97,6 @@ We can extract the expectation value of the number operator ``\hat{a}^\dagger \h
 
 We can easily pass the computation to the GPU, by simply passing all the `Qobj`s to the GPU:
 
-!!! compat "Compat"
-    The described feature requires `Julia 1.9+`. See [CUDA extension](@ref doc:CUDA) for more details.
-
 ```julia
 using QuantumToolbox
 using CUDA
diff --git a/docs/src/users_guide/type_stability.md b/docs/src/type_stability.md
similarity index 100%
rename from docs/src/users_guide/type_stability.md
rename to docs/src/type_stability.md
diff --git a/docs/src/users_guide/QuantumObject/QuantumObject.md b/docs/src/users_guide/QuantumObject/QuantumObject.md
index 61e9569f..01a86914 100644
--- a/docs/src/users_guide/QuantumObject/QuantumObject.md
+++ b/docs/src/users_guide/QuantumObject/QuantumObject.md
@@ -14,6 +14,9 @@ The key difference between classical and quantum mechanics is the use of operato
 - `CUDA.CUSPARSE.CuSparseMatrixCSR` (sparse GPU matrix)
 - and even more ...
 
+!!! note "Support for GPU arrays"
+    See [CUDA extension](@ref doc:CUDA) for more details.
+
 We can create a [`QuantumObject`](@ref) with a user defined data set by passing an array of data into the [`QuantumObject`](@ref):
 
 ```@setup Qobj
@@ -53,7 +56,7 @@ Qobj(rand(4, 4), type = SuperOperator)
 ```
 
 !!! note "Difference between `dims` and `size`"
-    Notice that `type`, `dims`, and `size` will change according to the input `data`. Although `dims` and `size` appear to be the same, `dims` keep tracking the dimension of individual Hilbert spaces of a multipartite system, while `size` does not. We refer the reader to the section [tensor products and partial traces](@ref doc:Tensor-products) for more information.
+    Notice that `type`, `dims`, and `size` will change according to the input `data`. Although `dims` and `size` appear to be the same, `dims` keep tracking the dimension of individual Hilbert spaces of a multipartite system, while `size` does not. We refer the reader to the section [Tensor Products and Partial Traces](@ref doc:Tensor-products-and-Partial-Traces) for more information.
 
 ## States and operators
 
diff --git a/docs/src/users_guide/extensions/cuda.md b/docs/src/users_guide/extensions/cuda.md
index e40c9091..11290c6f 100644
--- a/docs/src/users_guide/extensions/cuda.md
+++ b/docs/src/users_guide/extensions/cuda.md
@@ -4,9 +4,6 @@
 
 This is an extension to support `QuantumObject.data` conversion from standard dense and sparse CPU arrays to GPU ([`CUDA.jl`](https://github.com/JuliaGPU/CUDA.jl)) arrays.
 
-!!! note "Requirements"
-    The [`CUDA.jl`](https://github.com/JuliaGPU/CUDA.jl) extension for `QuantumToolbox.jl` requires `Julia 1.9+`.
-
 This extension will be automatically loaded if user imports both `QuantumToolbox` and [`CUDA.jl`](https://github.com/JuliaGPU/CUDA.jl):
 
 ```julia
diff --git a/docs/src/users_guide/states_and_operators.md b/docs/src/users_guide/states_and_operators.md
index 565aece2..afed566d 100644
--- a/docs/src/users_guide/states_and_operators.md
+++ b/docs/src/users_guide/states_and_operators.md
@@ -1,17 +1,407 @@
 # [States and Operators](@id doc:States-and-Operators)
 
-This page is still under construction, please visit [API](@ref doc-API) first.
-
 ## Introduction
+In the previous guide section [Basic Operations on Quantum Objects](@ref doc:Qobj), we saw how to create states and operators, using the functions built into `QuantumToolbox`. In this portion of the guide, we will look at performing basic operations with states and operators. For more detailed demonstrations on how to use and manipulate these objects, see the examples given in the tutorial section.
+
+```@setup states_and_operators
+using QuantumToolbox
+```
+
+## [State Vectors (kets or bras)](@id doc:State-vectors)
+Here we begin by creating a Fock [`basis`](@ref) (or [`fock`](@ref)) vacuum state vector ``|0\rangle`` with in a Hilbert space with `5` number states, from `0` to `4`:
+
+```@example states_and_operators
+vac = basis(5, 0)
+```
+
+and then create a lowering operator ``\hat{a}`` corresponding to `5` number states using the [`destroy`](@ref) function:
+
+```@example states_and_operators
+a = destroy(5)
+```
+
+Now lets apply the lowering operator `a` to our vacuum state `vac`:
+
+```@example states_and_operators
+a * vac
+```
+
+We see that, as expected, the vacuum is transformed to the zero vector. A more interesting example comes from using the `adjoint` of the lowering operator ``\hat{a}``, the raising operator ``\hat{a}^\dagger``:
+
+```@example states_and_operators
+a' * vac
+```
+
+The raising operator has in indeed raised the state `vac` from the vacuum to the ``|1\rangle`` state. Instead of using the `adjoint` method to raise the state, we could have also used the built-in [`create`](@ref) function to make a raising operator:
+
+```@example states_and_operators
+ad = create(5)
+ad * vac
+```
+
+which does the same thing. We can raise the vacuum state more than once by successively apply the raising operator:
+
+```@example states_and_operators
+ad * ad * vac
+```
+
+or just taking the square of the raising operator ``\left(\hat{a}^\dagger\right)^2``:
+
+```@example states_and_operators
+ad^2 * vac
+```
+
+Applying the raising operator twice gives the expected ``\sqrt{n+1}`` dependence. We can use the product of ``a^\dagger a`` to also apply the number operator to the state vector `vac`:
+
+```@example states_and_operators
+ad * a * vac
+```
+
+or on the ``|1\rangle`` state:
+
+```@example states_and_operators
+ad * a * (ad * vac)
+```
+
+or on the ``|2\rangle`` state:
+
+```@example states_and_operators
+ad * a * (ad^2 * vac)
+```
+
+Notice how in this last example, application of the number operator does not give the expected value ``n=2``, but rather ``2\sqrt{2}``. This is because this last state is not normalized to unity as ``a^\dagger|n\rangle=\sqrt{n+1}|n+1\rangle``. Therefore, we should [`normalize`](@ref) (or use [`unit`](@ref)) our vector first:
+
+```@example states_and_operators
+ad * a * normalize(ad^2 * vac)
+```
+
+Since we are giving a demonstration of using states and operators, we have done a lot more work than we should have. For example, we do not need to operate on the vacuum state to generate a higher number Fock state. Instead we can use the [`basis`](@ref) (or [`fock`](@ref)) function to directly obtain the required state:
+
+```@example states_and_operators
+ket = basis(5, 2)
+```
+
+Notice how it is automatically normalized. We can also use the built in number operator [`num`](@ref):
+
+```@example states_and_operators
+n = num(5)
+```
+
+Therefore, instead of `ad * a * normalize(ad^2 * vac)`, we have:
+
+```@example states_and_operators
+n * ket
+```
+
+We can also create superpositions of states:
+
+```@example states_and_operators
+ket = normalize(basis(5, 0) + basis(5, 1))
+```
+
+where we have used the `normalize` function again to normalize the state. Apply the number opeartor again:
+
+```@example states_and_operators
+n * ket
+```
+
+We can also create coherent states and squeezed states by applying the [`displace`](@ref) and [`squeeze`](@ref) functions to the vacuum state:
+
+```@example states_and_operators
+vac = basis(5, 0)
+
+d = displace(5, 1im)
+
+s = squeeze(5, 0.25 + 0.25im)
+
+d * vac
+```
+
+```@example states_and_operators
+d * s * vac
+```
+
+Of course, displacing the vacuum gives a coherent state, which can also be generated using the built in [`coherent`](@ref) function.
+
+## [Density matrices](@id doc:Density-matrices)
+
+One of the main purpose of `QuantumToolbox` is to explore the dynamics of open quantum systems, where the most general state of a system is no longer a state vector, but rather a density matrix. Since operations on density matrices operate identically to those of vectors, we will just briefly highlight creating and using these structures.
+
+The simplest density matrix is created by forming the outer-product ``|\psi\rangle\langle\psi|`` of a ket vector:
+
+```@example states_and_operators
+ket = basis(5, 2)
+ket * ket'
+```
+
+A similar task can also be accomplished via the [`fock_dm`](@ref) or [`ket2dm`](@ref) functions:
+
+```@example states_and_operators
+fock_dm(5, 2)
+```
+
+```@example states_and_operators
+ket2dm(ket)
+```
+
+If we want to create a density matrix with equal classical probability of being found in the ``|2\rangle`` or ``|4\rangle`` number states, we can do the following:
+
+```@example states_and_operators
+0.5 * fock_dm(5, 2) + 0.5 * fock_dm(5, 4) # with fock_dm
+0.5 * ket2dm(basis(5, 2)) + 0.5 * ket2dm(basis(5, 4)) # with ket2dm
+```
+
+There are also several other built-in functions for creating predefined density matrices, for example [`coherent_dm`](@ref) and [`thermal_dm`](@ref) which create coherent state and thermal state density matrices, respectively.
+
+```@example states_and_operators
+coherent_dm(5, 1.25)
+```
+
+```@example states_and_operators
+thermal_dm(5, 1.25)
+```
+
+`QuantumToolbox` also provides a set of distance metrics for determining how close two density matrix distributions are to each other. Included are the [`fidelity`](@ref), and trace distance ([`tracedist`](@ref)).
+
+```@example states_and_operators
+x = coherent_dm(5, 1.25)
+
+y = coherent_dm(5, 1.25im)
+
+z = thermal_dm(5, 0.125)
+
+fidelity(x, y)
+```
+Note that the definition of [`fidelity`](@ref) here is from **Nielsen & Chuang, "Quantum Computation and Quantum Information"**. It is the square root of the fidelity defined in **R. Jozsa, Journal of Modern Optics, 41:12, 2315 (1994)**. We also know that for two pure states, the trace distance (``T``) and the fidelity (``F``) are related by ``T = \sqrt{1-F^2}``:
+
+```@example states_and_operators
+tracedist(x, y) ≈ sqrt(1 - (fidelity(x, y))^2)
+```
+
+For a pure state and a mixed state, ``1 - F \leq T`` which can also be verified:
+
+```@example states_and_operators
+1 - fidelity(x, z) < tracedist(x, z)
+```
+
+## [Two-level systems (Qubits)](@id doc:Two-level-systems)
+
+Having spent a fair amount of time on basis states that represent harmonic oscillator states, we now move on to qubit, or two-level quantum systems (for example a spin-``1/2``). To create a state vector corresponding to a qubit system, we use the same basis, or fock, function with only two levels:
+
+```@example states_and_operators
+spin = basis(2, 0)
+```
+
+Now at this point one may ask how this state is different than that of a harmonic oscillator in the vacuum state truncated to two energy levels?
+
+```@example states_and_operators
+vac = basis(2, 0)
+```
+
+At this stage, there is no difference. This should not be surprising as we called the exact same function twice. The difference between the two comes from the action of the spin operators [`sigmax`](@ref), [`sigmay`](@ref), [`sigmaz`](@ref), [`sigmap`](@ref), and [`sigmam`](@ref) on these two-level states. For example, if `vac` corresponds to the vacuum state of a harmonic oscillator, then, as we have already seen, we can use the raising operator ([`create`](@ref)) to get the ``|1\rangle`` state:
+
+```@example states_and_operators
+create(2) * vac
+```
+
+For a spin system, the operator analogous to the raising operator is the ``\sigma_+`` operator [`sigmap`](@ref). Applying on the spin state gives:
+
+```@example states_and_operators
+sigmap() * spin
+```
+
+Now we see the difference! The [`sigmap`](@ref) operator acting on the spin state returns the zero vector. Why is this? To see what happened, let us use the [`sigmaz`](@ref) operator:
+
+```@example states_and_operators
+sigmaz()
+```
+
+```@example states_and_operators
+sigmaz() * spin
+```
+
+```@example states_and_operators
+spin2 = basis(2, 1)
+```
+
+```@example states_and_operators
+sigmaz() * spin2
+```
+
+The answer is now apparent. Since the `QuantumToolbox` [`sigmaz`](@ref) function uses the standard ``Z``-basis representation of the ``\sigma_z`` spin operator, the `spin` state corresponds to the ``|\uparrow\rangle`` state of a two-level spin system while `spin2` gives the ``|\downarrow\rangle`` state. Therefore, in our previous example `sigmap() * spin`, we raised the qubit state out of the truncated two-level Hilbert space resulting in the zero state.
+
+While at first glance this convention might seem somewhat odd, it is in fact quite handy. For one, the spin operators remain in the conventional form. Second, this corresponds nicely with the quantum information definitions of qubit states, where the excited ``|\uparrow\rangle`` state is label as ``|0\rangle``, and the ``|\downarrow\rangle`` state by ``|1\rangle``.
+
+If one wants to create spin operators for higher spin systems, then the [`jmat`](@ref) function comes in handy.
+
+## [Expectation values](@id doc:Expectation-values)
+
+Some of the most important information about quantum systems comes from calculating the expectation value of operators, both Hermitian and non-Hermitian, as the state or density matrix of the system varies in time. Therefore, in this section we demonstrate the use of the [`expect`](@ref) function. To begin:
+
+```@example states_and_operators
+vac = basis(5, 0)
+
+one = basis(5, 1)
+
+c = create(5)
+
+N = num(5)
+
+coh = coherent_dm(5, 1.0im)
+
+cat = normalize(basis(5, 4) + 1.0im * basis(5, 3))
+
+println(expect(N, vac) ≈ 0)
+println(expect(N, one) ≈ 1)
+println(expect(N, coh) ≈ 0.9970555745806597)
+println(expect(c, cat) ≈ 1im)
+```
+
+The [`expect`](@ref) function also accepts lists or arrays of state vectors or density matrices for the second input:
+
+```@example states_and_operators
+states = [normalize(c^k * vac) for k in 0:4]
+
+expect(N, states)
+```
+
+```@example states_and_operators
+cat_list = [normalize(basis(5, 4) + x * basis(5, 3)) for x in [0, 1.0im, -1.0, -1.0im]]
+
+expect(c, cat_list)
+```
+
+Notice how in this last example, all of the return values are complex numbers. This is because the expect function looks to see whether the operator is Hermitian or not. If the operator is Hermitian, then the output will always be real. In the case of non-Hermitian operators, the return values may be complex. Therefore, the expect function will return an array of complex values for non-Hermitian operators when the input is a list/array of states or density matrices.
+
+Of course, the expect function works for spin states and operators:
+
+```@example states_and_operators
+up = basis(2, 0)
+
+dn = basis(2, 1)
+
+println(expect(sigmaz(), up) ≈ 1)
+println(expect(sigmaz(), dn) ≈ -1)
+```
+
+as well as the composite objects discussed in the next section [Tensor Products and Partial Traces](@ref doc:Tensor-products-and-Partial-Traces):
+
+```@example states_and_operators
+spin1 = basis(2, 0)
+
+spin2 = basis(2, 1)
+
+two_spins = tensor(spin1, spin2)
+
+sz1 = tensor(sigmaz(), qeye(2))
+
+sz2 = tensor(qeye(2), sigmaz())
+
+println(expect(sz1, two_spins) ≈ 1)
+println(expect(sz2, two_spins) ≈ -1)
+```
+
+## [Superoperators and Vectorized Operators](@id doc:Superoperators-and-Vectorized-Operators)
+
+In addition to state vectors and density operators, `QuantumToolbox` allows for representing maps that act linearly on density operators using the Liouville supermatrix formalisms.
+
+This support is based on the correspondence between linear operators acting on a Hilbert space, and vectors in two copies of that Hilbert space (which is also called the Fock-Liouville space), 
+```math
+\textrm{vec} : \mathcal{L}(\mathcal{H}) \rightarrow \mathcal{H}\otimes\mathcal{H}.
+```
+Therefore, a given density matrix ``\rho`` can then be vectorized, denoted as 
+```math
+|\rho\rangle\rangle = \textrm{vec}(\rho).
+```
+
+`QuantumToolbox` uses the column-stacking convention for the isomorphism between ``\mathcal{L}(\mathcal{H})`` and ``\mathcal{H}\otimes\mathcal{H}``. This isomorphism is implemented by the functions [`mat2vec`](@ref) and [`vec2mat`](@ref):
+
+```@example states_and_operators
+rho = Qobj([1 2; 3 4])
+```
+
+```@example states_and_operators
+vec_rho = mat2vec(rho)
+```
+
+```@example states_and_operators
+rho2 = vec2mat(vec_rho)
+```
+
+The `QuantumObject.type` attribute indicates whether a quantum object is a vector corresponding to an [`OperatorKet`](@ref), or its Hermitian conjugate [`OperatorBra`](@ref). One can also use [`isoperket`](@ref) and [`isoperbra`](@ref) to check the type:
+
+```@example states_and_operators
+println(isoper(vec_rho))
+println(isoperket(vec_rho))
+println(isoperbra(vec_rho))
+println(isoper(vec_rho'))
+println(isoperket(vec_rho'))
+println(isoperbra(vec_rho'))
+```
+
+Because `Julia` is a column-oriented languages (like `Fortran` and `MATLAB`), in `QuantumToolbox`, we define the [`spre`](@ref) (left), [`spost`](@ref) (right), and [`sprepost`](@ref) (left-and-right) multiplication superoperators as follows:
+
+```math
+\begin{align}
+A\rho~~~ &\rightarrow \textrm{spre}(A) * \textrm{vec}(\rho) = \mathbb{1}\otimes A ~ |\rho\rangle\rangle,\notag\\
+\rho B &\rightarrow \textrm{spost}(B) * \textrm{vec}(\rho) = B^T\otimes \mathbb{1} ~ |\rho\rangle\rangle,\notag\\
+A \rho B &\rightarrow \textrm{sprepost}(A,B) * \textrm{vec}(\rho) = B^T\otimes A ~ |\rho\rangle\rangle,\notag
+\end{align}
+```
+where ``\mathbb{1}`` represents the identity operator with Hilbert space dimension equal to ``\rho``.
+
+```@example states_and_operators
+A = Qobj([1 2; 3 4])
+S_A = spre(A)
+```
+
+```@example states_and_operators
+B = Qobj([5 6; 7 8])
+S_B = spost(B)
+```
+
+```@example states_and_operators
+S_AB = sprepost(A, B)
+```
+
+```@example states_and_operators
+S_AB ≈ S_A * S_B ≈ S_B * S_A
+```
+
+One can also use [`issuper`](@ref) to check the type:
+
+```@example states_and_operators
+println(isoper(S_AB))
+println(issuper(S_AB))
+```
+
+With the above definitions, the following equality holds in `Julia`:
+
+```math
+\textrm{vec}(A \rho B) = \textrm{spre}(A) * \textrm{spre}(B) * \textrm{vec}(\rho) = \textrm{sprepost}(A,B) * \textrm{vec}(\rho) ~~\forall~~A, B, \rho
+```
+
+```@example states_and_operators
+N  = 10
+A = Qobj(rand(ComplexF64, N, N))
+B = Qobj(rand(ComplexF64, N, N))
+ρ = rand_dm(N) # random density matrix
+mat2vec(A * ρ * B) ≈ spre(A) * spost(B) * mat2vec(ρ) ≈ sprepost(A, B) * mat2vec(ρ)
+```
 
-## [State Vectors (kets or bras)](@id doc: State vectors)
+In addition, dynamical generators on this extended space, often called Liouvillian superoperators, can be created using the [`liouvillian`](@ref) function. Each of these takes a Hamiltonian along with a list of collapse operators, and returns a [`type=SuperOperator`](@ref SuperOperator) object that can be exponentiated to find the superoperator for that evolution.
 
-## [Density matrices](@id doc: Density matrices)
+```@example states_and_operators
+H = 10 * sigmaz()
 
-## [Two-level systems (qubits)](@id doc: Two-level systems)
+c = destroy(2)
 
-## [Expectation values](@id doc: Expectation values)
+L = liouvillian(H, [c])
+```
 
-## [Superoperators and Vectorized Operators](@id doc: Superoperators and Vectorized Operators)
+```@example states_and_operators
+t = 0.8
+exp(L * t)
+```
 
-## [Generating Random States and Operators](@id doc: Generating Random States and Operators)
+See the section [Time Evolution and Quantum System Dynamics](@ref doc:Time-Evolution-and-Quantum-System-Dynamics) for more details.
diff --git a/docs/src/users_guide/tensor.md b/docs/src/users_guide/tensor.md
index 23d89e0a..faf4ae7d 100644
--- a/docs/src/users_guide/tensor.md
+++ b/docs/src/users_guide/tensor.md
@@ -1,3 +1,3 @@
-# [Tensor products](@id doc:Tensor-products)
+# [Tensor Products and Partial Traces](@id doc:Tensor-products-and-Partial-Traces)
 
 This page is still under construction, please visit [API](@ref doc-API) first.
\ No newline at end of file

From eb4a3490ba19f2515a1e8f0407ddd2eb2f1baa3f Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Sun, 15 Sep 2024 21:03:09 +0800
Subject: [PATCH 04/10] format files

---
 src/qobj/functions.jl | 11 +++--------
 1 file changed, 3 insertions(+), 8 deletions(-)

diff --git a/src/qobj/functions.jl b/src/qobj/functions.jl
index 58bfd4f7..1c506aa2 100644
--- a/src/qobj/functions.jl
+++ b/src/qobj/functions.jl
@@ -79,10 +79,7 @@ function expect(
 ) where {TF<:Number,TR<:Real,T2}
     return real(tr(O * ρ))
 end
-function expect(
-    O::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
-    ρ::Vector{<:QuantumObject},
-) where {T1}
+function expect(O::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject}, ρ::Vector{<:QuantumObject}) where {T1}
     _expect = _ρ -> expect(O, _ρ)
     return _expect.(ρ)
 end
@@ -102,10 +99,8 @@ variance(
     O::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
     ψ::QuantumObject{<:AbstractArray{T2}},
 ) where {T1,T2} = expect(O^2, ψ) - expect(O, ψ)^2
-variance(
-    O::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
-    ψ::Vector{<:QuantumObject},
-) where {T1} = expect(O^2, ψ) .- expect(O, ψ).^2
+variance(O::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject}, ψ::Vector{<:QuantumObject}) where {T1} =
+    expect(O^2, ψ) .- expect(O, ψ) .^ 2
 
 @doc raw"""
     sparse_to_dense(A::QuantumObject)

From 6fb83bca142a19856c3dcebfc7c3d14c19035adc Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Sun, 15 Sep 2024 21:05:42 +0800
Subject: [PATCH 05/10] change section title in docs

---
 docs/src/users_guide/states_and_operators.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/users_guide/states_and_operators.md b/docs/src/users_guide/states_and_operators.md
index afed566d..5849bf07 100644
--- a/docs/src/users_guide/states_and_operators.md
+++ b/docs/src/users_guide/states_and_operators.md
@@ -1,4 +1,4 @@
-# [States and Operators](@id doc:States-and-Operators)
+# [Manipulating States and Operators](@id doc:Manipulating-States-and-Operators)
 
 ## Introduction
 In the previous guide section [Basic Operations on Quantum Objects](@ref doc:Qobj), we saw how to create states and operators, using the functions built into `QuantumToolbox`. In this portion of the guide, we will look at performing basic operations with states and operators. For more detailed demonstrations on how to use and manipulate these objects, see the examples given in the tutorial section.

From 02658fdbea8eebea08c6c5ce1d0f3f0a2c552554 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Mon, 16 Sep 2024 00:01:23 +0800
Subject: [PATCH 06/10] extended methods for `tensor` and `kron`

---
 src/qobj/functions.jl   | 10 +++++++---
 src/qobj/synonyms.jl    |  2 +-
 test/quantum_objects.jl | 13 +++++++++++++
 3 files changed, 21 insertions(+), 4 deletions(-)

diff --git a/src/qobj/functions.jl b/src/qobj/functions.jl
index 1c506aa2..a69d8c10 100644
--- a/src/qobj/functions.jl
+++ b/src/qobj/functions.jl
@@ -187,11 +187,15 @@ Quantum Object:   type=Operator   dims=[20, 20]   size=(400, 400)   ishermitian=
 ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠦
 ```
 """
-function LinearAlgebra.kron(
+LinearAlgebra.kron(
     A::QuantumObject{<:AbstractArray{T1},OpType},
     B::QuantumObject{<:AbstractArray{T2},OpType},
-) where {T1,T2,OpType<:Union{KetQuantumObject,BraQuantumObject,OperatorQuantumObject}}
-    return QuantumObject(kron(A.data, B.data), A.type, vcat(A.dims, B.dims))
+) where {T1,T2,OpType<:Union{KetQuantumObject,BraQuantumObject,OperatorQuantumObject}} =
+    QuantumObject(kron(A.data, B.data), A.type, vcat(A.dims, B.dims))
+LinearAlgebra.kron(A::QuantumObject) = A
+function LinearAlgebra.kron(A::Vector{<:QuantumObject})
+    @warn "`tensor(A)` or `kron(A)` with `A` is a `Vector` can hurt performance. Try to use `tensor(A...)` or `kron(A...)` instead."
+    return kron(A...)
 end
 
 @doc raw"""
diff --git a/src/qobj/synonyms.jl b/src/qobj/synonyms.jl
index 5b991d5c..9457bbb6 100644
--- a/src/qobj/synonyms.jl
+++ b/src/qobj/synonyms.jl
@@ -184,7 +184,7 @@ Quantum Object:   type=Operator   dims=[2, 2, 2]   size=(8, 8)   ishermitian=tru
  1.0+0.0im      ⋅          ⋅             ⋅          ⋅          ⋅
 ```
 """
-tensor(A::QuantumObject...) = kron(A...)
+tensor(A...) = kron(A...)
 
 @doc raw"""
     ⊗(A::QuantumObject, B::QuantumObject)
diff --git a/test/quantum_objects.jl b/test/quantum_objects.jl
index cbf54667..64d1433f 100644
--- a/test/quantum_objects.jl
+++ b/test/quantum_objects.jl
@@ -309,11 +309,24 @@
             @inferred a .^ 2
             @inferred a * a
             @inferred a * a'
+            @inferred kron(a)
             @inferred kron(a, σx)
             @inferred kron(a, eye(2))
         end
     end
 
+    @testset "tensor" begin
+        σx = sigmax()
+        X3 = kron(σx, σx, σx)
+        @test tensor(σx) == kron(σx)
+        @test tensor(fill(σx, 3)...) == X3
+        X_warn = @test_logs (
+            :warn,
+            "`tensor(A)` or `kron(A)` with `A` is a `Vector` can hurt performance. Try to use `tensor(A...)` or `kron(A...)` instead.",
+        ) tensor(fill(σx, 3))
+        @test X_warn == X3
+    end
+
     @testset "projection" begin
         N = 10
         ψ = fock(N, 3)

From c61937960c5e92f5aab5addf9a0a6dbd8a0adb1f Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Mon, 16 Sep 2024 00:43:15 +0800
Subject: [PATCH 07/10] add new section to docs

---
 docs/src/users_guide/tensor.md | 151 ++++++++++++++++++++++++++++++++-
 1 file changed, 150 insertions(+), 1 deletion(-)

diff --git a/docs/src/users_guide/tensor.md b/docs/src/users_guide/tensor.md
index faf4ae7d..4c7bdfbe 100644
--- a/docs/src/users_guide/tensor.md
+++ b/docs/src/users_guide/tensor.md
@@ -1,3 +1,152 @@
 # [Tensor Products and Partial Traces](@id doc:Tensor-products-and-Partial-Traces)
 
-This page is still under construction, please visit [API](@ref doc-API) first.
\ No newline at end of file
+```@setup tensor_products
+using QuantumToolbox
+```
+
+## [Tensor products](@id doc:Tensor-products)
+
+To describe the states of multipartite quantum systems (such as two coupled qubits, a qubit coupled to an oscillator, etc.) we need to expand the Hilbert space by taking the tensor product of the state vectors for each of the system components. Similarly, the operators acting on the state vectors in the combined Hilbert space (describing the coupled system) are formed by taking the tensor product of the individual operators.
+
+In `QuantumToolbox`, the function [`tensor`](@ref) (or [`kron`](@ref)) is used to accomplish this task. This function takes a collection of [`Ket`](@ref) or [`Operator`](@ref) as argument and returns a composite [`QuantumObject`](@ref) for the combined Hilbert space. The function accepts an arbitrary number of [`QuantumObject`](@ref) as argument. The `type` of returned [`QuantumObject`](@ref) is the same as that of the input(s).
+
+A collection of [`QuantumObject`](@ref):
+```@example tensor_products
+tensor(sigmax(), sigmax(), sigmax())
+```
+
+or a `Vector{QuantumObject}`:
+
+```@example tensor_products
+op_list = fill(sigmax(), 3)
+tensor(op_list)
+```
+
+!!! warning "Beware of type-stability!"
+    Please note that `tensor(op_list)` or `kron(op_list)` with `op_list` is a `Vector` is type-instable and can hurt performance. It is recommended to use `tensor(op_list...)` or `kron(op_list...)` instead. See the Section [The Importance of Type-Stability](@ref doc:Type-Stability) for more details.
+
+For example, the state vector describing two qubits in their ground states is formed by taking the tensor product of the two single-qubit ground state vectors:
+
+```@example tensor_products
+tensor(basis(2, 0), basis(2, 0))
+```
+
+One can generalize to more qubits by adding more component state vectors in the argument list to the [`tensor`](@ref) (or [`kron`](@ref)) function, as illustrated in the following example:
+
+```@example tensor_products
+states = QuantumObject[
+    normalize(basis(2, 0) + basis(2, 1)),
+    normalize(basis(2, 0) + basis(2, 1)),
+    basis(2, 0)
+]
+tensor(states...)
+```
+This state is slightly more complicated, describing two qubits in a superposition between the up and down states, while the third qubit is in its ground state.
+
+To construct operators that act on an extended Hilbert space of a combined system, we similarly pass a list of operators for each component system to the [`tensor`](@ref) (or [`kron`](@ref)) function. For example, to form the operator that represents the simultaneous action of the ``\sigma_x`` operator on two qubits:
+
+```@example tensor_products
+tensor(sigmax(), sigmax())
+```
+
+To create operators in a combined Hilbert space that only act on a single component, we take the tensor product of the operator acting on the subspace of interest, with the identity operators corresponding to the components that are to be unchanged. For example, the operator that represents ``\sigma_z`` on the first qubit in a two-qubit system, while leaving the second qubit unaffected:
+
+```@example tensor_products
+tensor(sigmaz(), qeye(2))
+```
+
+## Example: Constructing composite Hamiltonians
+
+The [`tensor`](@ref) (or [`kron`](@ref)) function is extensively used when constructing Hamiltonians for composite systems. Here we’ll look at some simple examples.
+
+### Two coupled qubits
+
+First, let’s consider a system of two coupled qubits. Assume that both the qubits have equal energy splitting, and that the qubits are coupled through a ``\sigma_x \otimes \sigma_x`` interaction with strength ``g = 0.05`` (in units where the bare qubit energy splitting is unity). The Hamiltonian describing this system is:
+
+```@example tensor_products
+H = tensor(sigmaz(), qeye(2)) + 
+    tensor(qeye(2), sigmaz()) + 
+    0.05 * tensor(sigmax(), sigmax())
+```
+
+### Three coupled qubits
+
+The two-qubit example is easily generalized to three coupled qubits:
+
+```@example tensor_products
+H = tensor(sigmaz(), qeye(2), qeye(2)) + 
+    tensor(qeye(2), sigmaz(), qeye(2)) + 
+    tensor(qeye(2), qeye(2), sigmaz()) + 
+    0.5  * tensor(sigmax(), sigmax(), qeye(2)) + 
+    0.25 * tensor(qeye(2), sigmax(), sigmax())
+```
+
+### A two-level system coupled to a cavity: The Jaynes-Cummings model
+
+The simplest possible quantum mechanical description for light-matter interaction is encapsulated in the Jaynes-Cummings model, which describes the coupling between a two-level atom and a single-mode electromagnetic field (a cavity mode). Denoting the energy splitting of the atom and cavity ``\omega_a`` and ``\omega_c``, respectively, and the atom-cavity interaction strength ``g``, the Jaynes-Cummings Hamiltonian can be constructed as:
+
+```math
+H = \frac{\omega_a}{2}\sigma_z + \omega_c a^\dagger a + g (a^\dagger \sigma_- + a \sigma_+)
+```
+
+```@example tensor_products
+N = 6     # cavity fock space truncation
+ωc = 1.25 # frequency of cavity
+ωa = 1.0  # frequency of two-level atom
+g = 0.75  # interaction strength
+
+a = tensor(qeye(2), destroy(N)) # cavity annihilation operator
+
+# two-level atom operators
+σm = tensor(destroy(2), qeye(N))
+σz = tensor(sigmaz(), qeye(N))
+
+H = 0.5 * ωa * σz + ωc * a' * a + g * (a' * σm + a * σm')
+```
+
+## [Partial trace](@id doc:Partial-trace)
+
+The partial trace is an operation that reduces the dimension of a Hilbert space by eliminating some degrees of freedom by averaging (tracing). In this sense it is therefore the converse of the tensor product. It is useful when one is interested in only a part of a coupled quantum system. For open quantum systems, this typically involves tracing over the environment leaving only the system of interest. In `QuantumToolbox` the function [`ptrace`](@ref) is used to take partial traces. [`ptrace`](@ref) takes one [`QuantumObject`](@ref) as an input, and also one argument `sel`, which marks the component systems that should be kept, and all the other components are traced out. 
+
+Remember that the index of `Julia` starts from `1`, and all the elements in `sel` should be positive `Integer`. Therefore, the type of `sel` can be either `Integer`, `Tuple`, `SVector`, or `Vector`.
+
+!!! warning "Beware of type-stability!"
+    Although it supports also `Vector` type, it is recommended to use `Tuple` or `SVector` from [`StaticArrays.jl`](https://github.com/JuliaArrays/StaticArrays.jl) to improve performance. For a brief explanation on the impact of the type of `sel`, see the section [The Importance of Type-Stability](@ref doc:Type-Stability).
+
+For example, the density matrix describing a single qubit obtained from a coupled two-qubit system is obtained via:
+
+```@example tensor_products
+ψ = tensor(
+    basis(2, 0), 
+    basis(2, 1), 
+    normalize(basis(2, 0) + basis(2, 1))
+)
+```
+
+```@example tensor_products
+ptrace(ψ, 1) # trace out 2nd and 3rd systems
+```
+
+```@example tensor_products
+ptrace(ψ, (1, 3)) # trace out 2nd system
+```
+
+Note that the partial trace always results in a [`Operator`](@ref) (density matrix), regardless of whether the composite system is a pure state (described by a [`Ket`](@ref)) or a mixed state (described by a [`Operator`](@ref)):
+
+```@example tensor_products
+ψ1 = normalize(basis(2, 0) + basis(2, 1))
+ψ2 = basis(2, 0)
+ψT = tensor(ψ1, ψ2)
+```
+
+```@example tensor_products
+ptrace(ψT, 1)
+```
+
+```@example tensor_products
+ρT = tensor(ket2dm(ψ1), ket2dm(ψ1))
+```
+
+```@example tensor_products
+ptrace(ρT, 1)
+```

From 61db86c29df2cf1982d3a569cf6fe679deb09751 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Mon, 16 Sep 2024 17:55:34 +0800
Subject: [PATCH 08/10] fix `ptrace` and add extended methods for `tensor` and
 `kron`

---
 src/qobj/arithmetic_and_attributes.jl | 71 ++++++++++++++++++---------
 src/qobj/functions.jl                 |  5 ++
 src/qobj/synonyms.jl                  |  2 +-
 test/quantum_objects.jl               | 57 +++++++++++++++++++++
 4 files changed, 110 insertions(+), 25 deletions(-)

diff --git a/src/qobj/arithmetic_and_attributes.jl b/src/qobj/arithmetic_and_attributes.jl
index 4ba7d4bc..66ac02c8 100644
--- a/src/qobj/arithmetic_and_attributes.jl
+++ b/src/qobj/arithmetic_and_attributes.jl
@@ -514,18 +514,38 @@ Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
 ```
 """
 function ptrace(QO::QuantumObject{<:AbstractArray,KetQuantumObject}, sel::Union{AbstractVector{Int},Tuple})
-    length(QO.dims) == 1 && return QO
+    ns = length(sel)
+    if ns == 0
+        return tr(QO * QO')
+    else
+        nd = length(QO.dims)
+        !all(nd .>= sel .> 0) && throw(
+            ArgumentError("Invalid indices in `sel`: $(sel), the given QuantumObject only have $(nd) sub-systems"),
+        )
+        ns != length(unique(sel)) && throw(ArgumentError("Duplicate selection indices in `sel`: $(sel)"))
+        nd == 1 && return QO * QO'  # ptrace should always return Operator
+    end
 
-    ρtr, dkeep = _ptrace_ket(QO.data, QO.dims, SVector(sel))
+    ρtr, dkeep = _ptrace_ket(QO.data, QO.dims, sort(SVector(sel)))
     return QuantumObject(ρtr, type = Operator, dims = dkeep)
 end
 
 ptrace(QO::QuantumObject{<:AbstractArray,BraQuantumObject}, sel::Union{AbstractVector{Int},Tuple}) = ptrace(QO', sel)
 
 function ptrace(QO::QuantumObject{<:AbstractArray,OperatorQuantumObject}, sel::Union{AbstractVector{Int},Tuple})
-    length(QO.dims) == 1 && return QO
+    ns = length(sel)
+    if ns == 0
+        return tr(QO)
+    else
+        nd = length(QO.dims)
+        !all(nd .>= sel .> 0) && throw(
+            ArgumentError("Invalid indices in `sel`: $(sel), the given QuantumObject only have $(nd) sub-systems"),
+        )
+        ns != length(unique(sel)) && throw(ArgumentError("Duplicate selection indices in `sel`: $(sel)"))
+        nd == 1 && return QO
+    end
 
-    ρtr, dkeep = _ptrace_oper(QO.data, QO.dims, SVector(sel))
+    ρtr, dkeep = _ptrace_oper(QO.data, QO.dims, sort(SVector(sel)))
     return QuantumObject(ρtr, type = Operator, dims = dkeep)
 end
 ptrace(QO::QuantumObject, sel::Int) = ptrace(QO, SVector(sel))
@@ -538,17 +558,20 @@ function _ptrace_ket(QO::AbstractArray, dims::Union{SVector,MVector}, sel)
     qtrace = filter(i -> i ∉ sel, 1:nd)
     dkeep = dims[sel]
     dtrace = dims[qtrace]
-    # Concatenate sel and qtrace without loosing the length information
-    sel_qtrace = ntuple(Val(nd)) do i
-        if i <= length(sel)
-            @inbounds sel[i]
+    nt = length(dtrace)
+
+    # Concatenate qtrace and sel without loosing the length information
+    # Tuple(qtrace..., sel...)
+    qtrace_sel = ntuple(Val(nd)) do i
+        if i <= nt
+            @inbounds qtrace[i]
         else
-            @inbounds qtrace[i-length(sel)]
+            @inbounds sel[i-nt]
         end
     end
 
     vmat = reshape(QO, reverse(dims)...)
-    topermute = nd + 1 .- sel_qtrace
+    topermute = reverse(nd + 1 .- qtrace_sel)
     vmat = permutedims(vmat, topermute) # TODO: use PermutedDimsArray when Julia v1.11.0 is released
     vmat = reshape(vmat, prod(dkeep), prod(dtrace))
 
@@ -563,27 +586,27 @@ function _ptrace_oper(QO::AbstractArray, dims::Union{SVector,MVector}, sel)
     qtrace = filter(i -> i ∉ sel, 1:nd)
     dkeep = dims[sel]
     dtrace = dims[qtrace]
-    # Concatenate sel and qtrace without loosing the length information
+    nk = length(dkeep)
+    nt = length(dtrace)
+    _2_nt = 2 * nt
+
+    # Concatenate qtrace and sel without loosing the length information
+    # Tuple(qtrace..., sel...)
     qtrace_sel = ntuple(Val(2 * nd)) do i
-        if i <= length(qtrace)
+        if i <= nt
             @inbounds qtrace[i]
-        elseif i <= 2 * length(qtrace)
-            @inbounds qtrace[i-length(qtrace)] + nd
-        elseif i <= 2 * length(qtrace) + length(sel)
-            @inbounds sel[i-length(qtrace)-length(sel)]
+        elseif i <= _2_nt
+            @inbounds qtrace[i-nt] + nd
+        elseif i <= _2_nt + nk
+            @inbounds sel[i-_2_nt]
         else
-            @inbounds sel[i-length(qtrace)-2*length(sel)] + nd
+            @inbounds sel[i-_2_nt-nk] + nd
         end
     end
 
     ρmat = reshape(QO, reverse(vcat(dims, dims))...)
-    topermute = 2 * nd + 1 .- qtrace_sel
-    ρmat = permutedims(ρmat, reverse(topermute)) # TODO: use PermutedDimsArray when Julia v1.11.0 is released
-
-    ## TODO: Check if it works always
-
-    # ρmat = row_major_reshape(ρmat, prod(dtrace), prod(dtrace), prod(dkeep), prod(dkeep))
-    # res = dropdims(mapslices(tr, ρmat, dims=(1,2)), dims=(1,2))
+    topermute = reverse(2 * nd + 1 .- qtrace_sel)
+    ρmat = permutedims(ρmat, topermute) # TODO: use PermutedDimsArray when Julia v1.11.0 is released
     ρmat = reshape(ρmat, prod(dkeep), prod(dkeep), prod(dtrace), prod(dtrace))
     res = map(tr, eachslice(ρmat, dims = (1, 2)))
 
diff --git a/src/qobj/functions.jl b/src/qobj/functions.jl
index 667f7930..19f41629 100644
--- a/src/qobj/functions.jl
+++ b/src/qobj/functions.jl
@@ -183,6 +183,11 @@ function LinearAlgebra.kron(
 ) where {T1,T2,OpType<:Union{KetQuantumObject,BraQuantumObject,OperatorQuantumObject}}
     return QuantumObject(kron(A.data, B.data), A.type, vcat(A.dims, B.dims))
 end
+LinearAlgebra.kron(A::QuantumObject) = A
+function LinearAlgebra.kron(A::Vector{<:QuantumObject})
+    @warn "`tensor(A)` or `kron(A)` with `A` is a `Vector` can hurt performance. Try to use `tensor(A...)` or `kron(A...)` instead."
+    return kron(A...)
+end
 
 @doc raw"""
     vec2mat(A::AbstractVector)
diff --git a/src/qobj/synonyms.jl b/src/qobj/synonyms.jl
index 5b991d5c..9457bbb6 100644
--- a/src/qobj/synonyms.jl
+++ b/src/qobj/synonyms.jl
@@ -184,7 +184,7 @@ Quantum Object:   type=Operator   dims=[2, 2, 2]   size=(8, 8)   ishermitian=tru
  1.0+0.0im      ⋅          ⋅             ⋅          ⋅          ⋅
 ```
 """
-tensor(A::QuantumObject...) = kron(A...)
+tensor(A...) = kron(A...)
 
 @doc raw"""
     ⊗(A::QuantumObject, B::QuantumObject)
diff --git a/test/quantum_objects.jl b/test/quantum_objects.jl
index e61ff996..2755257c 100644
--- a/test/quantum_objects.jl
+++ b/test/quantum_objects.jl
@@ -309,11 +309,24 @@
             @inferred a .^ 2
             @inferred a * a
             @inferred a * a'
+            @inferred kron(a)
             @inferred kron(a, σx)
             @inferred kron(a, eye(2))
         end
     end
 
+    @testset "tensor" begin
+        σx = sigmax()
+        X3 = kron(σx, σx, σx)
+        @test tensor(σx) == kron(σx)
+        @test tensor(fill(σx, 3)...) == X3
+        X_warn = @test_logs (
+            :warn,
+            "`tensor(A)` or `kron(A)` with `A` is a `Vector` can hurt performance. Try to use `tensor(A...)` or `kron(A...)` instead.",
+        ) tensor(fill(σx, 3))
+        @test X_warn == X3
+    end
+
     @testset "projection" begin
         N = 10
         ψ = fock(N, 3)
@@ -576,6 +589,50 @@
         @test ρ1.data ≈ ρ1_ptr.data atol = 1e-10
         @test ρ2.data ≈ ρ2_ptr.data atol = 1e-10
 
+        ψlist = [rand_ket(2), rand_ket(3), rand_ket(4), rand_ket(5)]
+        ρlist = [rand_dm(2), rand_dm(3), rand_dm(4), rand_dm(5)]
+        ψtotal = tensor(ψlist...)
+        ρtotal = tensor(ρlist...)
+        sel_tests = [
+            Int64[],
+            1,
+            2,
+            3,
+            4,
+            (1, 2),
+            (1, 3),
+            (1, 4),
+            (2, 3),
+            (2, 4),
+            (3, 4),
+            (1, 2, 3),
+            (1, 2, 4),
+            (1, 3, 4),
+            (2, 3, 4),
+            (1, 2, 3, 4),
+        ]
+        for sel in sel_tests
+            if length(sel) == 0
+                @test ptrace(ψtotal, sel) ≈ 1.0
+                @test ptrace(ρtotal, sel) ≈ 1.0
+            else
+                @test ptrace(ψtotal, sel) ≈ tensor([ket2dm(ψlist[i]) for i in sel]...)
+                @test ptrace(ρtotal, sel) ≈ tensor([ρlist[i] for i in sel]...)
+            end
+        end
+        @test ptrace(ψtotal, (1, 3, 4)) ≈ ptrace(ψtotal, (4, 3, 1))
+        @test ptrace(ρtotal, (1, 3, 4)) ≈ ptrace(ρtotal, (3, 1, 4))
+        @test_throws ArgumentError ptrace(ψtotal, 0)
+        @test_throws ArgumentError ptrace(ψtotal, 5)
+        @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(ρtotal, 0)
+        @test_throws ArgumentError ptrace(ρtotal, 5)
+        @test_throws ArgumentError ptrace(ρtotal, (0, 2))
+        @test_throws ArgumentError ptrace(ρtotal, (2, 5))
+        @test_throws ArgumentError ptrace(ρtotal, (2, 2, 3))
+
         @testset "Type Inference (ptrace)" begin
             @inferred ptrace(ρ, 1)
             @inferred ptrace(ρ, 2)

From ef5cbc3e4672e17a507ff7f9e92e75ac1ea63b98 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Tue, 17 Sep 2024 09:00:57 +0800
Subject: [PATCH 09/10] remove duplicated definition of `kron`

---
 src/qobj/functions.jl | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/src/qobj/functions.jl b/src/qobj/functions.jl
index 60cfd578..a69d8c10 100644
--- a/src/qobj/functions.jl
+++ b/src/qobj/functions.jl
@@ -197,11 +197,6 @@ function LinearAlgebra.kron(A::Vector{<:QuantumObject})
     @warn "`tensor(A)` or `kron(A)` with `A` is a `Vector` can hurt performance. Try to use `tensor(A...)` or `kron(A...)` instead."
     return kron(A...)
 end
-LinearAlgebra.kron(A::QuantumObject) = A
-function LinearAlgebra.kron(A::Vector{<:QuantumObject})
-    @warn "`tensor(A)` or `kron(A)` with `A` is a `Vector` can hurt performance. Try to use `tensor(A...)` or `kron(A...)` instead."
-    return kron(A...)
-end
 
 @doc raw"""
     vec2mat(A::AbstractVector)

From 8d49dc0b7759f19ad10de0a7c0824c53b19db568 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Tue, 17 Sep 2024 10:02:04 +0800
Subject: [PATCH 10/10] add non static array warning for `ptrace`

---
 src/qobj/arithmetic_and_attributes.jl | 37 ++++++++++++++++-----------
 test/quantum_objects.jl               | 12 +++++++--
 2 files changed, 32 insertions(+), 17 deletions(-)

diff --git a/src/qobj/arithmetic_and_attributes.jl b/src/qobj/arithmetic_and_attributes.jl
index 66ac02c8..71004632 100644
--- a/src/qobj/arithmetic_and_attributes.jl
+++ b/src/qobj/arithmetic_and_attributes.jl
@@ -475,8 +475,9 @@ proj(ψ::QuantumObject{<:AbstractArray{T},BraQuantumObject}) where {T} = ψ' * 
 @doc raw"""
     ptrace(QO::QuantumObject, sel)
 
-[Partial trace](https://en.wikipedia.org/wiki/Partial_trace) of a quantum state `QO` leaving only the dimensions
-with the indices present in the `sel` vector.
+[Partial trace](https://en.wikipedia.org/wiki/Partial_trace) of a quantum state `QO` leaving only the dimensions with the indices present in the `sel` vector.
+
+Note that this function will always return [`Operator`](@ref). No matter the input [`QuantumObject`](@ref) is a [`Ket`](@ref), [`Bra`](@ref), or [`Operator`](@ref)
 
 # Examples
 Two qubits in the state ``\ket{\psi} = \ket{e,g}``:
@@ -515,18 +516,21 @@ Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
 """
 function ptrace(QO::QuantumObject{<:AbstractArray,KetQuantumObject}, sel::Union{AbstractVector{Int},Tuple})
     ns = length(sel)
-    if ns == 0
-        return tr(QO * QO')
+    if ns == 0 # return full trace for empty sel
+        return tr(ket2dm(QO))
     else
         nd = length(QO.dims)
-        !all(nd .>= sel .> 0) && throw(
+
+        (any(>(nd), sel) || any(<(1), sel)) && throw(
             ArgumentError("Invalid indices in `sel`: $(sel), the given QuantumObject only have $(nd) sub-systems"),
         )
-        ns != length(unique(sel)) && throw(ArgumentError("Duplicate selection indices in `sel`: $(sel)"))
-        nd == 1 && return QO * QO'  # ptrace should always return Operator
+        (ns != length(unique(sel))) && throw(ArgumentError("Duplicate selection indices in `sel`: $(sel)"))
+        (nd == 1) && return ket2dm(QO) # ptrace should always return Operator
     end
 
-    ρtr, dkeep = _ptrace_ket(QO.data, QO.dims, sort(SVector(sel)))
+    _non_static_array_warning("sel", sel)
+    _sort_sel = sort(SVector{length(sel),Int}(sel))
+    ρtr, dkeep = _ptrace_ket(QO.data, QO.dims, _sort_sel)
     return QuantumObject(ρtr, type = Operator, dims = dkeep)
 end
 
@@ -534,18 +538,21 @@ ptrace(QO::QuantumObject{<:AbstractArray,BraQuantumObject}, sel::Union{AbstractV
 
 function ptrace(QO::QuantumObject{<:AbstractArray,OperatorQuantumObject}, sel::Union{AbstractVector{Int},Tuple})
     ns = length(sel)
-    if ns == 0
+    if ns == 0 # return full trace for empty sel
         return tr(QO)
     else
         nd = length(QO.dims)
-        !all(nd .>= sel .> 0) && throw(
+
+        (any(>(nd), sel) || any(<(1), sel)) && throw(
             ArgumentError("Invalid indices in `sel`: $(sel), the given QuantumObject only have $(nd) sub-systems"),
         )
-        ns != length(unique(sel)) && throw(ArgumentError("Duplicate selection indices in `sel`: $(sel)"))
-        nd == 1 && return QO
+        (ns != length(unique(sel))) && throw(ArgumentError("Duplicate selection indices in `sel`: $(sel)"))
+        (nd == 1) && return QO
     end
 
-    ρtr, dkeep = _ptrace_oper(QO.data, QO.dims, sort(SVector(sel)))
+    _non_static_array_warning("sel", sel)
+    _sort_sel = sort(SVector{length(sel),Int}(sel))
+    ρtr, dkeep = _ptrace_oper(QO.data, QO.dims, _sort_sel)
     return QuantumObject(ρtr, type = Operator, dims = dkeep)
 end
 ptrace(QO::QuantumObject, sel::Int) = ptrace(QO, SVector(sel))
@@ -560,7 +567,7 @@ function _ptrace_ket(QO::AbstractArray, dims::Union{SVector,MVector}, sel)
     dtrace = dims[qtrace]
     nt = length(dtrace)
 
-    # Concatenate qtrace and sel without loosing the length information
+    # Concatenate qtrace and sel without losing the length information
     # Tuple(qtrace..., sel...)
     qtrace_sel = ntuple(Val(nd)) do i
         if i <= nt
@@ -590,7 +597,7 @@ function _ptrace_oper(QO::AbstractArray, dims::Union{SVector,MVector}, sel)
     nt = length(dtrace)
     _2_nt = 2 * nt
 
-    # Concatenate qtrace and sel without loosing the length information
+    # Concatenate qtrace and sel without losing the length information
     # Tuple(qtrace..., sel...)
     qtrace_sel = ntuple(Val(2 * nd)) do i
         if i <= nt
diff --git a/test/quantum_objects.jl b/test/quantum_objects.jl
index f7e03734..365082d5 100644
--- a/test/quantum_objects.jl
+++ b/test/quantum_objects.jl
@@ -627,8 +627,16 @@
                 @test ptrace(ρtotal, sel) ≈ tensor([ρlist[i] for i in sel]...)
             end
         end
-        @test ptrace(ψtotal, (1, 3, 4)) ≈ ptrace(ψtotal, (4, 3, 1))
-        @test ptrace(ρtotal, (1, 3, 4)) ≈ ptrace(ρtotal, (3, 1, 4))
+        @test ptrace(ψtotal, (1, 3, 4)) ≈ ptrace(ψtotal, (4, 3, 1)) # check sort of sel
+        @test ptrace(ρtotal, (1, 3, 4)) ≈ ptrace(ρtotal, (3, 1, 4)) # check sort of sel
+        @test_logs (
+            :warn,
+            "The argument sel should be a Tuple or a StaticVector for better performance. Try to use `sel = (1, 2)` or `sel = SVector(1, 2)` instead of `sel = [1, 2]`.",
+        ) ptrace(ψtotal, [1, 2])
+        @test_logs (
+            :warn,
+            "The argument sel should be a Tuple or a StaticVector for better performance. Try to use `sel = (1, 2)` or `sel = SVector(1, 2)` instead of `sel = [1, 2]`.",
+        ) ptrace(ρtotal, [1, 2])
         @test_throws ArgumentError ptrace(ψtotal, 0)
         @test_throws ArgumentError ptrace(ψtotal, 5)
         @test_throws ArgumentError ptrace(ψtotal, (0, 2))