Generalized Bloch vector¶
import itertools
import numpy as np
import matplotlib.pyplot as plt
try:
import numqi
except ImportError:
%pip install numqi
import numqi
np_rng = np.random.default_rng()
Density matrix is a Hermitian matrix with trace $1$. The density matrix can be expanded as a linear combination of the Gell-Mann matrices $\{M_i\}$ with real coefficients. For more details, please check doc-link.
For a $d$-dimensional system,
$$ \rho=\rho_0+\sum_{i=0}^{d^2-2}\vec{\rho}_iM_i $$
where $\vec{\rho}$ is a real vector with $d^2-1$ components, called generalized Bloch vector and $\rho_0$ is the maximally mixed state. Apparently, the maximally mixed state $\rho_0$ will correspond to the zero vector $\vec{\rho}=0$.
Comment: the definition of generalized Bloch vector may be slightly different from the orginal one.
rho0 = np.eye(3)/3
vec0 = numqi.gellmann.dm_to_gellmann_basis(rho0)
print(f'vec(maximally mixed state)=\n{vec0}')
rho = numqi.random.rand_density_matrix(3)
vec = numqi.gellmann.dm_to_gellmann_basis(rho)
print(f'vec(random density matrix)=\n{vec}')
vec(maximally mixed state)= [0. 0. 0. 0. 0. 0. 0. 0.] vec(random density matrix)= [-0.02981854 -0.00295902 0.0558201 0.16223101 0.18609394 -0.10494997 -0.12945388 -0.0673531 ]
Given a density matrix $\rho$, we can convert it into Gell-Mann basis. Its Bloch vector form $\vec{\rho}$, unit direction $\hat{\rho}$ and its vector norm $\lVert \vec{\rho}\rVert_2$.
$$ \rho=\rho_{0}+\vec{\rho}\cdot\vec{M}=\rho_{0}+\left\Vert \vec{\rho}\right\Vert _{2}\hat{\rho}\cdot\vec{M} $$
Two interpolation schemes might be useful
$$ \rho\left(\alpha\right)=\left(1-\alpha\right)\rho_{0}+\alpha\rho,\rho\left(\alpha=1\right)=\rho $$
$$ \rho\left(\beta,\hat{\rho}\right)=\rho_{0}+\beta\hat{\rho}\cdot\vec{M},\rho\left(\beta=\left\Vert \vec{\rho}\right\Vert _{2},\hat{\rho}\right)=\rho $$
Their connection is $\beta=\alpha \lVert \vec{\rho}\rVert_2$. In numqi, the second form is used in most cases to avoid confusion but you can always convert between them.
The purity $\gamma_{\rho}$ of a density matrix $\rho$ is closely related to the norm of the Bloch vector $\vec{\rho}$
$$\gamma_{\rho}=\mathrm{Tr}\left[ \rho ^2 \right] =\frac{1}{d}+2\lVert\vec{\rho}\rVert^{2}_2$$
It's easy to see that the purity is always larger than $1/d$ and the maximum mixed state $\rho=I/d$ gives the minimum purity $1/d$.
d = 4
rho = numqi.random.rand_density_matrix(d)
vec = numqi.gellmann.dm_to_gellmann_basis(rho)
purity = np.trace(rho @ rho).real
purity_from_vec = 1/d + 2*np.vdot(vec, vec).real
print(f'purity={purity}, purity_from_vec={purity_from_vec}')
purity=0.40977900973389725, purity_from_vec=0.40977900973389736
Easy to prove (see paper The Bloch Vector for N-Level Systems arxiv-link):
For pure states (rank 1):
$$ \mathrm{rank}(\rho)=1\Rightarrow \lVert \vec{\rho}\rVert_2^2=\frac{d-1}{2d} $$
For low-rank states (not full rank):
$$ \mathrm{rank}(\rho)<d\Rightarrow \max_{\rho}\lVert \vec{\rho}\rVert_2^2=\frac{d-1}{2d} $$
$$ \mathrm{rank}(\rho)<d\Rightarrow \min_{\rho}\lVert \vec{\rho}\rVert_2^2=\frac{1}{2d^2-2d} $$
where the maximum or minimum is taken over all density matrices with rank lower than $d$.
The set of density matrices is actually convex, where the low-rank states and full-rank states are on the boundary and interior seperately. Based on the above conculsion, we realize that in such a high-dimensional vector space, there exist a inscribed ball and a bounding ball with radius
$$ r_i^2=\frac{1}{2d^2-2d},\; r_o^2=\frac{d-1}{2d},$$
where the boundary of density matrices is between them.
However, imaging such a high-dimensional oject is quite difficult. Instead, we can study the cross sections of this oject. In previous tutorials, we have seen joint algebraic numerical range (JANR) of density matrix set $\{\rho\}$ for the given operators. However, numerical range is equivalent to a projection on the desired plane, which is different from cross sections. In the following, we will see how to plot the cross section of the density matrix set. Both the JANR and the cross section are useful for visualizing the density matrix set.
Example: 3-dimensional ball set¶
Before diving into the complicated density matrix set, let's begining with a 3-dimensional ball set
$$ B=\{(x,y,z): x^2+y^2+(z-0.8)^2\leq 1\}\subset \mathbb{R}^3 $$
which you may think as the Bloch Ball but the center is shifted to the position $(0,0,0.8)$. As for the joint algebraic numerical to the density matrix set, it's equivalent to make a projection along some plane for the ball $B$ here. As for the cross section, we can think about slicing the ball $B$ along some plane passing through the origin point $(0,0,0)$, and the intersection of the plane and the ball is the cross section of the ball $B$.
Below, we randomly pick a plane passing through the origin point $(0,0,0)$, and plot the cross section of the ball $B$. And the basis of the plane can be used to make the projection. And you can compare the difference between the cross section and the projection.
# generate two orthogonal vectors
np_rng1 = np.random.default_rng(233) #to make difference between projection and cross section significant
tmp0 = np_rng1.normal(size=3)
vec0 = tmp0 / np.linalg.norm(tmp0)
tmp0 = np_rng1.normal(size=3)
tmp0 = tmp0 - np.dot(tmp0, vec0) * vec0
vec1 = tmp0 / np.linalg.norm(tmp0)
# vec0 = np.array([1,0,0])
# vec1 = np.array([0,1,0])
radius = 1
center = np.array([0,0,0.8])
center_proj = np.array([np.dot(center,vec0), np.dot(center,vec1)])
radius_proj = radius
radius_cs = np.sqrt(radius**2 - (np.dot(center,center)-np.dot(center_proj,center_proj)))
fig,ax = plt.subplots()
theta_list = np.linspace(0,2*np.pi,101)
tmp0 = center_proj[0] + radius_proj * np.cos(theta_list)
tmp1 = center_proj[1] + radius_proj * np.sin(theta_list)
ax.plot(tmp0, tmp1, label='projection')
tmp0 = center_proj[0] + radius_cs * np.cos(theta_list)
tmp1 = center_proj[1] + radius_cs * np.sin(theta_list)
ax.plot(tmp0, tmp1, label='cross section')
ax.legend()
ax.set_xlabel('vec0')
ax.set_ylabel('vec1')
ax.set_aspect('equal')
fig.tight_layout()
print(f'vec0={vec0}, vec1={vec1}')
vec0=[0.59040693 0.54196798 0.59807221], vec1=[ 0.52966774 -0.8192958 0.21955975]
Usually, the difference between the projection and the cross section should be remarkable (if not, re-run it to randomize the two vectors $v_0$ and $v_1$). The projection is always of randius $1$ while the cross section will be smaller depends on the chosen plane.
task: at what plane will the difference between the projection and the cross section be the largest?
task: What if the center is not shifted and stays at the origin point $(0,0,0)$? can one still find the difference between the projection and the cross section?
For the one qubit case, the density matrix set, the Bloch ball, is a 3-dimensional ball with the center at the origin, so one cannot see the difference between the projection and the cross section. But for the qudit with $d>2$ (e.g. qutrit), the situation will be quite different.
Projection and cross section on random plane¶
Let us select two random operators, and see the difference between JANR and cross section.
dim = 4
tmp0 = np_rng.normal(size=dim*dim-1)
vec0 = tmp0 / np.linalg.norm(tmp0)
tmp0 = np_rng.normal(size=dim*dim-1)
tmp0 = tmp0 - np.dot(tmp0, vec0) * vec0
vec1 = tmp0 / np.linalg.norm(tmp0)
op0 = numqi.gellmann.gellmann_basis_to_dm(vec0) - np.eye(dim)/dim #make it traceless
op1 = numqi.gellmann.gellmann_basis_to_dm(vec1) - np.eye(dim)/dim
# op0 = np.diag([1,-1,0,0])
# op1 = np.diag([0,0,1,-1])
num_theta = 201
theta_list = np.linspace(0, 2*np.pi, num_theta)
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
angle_op0_op1,hf_plane = numqi.entangle.get_density_matrix_plane(op0, op1)
assert abs(angle_op0_op1-np.pi/2) < 1e-7 #as we have chosen orthogonal vectors (vec0 orth to vec1)
# factor 2 is for the normalization factor of the Gell-Mann matrices
beta_dm_cs = np.array([2*numqi.entangle.get_density_matrix_boundary(hf_plane(x))[1] for x in theta_list])
beta_dm_cs = beta_dm_cs.reshape(-1,1) * direction
fig,ax = plt.subplots(figsize=(8,6))
ax.plot(nr_dm[:,0], nr_dm[:,1], label='projection (JANR)')
ax.plot(beta_dm_cs[:,0], beta_dm_cs[:,1], label='cross section')
ax.legend()
ax.set_label('$O_0$')
ax.set_label('$O_1$')
Apparently, the cross section will always fall within the JANR, and both the cross section and the JANR are convex sets. When the number of the operators $\lvert \{A\}\rvert$ is equal to $d^2-1$, then the projection is simply a none operation, and the cross section will be the same as the JANR. However, we cannot visualize such high dimensional space, so we will only consider the case with two operators $\lvert \{A\}\rvert=2$.
task: can you find two operators that the cross section will touch the JANR at some points?
Cross sections on Gell-Mann planes¶
We can consider the cross sections of the planes spanned by different pairs of Gell-Mann matrices.
In the following figures, the orange dots are the inscribed ball and bounding ball, while the blue solid lines are the boundary of the density matrix set. Each figure is for one equivalent class of the cross section.
def get_ceil_square_root(N0:int):
# 2=1x2 3=2x2 4=2x2 5=2x3 6=2x3 7=3x3 9=3x3 10=3x4 12=3x4
assert N0>=1
tmp0 = int(np.ceil((np.sqrt(4*N0+1) - 1)/2))
tmp1 = int(np.ceil(np.sqrt(N0)))
ret = tmp0, (tmp0 if (tmp0==tmp1) else tmp0+1)
return ret
3-level system¶
dim = 3
gm_list = numqi.gellmann.all_gellmann_matrix(dim)
index01_set = list(itertools.combinations(list(range(dim*dim-1)), 2))
moment_list = [numqi.entangle.get_dm_cross_section_moment(gm_list[x], gm_list[y], order=1) for x,y in index01_set]
tmp0 = numqi.entangle.group_dm_cross_section_moment(moment_list, zero_eps=1e-4)
group_list = [[index01_set[y] for y in x] for x in tmp0]
index_to_cross_section = {x:numqi.entangle.get_dm_cross_section_boundary(gm_list[x[0]], gm_list[x[1]], num_point=201)['beta_dm'] for x in index01_set}
for ind0,group in enumerate(group_list):
print(f'[group-{ind0}]', group)
# (dim, len(group)): (3,4), (4,10), (5,17)
[group-0] [(0, 1), (0, 2), (0, 3), (0, 4), (0, 5), (0, 6), (1, 2), (1, 3), (1, 4), (1, 5), (2, 3), (2, 4), (2, 5), (3, 4), (3, 5), (3, 6), (4, 5)] [group-1] [(0, 7), (3, 7), (6, 7)] [group-2] [(1, 6), (2, 6), (4, 6), (5, 6)] [group-3] [(1, 7), (2, 7), (4, 7), (5, 7)]
fig,tmp0 = plt.subplots(*get_ceil_square_root(len(group_list)), figsize=(9,7))
ax_list = [y for x in tmp0 for y in x]
for ax,group in zip(ax_list, group_list):
for ind0,ind1 in group:
numqi.entangle.plot_dm_cross_section(index_to_cross_section[(ind0,ind1)], dim=dim, ax=ax, tag_show_legend=False)
ax.set_aspect('equal')
for x in ax_list[len(group_list):]:
x.axis('off')
fig.suptitle(f'Cross Section of {dim}-Level Systems')
fig.tight_layout()
In 3rd figure, there are two solid lines and they are equal up to a reflection. Qutrit ($d=3$) has 4 inequivalent cross sections. Similarly, we can perform the calculation for qudit $d=4$ which has 17 equivalent class as follows.
4-level system¶
dim = 4
gm_list = numqi.gellmann.all_gellmann_matrix(dim)
index01_set = list(itertools.combinations(list(range(dim*dim-1)), 2))
moment_list = [numqi.entangle.get_dm_cross_section_moment(gm_list[x], gm_list[y], order=1) for x,y in index01_set]
tmp0 = numqi.entangle.group_dm_cross_section_moment(moment_list, zero_eps=1e-4)
group_list = [[index01_set[y] for y in x] for x in tmp0]
index_to_cross_section = {x:numqi.entangle.get_dm_cross_section_boundary(gm_list[x[0]], gm_list[x[1]], num_point=201)['beta_dm'] for x in index01_set}
fig,tmp0 = plt.subplots(*get_ceil_square_root(len(group_list)), figsize=(9,7))
ax_list = [y for x in tmp0 for y in x]
for ax,group in zip(ax_list, group_list):
for ind0,ind1 in group:
numqi.entangle.plot_dm_cross_section(index_to_cross_section[(ind0,ind1)], dim=dim, ax=ax, tag_show_legend=False)
ax.set_aspect('equal')
for x in ax_list[len(group_list):]:
x.axis('off')
fig.suptitle(f'Cross Section of {dim}-Level Systems')
fig.tight_layout()
