8  Representing Quantum States and Operators with NumPy

In quantum mechanics, states and observables are represented using the algebra of Hilbert spaces. However, their infinite dimensions is incompatible with numerical simulations, that always requires finite elements. Hence, we truncate Hilbert spaces to a finite size, allowing us to run the quantum calculation on a computer. We can thus say that the whole problem of numerical quantum mechanics is then reduced to a problem of linear algebra. However, the intricated tensor structures of many-body Hilbert spaces requires also a powerful organization of the code and an easy way to access relevant information.

In the following we consider a system with a Hilbert space of dimension \(d\). The set of basis states \(\{|k\rangle: k=1, \ldots, d\}\) form an orthonormal basis, i.e., \(\left\langle k \mid k^{\prime}\right\rangle=\delta_{k, k^{\prime}}\). In general there are systems with an infinite dimensional Hilbert space, or systems, where the dimension is too large to be tractable on a computer. In this case \(d\) denotes the number of truncated basis states, which is used in the numerical simulation. For a given choice of basis states we can express any state vector and any operator as

\[ |\psi\rangle = \sum_{k=1}^{d} c_{k} |k\rangle, \quad \hat{A}=\sum_{k, l} A_{k l}|k\rangle\langle l|, \]

where \(c_{k} = \langle k \mid \psi\rangle\) and \(A_{k l}=\langle k| \hat{A}|l\rangle\). Therefore, in numerical simulations we represent states by vectors and operators by matrices according to the mapping

\[ |\psi\rangle \mapsto \vec{\psi}=\left(\begin{array}{c} c_{1} \\ c_{2} \\ \vdots \\ c_{d} \end{array}\right), \quad \hat{A} \mapsto A=\left(\begin{array}{cccc} A_{11} & A_{12} & \ldots & A_{1 d} \\ A_{21} & A_{22} & \ldots & A_{2 d} \\ \vdots & \vdots & \ddots & \vdots \\ A_{d 1} & A_{d 2} & \ldots & A_{d d} \end{array}\right) . \]

The left and right operations of an operator on a vector then simply translate into matrix vector multiplications,

\[ \hat{A}|\psi\rangle \mapsto \texttt{np.dot(A, psi)}, \quad \langle \psi| \hat{A} \mapsto \texttt{np.dot(np.conj(psi.T), A)}, \]

where in Numpy np.conj(psi.T) is the hermitian transpose of a matrix or vector.

8.1 Pauli Operators

The Pauli operators are fundamental in quantum mechanics, especially in the context of qubits. They are represented as matrices in a two-dimensional Hilbert space, which is the simplest non-trivial quantum system.

\[ \hat{\sigma}_{x} \mapsto\left(\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right), \quad \hat{\sigma}_{y} \mapsto\left(\begin{array}{cc} 0 & i \\ -i & 0 \end{array}\right), \quad \hat{\sigma}_{z} \mapsto\left(\begin{array}{cc} 1 & 0 \\ 0 & -1 \end{array}\right) . \]

In Numpy we simply define the corresponding matrices

import numpy as np

sx = np.array([[0, 1], [1, 0]])
sy = np.array([[0, 1j], [-1j, 0]])
sz = np.array([[1, 0], [0, -1]])

8.2 Spin S Systems

A system with a fixed total spin \(S\) is described by the three collective spin operators \(S_{x}, S_{y}, S_{z}\), which obey \(\left[S_{x}, S_{y}\right]=i S_{z}\). The \(d=2 S+1\) basis states can be labeled as \(|s, m\rangle\), where \(m=S, S-1, \ldots,-S\). By introducing spin raising and lowering operators \(S_{ \pm}=S_{x} \pm i S_{y}\), where \(S_{-}^{\dagger}=S_{+}\), all matrix elements can be obtained from the two relations

\[ S_{z}|s, m\rangle=m|s, m\rangle, \quad S_{-}|s, m\rangle=\sqrt{S(S+1)-m(m-1)}|s, m-1\rangle=: s_{-}^{m}|s, m-1\rangle \, . \]

Explicitly,

\[ \hat{S}_{z} \mapsto S Z=\left(\begin{array}{cccc} S & 0 & \ldots & 0 \\ 0 & S-1 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & -S \end{array}\right), \quad \hat{S}_{-} \mapsto S M=\left(\begin{array}{cccc} 0 & 0 & \ldots & 0 \\ s_{-}^{S / 2} & 0 & \ldots & 0 \\ 0 & s_{-}^{S / 2-1} & \ddots & \vdots \\ \vdots & 0 & \ldots & 0 \end{array}\right) . \]

