Numerical range of density matrix¶
Comment: the numerical range in this page is actually "joint algebraic numerical range" to be precise, but we will just call it numerical range for simplicity.
For more mathematical details of numerical range, please check doc-link.
Density matrix is postive semi-definite, Hermitian, trace one matrix.
$$ \{ \rho\in\mathbb{C}^{d\times d}: \rho\succeq 0,\mathrm{Tr}[\rho]=1 \} $$
a positive semi-definite matrix is always Hermitian, so we don't need to specify it.
import numpy as np
import matplotlib.pyplot as plt
try:
import numqi
except ImportError:
%pip install numqi
import numqi
One qubit case¶
For one qubit case, density matrix is a $2\times 2$ matrix, and it can be written as
$$ \rho=\rho_{0}+\frac{1}{2}\vec{n}\cdot\vec{\sigma} $$
where $\vec{n}\in \mathbb{R}^3$ is a three dimensional vector. To make sure the expression above PSD, the vector is required to be inside the Bloch sphere, i.e., $\lVert\vec{n}\rVert _2\leq 1$ wiki-link. And the numerical range of two operators is simply the projection of the Bloch sphere onto the plane spanned by the two operators. For this Bloch sphere, the projection is a circle if these two operators are perpendicular to each other.
np_rng = np.random.default_rng()
# generate two random, orthogonal unit vectors
tmp0 = np_rng.uniform(-1, 1, size=3)
vec0 = tmp0 / np.linalg.norm(tmp0)
op0 = vec0[0] * numqi.gate.X + vec0[1] * numqi.gate.Y + vec0[2] * numqi.gate.Z
tmp1 = np_rng.uniform(-1, 1, size=3)
tmp1 = tmp1 - tmp1.dot(vec0) * vec0 # orthogonalize
vec1 = tmp1 / np.linalg.norm(tmp1)
op1 = vec1[0] * numqi.gate.X + vec1[1] * numqi.gate.Y + vec1[2] * numqi.gate.Z
fig,ax = plt.subplots()
theta_list = np.linspace(0, 2*np.pi, 200)
direction = np.stack([np.cos(theta_list), np.sin(theta_list)], axis=1)
beta_list = numqi.matrix_space.get_joint_algebraic_numerical_range([op0,op1], direction)
ax.plot(beta_list*np.cos(theta_list), beta_list*np.sin(theta_list))
ax.set_aspect('equal')
ax.set_xlabel('$O_0$')
ax.set_ylabel('$O_1$')
_ = ax.set_title('two random orthogonal operators')
Two qubits case¶
For two qubits case, the shape of numerical range is much more complicated. Let's start with two commutative operators,
$$ O_0=X\otimes X,O_1=Z\otimes Z $$
one can verify that $[O_0,O_1]=0$.
op0 = np.kron(numqi.gate.X, numqi.gate.X)
op1 = np.kron(numqi.gate.Z, numqi.gate.Z)
# op1 = np.kron(numqi.gate.Z, numqi.gate.I) + np.kron(numqi.gate.I, numqi.gate.Z)
theta_list = np.linspace(0, 2*np.pi, 400)
direction = np.stack([np.cos(theta_list), np.sin(theta_list)], axis=1)
beta_dm = numqi.matrix_space.get_joint_algebraic_numerical_range([op0,op1], direction)
nr_dm = beta_dm.reshape(-1,1)*direction
fig,ax = plt.subplots(figsize=(8,6))
ax.plot(nr_dm[:,0], nr_dm[:,1], label='JANR')
ax.legend()
ax.set_xlabel('$O_0$')
ax.set_ylabel('$O_1$')
Text(0, 0.5, '$O_1$')
The numerical range of these two operators is a square with 4 vertices. We can figure out the vertices $(1,1)$ corresponds to one Bell state
$$ |\psi\rangle=\frac{1}{\sqrt{2}}(|00\rangle + |11\rangle) $$
for that
$$ \langle \psi | O_0 | \psi\rangle= \langle \psi | O_1 | \psi\rangle = 1$$
Task: can you identify the other three vertices?
Task: what about the numerical range of the operators $O_0=X\otimes X,O_1=Z\otimes I+I\otimes Z $
Then we may consider two random Hermitian operators $O_0,O_1$.
op0 = numqi.random.rand_hermitian_matrix(4)
op1 = numqi.random.rand_hermitian_matrix(4)
theta_list = np.linspace(0, 2*np.pi, 400)
direction = np.stack([np.cos(theta_list), np.sin(theta_list)], axis=1)
beta_dm = numqi.matrix_space.get_joint_algebraic_numerical_range([op0,op1], direction)
nr_dm = beta_dm.reshape(-1,1)*direction
fig,ax = plt.subplots(figsize=(8,6))
ax.plot(nr_dm[:,0], nr_dm[:,1], label='JANR')
ax.legend()
ax.set_xlabel('$O_0$')
ax.set_ylabel('$O_1$')
Text(0, 0.5, '$O_1$')