All other operators can be obtained using \(S_{x}=\left(S_{+}+S_{-}\right) / 2\) and \(S_{y}=i\left(S_{+}-S_{-}\right) / 2\). Example for a spin \(S=1\) system

s = 1
d = int(2 * s + 1)

# Vector of the diagonal elements for the SZ operator
vecm = np.flip(np.arange(-s, s + 1))

SZ = np.diag(vecm, 0)

# Vector for SM operator (ladder down)
vec2 = vecm[:d-1]  # remove the last entry
vec3 = np.sqrt(s * (s + 1) - vec2 * (vec2 - 1))
SM = np.diag(vec3, k=-1)  # place vec3 in lower diagonal

# Construct SX and SY
SX = (SM + SM.T) / 2
SY = 1j * (SM.T - SM) / 2

Note that in Numpy the command A * B is a element-wise multiplication of two matrices, while A @ B would implement the usual matrix multiplication.

8.3 Harmonic Oscillator

In a Hilbert space of dimension \(N\), quantum states can be represented as vectors, and operators as matrices. Here we demonstrate the destroy operator, \(a\), which lowers the state by one quantum number. For a detailed discussion on the quantum harmonic oscillator and the bosonic annihilation operator, refer to Appendix A.

For a harmonic oscillator with number states \(|n\rangle\) the only nonzero matrix elements of the annihilation operator \(\hat{a}\) are given by \(\langle n-1| \hat{a}|n\rangle=\sqrt{n}\)

\[ \hat{a} \mapsto A=\left(\begin{array}{cccccc} 0 & 1 & \ldots & \cdots & \ldots & 0 \\ 0 & 0 & \sqrt{2} & \ldots & \cdots & 0 \\ 0 & 0 & 0 & \sqrt{3} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \ddots & \vdots \\ \vdots & \vdots & \vdots & \vdots & 0 & \sqrt{d-1} \\ 0 & 0 & 0 & \cdots & \cdots & 0 \end{array}\right) \]

This operator acts on Fock states to lower their quantum number by one, with a factor of \(\sqrt{n}\), where \(n\) is the quantum number of the initial state. In other words, \(\hat{a}|n\rangle = \sqrt{n}|n-1\rangle\). In the following code, we define the destroy operator by using NumPy, and we also define some Fock states for demonstration.

In Numpy we use the command np.diag(v, k=r), which creates a diagonal matrix with the elements of the vector v placed in the \(r\)-th diagonal \((r=0, \pm 1, \pm 2, \ldots)\).

def destroy(d):
    # creates a vector of the d-1 off-diagonal elements
    v=np.sqrt( np.arange(d-1) )
    # matrix with the elements of vec placed in the upper diagonal
    a=np.diag(v,k=1)
    return a

# Define the fock states
def fock(d, i):
    res = np.zeros(d)
    res[i] = 1
    return res

d = 7
zero_state = fock(d, 0)
one_state = fock(d, 1)
two_state = fock(d, 2)
three_state = fock(d, 3)

destroy_operator = destroy(d)
destroy_operator

array([[0.        , 0.        , 0.        , 0.        , 0.        ,
        0.        , 0.        ],
       [0.        , 0.        , 1.        , 0.        , 0.        ,
        0.        , 0.        ],
       [0.        , 0.        , 0.        , 1.41421356, 0.        ,
        0.        , 0.        ],
       [0.        , 0.        , 0.        , 0.        , 1.73205081,
        0.        , 0.        ],
       [0.        , 0.        , 0.        , 0.        , 0.        ,
        2.        , 0.        ],
       [0.        , 0.        , 0.        , 0.        , 0.        ,
        0.        , 2.23606798],
       [0.        , 0.        , 0.        , 0.        , 0.        ,
        0.        , 0.        ]])

Other operators (e.g., \(\hat{a}^{\dagger}, \hat{a}^{\dagger} \hat{a}\)) can be obtained by a hermitian transpose

\[ \hat{a}^{\dagger} \mapsto\texttt{np.conj(a.T)} \]

and matrix multiplications

\[ \hat{a}^{\dagger} \hat{a} \mapsto \texttt{np.matmul(np.conj(a.T) , a)} \, . \]

Note that in some cases this introduces truncation artifacts. For example, the matrix for the operator \(M=\texttt{np.matmul(a, np.conj(a.T))}\) has a zero diagonal element \(M[d, d]=0\) inherited from the matrix \(\texttt{np.conj(a.T))}\), while the same operator constructed in a different way, \(M_2=\texttt{np.conj(a.T) * a + np.eye(d)}\), does not. This can be avoided by constructing this operator explicitly. Note that this type of truncation artifacts are related to the fact that in a infinite Hilbert space \(\mathrm{Tr}([a, a^{\dagger}])\neq 0\) (actually, striclty speaking, \(=\infty\)) as a consequence of the canonical commutation relation. On the contrary, in a finite Hilbert space for any two operators \(O_1, O_2\), \(\mathrm{Tr}([O_1, O_2])= 0\). Taking a dimension \(d\) large enough allows to make these artifacts a negligible error in the whole computation.

8.3.1 Action of the Destroy Operator on a Fock State

The action of the destroy operator \(a\) on a Fock state \(|n\rangle\) lowers the state by one quantum number, multiplied by a factor \(\sqrt{n}\). For example, applying \(a\) to the state \(|3\rangle\) yields:

\[ \hat{a}|3\rangle = \sqrt{3}|2\rangle \]

This demonstrates the lowering action of the destroy operator with a specific factor, dependent on the quantum number of the state being acted upon.

# Apply the destroy operator on the one state
result_state = np.dot(destroy_operator, three_state)

print("Resulting State:")
result_state

Resulting State:

array([0.        , 0.        , 1.41421356, 0.        , 0.        ,
       0.        , 0.        ])

8.4 Partial Trace

In Section 4.1.4, we have already discussed the concept of tensor products. Here we will introduce the partial trace, a crucial operation in quantum mechanics that allows us to focus on a subsystem of a larger composite system.

The partial trace over a subsystem, say \(B\), of a composite system \(AB\), mathematically expresses as “tracing out” \(B\), leaving the reduced state of \(A\). For a bipartite state \(\rho_{AB}\), the partial trace over \(B\) is:

\[ \text{Tr}_B(\hat{\rho}_{AB}) = \sum_{i \in \mathcal{H}_B} \langle i| \hat{\rho}_{AB} |i\rangle \]

where \(\{|i\rangle\}\) forms a complete basis for subsystem \(B\).

Let’s try it with an entangled Bell’s state between two qubits:

\[ \vert \phi^+ \rangle = \frac{1}{\sqrt{2}} \left( \vert 0, 0 \rangle + \vert 1, 1 \rangle \right) \]

def ptrace(psi, subspace_to_keep, dim_subspace):
    dim1, dim2 = dim_subspace

    rho = np.outer(psi, psi.conj())

    # Reshape rho to separate the subsystems' degrees of freedom
    rho_reshaped = rho.reshape(dim1, dim2, dim1, dim2)

    if subspace_to_keep == 1:
        # Perform the trace over the second subsystem
        traced_out = np.trace(rho_reshaped, axis1=1, axis2=3)
    elif subspace_to_keep == 2:
        # Perform the trace over the first subsystem
        traced_out = np.trace(rho_reshaped, axis1=0, axis2=2)
    else:
        raise ValueError("subspace_to_keep must be either 1 or 2.")

    return traced_out

# Bell state between two qubits
phi_plus = ( np.kron(fock(2, 1), fock(2, 1)) + np.kron(fock(2, 0), fock(2, 0)) ) / np.sqrt(2)

# Reduced density matrix of the first qubit
rho_1 = ptrace(phi_plus, 1, (2, 2))
rho_1

array([[0.5, 0. ],
       [0. , 0.5]])

8.5 Why QuTiP?

While NumPy and SciPy are powerful tools for numerical computations, they lack specific functionalities for efficiently handling complex quantum systems. QuTiP is designed to fill this gap, offering features such as:

• Easy manipulation and visualization of quantum objects.
• Support for operations on states and operators in different Hilbert spaces.
• Tools for dealing with composite systems, partial traces, and superoperators. It is like to have the book “Quantum noise” (by Gardiner and Zoller) already implemented in your laptop!

In the next chapters, we’ll explore how QuTiP simplifies these tasks, making it an invaluable tool for quantum optics simulations.