So far we've essentially discussed "bosonic" second quantization, where we assign a harmonic oscillator to each degree of freedom of some single particle Hilbert space, and the number operator of each oscillator keeps track of how many particles are in that state. It works out that this "Fock space" corresponds dimensionally to a tower of permutation symmetric tensor products of $0, 1, 2, 3, \dots$ particles. In other words, the use of harmonic oscillators automatically imposes for us that condition that the second quantized particles be permutation symmetric.
It's worth remember that although the creation and annihilation operators are not hermitian, they can be used as building blocks for hermitian operators. And if we have multiple oscillators, we can just tensor these creation and annihilation operators with identities so they act on the appropriate subspaces.
For bosons, if we have multiple creation and annihilation operators, they should obey the following commutation relations:
$$ [b_{i}, b_{j}] = [b^{\dagger}_{i}, b^{\dagger}_{j}] = 0 $$$$ [b_{i}, b^{\dagger}_{j}] = \delta^{i}_{j} $$Where $ b^{\dagger} $ is a creation operator, $ b $ is an annihilation operator, and: the commutator is $ [A, B] = AB - BA $.
In practice, if we work in finite dimensions, we can cap the maximum number of quanta in a given harmonic oscillator, and get finite dimensional "creation" and "annihilation" operators, which add and subtract quanta from each oscillator, although they don't satisfy perfectly the correct commutation relations.
Okay, so here's the thing:
If we want to second quantize fermionically, as opposed to bosonically, we have to make some changes. First, our oscillators can only have 0 or 1 excitations. In other words, at maximum only 1 fermion can be in a given state. This is the Pauli exclusion principle at work. The whole thing has to do with the fact while fermions and bosons are both "indistinguishable" (with respect to the tensor product) particles, fermions must always be in a permutation antisymmetric state, unlike bosons which must be in a permutation symmetric state. Look what happens if we try to antisymmetrize two identical states via $A, B \rightarrow A \otimes B - B \otimes A$: $\mid \uparrow \rangle \mid \uparrow \rangle - \mid \uparrow \rangle \mid \uparrow \rangle = 0$: we get the zero vector. The Pauli exclusion principle is just a fancy way of saying that fermions, by nature, are permutation antisymmetric particles. If we permute them, then the whole tensor state is multiplied by -1. This just changes the overall phase, and not the probabilities: indeed, the particles are still indistinguishable, despite them all picking up a sign flip when they are swapped in the tensor product.
For this reason, in terms of second quantization, it doesn't matter the order in which we create bosons, but it does matter for fermions: $ f^{\dagger}_{i} f^{\dagger}_{j} \mid 0 \rangle = -f^{\dagger}_{j} f^{\dagger}_{i} \mid 0 \rangle $, where $ \mid 0 \rangle $ is the fermion vacuum aka the state of all oscillators having 0 quanta.
The upshot is that the commutation relations between fermions involve the anticommutator $ \{A, B\} = AB + BA $, instead of the commutator.
$$ \{f_{i}, f_{j}\} = \{f^{\dagger}_{i}, f^{\dagger}_{j}\} = 0 $$$$ \{f_{i}, f^{\dagger}_{j}\} = \delta^{i}_{j} $$Luckily, with a trick, it's not hard to construct the right matrix representation for these operators. Indeed, it's kinda nice that the fermion excitations are capped at 1, since we can use $2$ x $2$ matrices for their creation and annihilation operators, and we can deal with a nice actually finite dimensional vector space.
For example, suppose we have 5 fermions. The standard $2$ x $2$ matrices for creation and annihilation operators are just:
$$ f^{\dagger} = \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix} $$$$ f = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} $$But here's the catch: to get 5 annihilation operators with the correct commutation relations, preserving the antisymmetry of the fermions, we take:
$$ f_{0} = f \otimes I \otimes I \otimes I \otimes I $$$$ f_{1} = Z \otimes f \otimes I \otimes I \otimes I $$$$ f_{2} = Z \otimes Z \otimes f \otimes I \otimes I $$$$ f_{3} = Z \otimes Z \otimes Z \otimes f \otimes I $$$$ f_{4} = Z \otimes Z \otimes Z \otimes Z \otimes f $$The idea is that there is a "normal ordering" for the fermions. When the operators are applied to the vacuum in the descending order, the $Z$'s in $ f_{4} $, don't matter since there are no excitations in those oscillators: we're in the vacuum after all. But when the operators are applied in the reverse order, they pick up a negative sign from the $Z$'s. You can check that this works:
import qutip as qt
import numpy as np
def anticommutator(a, b):
return a*b + b*a
#########################################################################################
def fermion_operators(n):
return [qt.tensor(*[qt.destroy(2) if i == j\
else (qt.sigmaz() if j < i\
else qt.identity(2))\
for j in range(n)])\
for i in range(n)]
def test_fermion_operators(f):
for i in range(len(f)):
for j in range(len(f)):
d = f[i].shape[0]
test1 = anticommutator(f[i], f[j]).full()
test2 = anticommutator(f[i], f[j].dag()).full()
if not \
(np.isclose(test1, np.zeros((d,d))).all()\
and \
((np.isclose(test2, np.zeros((d,d))).all() and i != j)\
or (np.isclose(test2, np.eye(d)).all() and i == j))):
return False
return True
#########################################################################################
n = 3
IDn = qt.identity(2**n)
IDn.dims = [[2]*n, [2]*n]
f = fermion_operators(n)
print(test_fermion_operators(f))
N = sum([a.dag()*a for a in f]) # number operator
I = qt.basis(2**n, 0) # vacuum state
I.dims = [[2]*n, [1]*n]
Okay, so we have $n$ fermionic oscillators to play with. And we want to second quantize an $n$ dimensional single particle Hilbert space fermionically. So now we have to talk about what the tower of antisymmetric multi-particle states looks like.
First, we observe that $sgn(\sigma) = (-1)^{T(\sigma)}$, where $\sigma$ is some permutation: e.g., $132$ would be a permutation of $123$. $T(\sigma)$ is the number of transpositions in the permutation. And $sgn(\sigma)$ is known as the parity of the permutation.
If we want to antisymmetrize a set of states, we basically want to sum over all permutations of the states, forming the tensor product of the pieces in each order, multiplying each tensor product by $-1$ for each swap in the permutation:
$ \sum_{\sigma} (-1)^{T(\sigma)} ( \sigma_{i} \otimes \sigma_{j } \otimes \dots)$
Using this, we can easily calculate the basis states for the permutation antisymmetric subspace of $n$ particles who each live in a $d$ dimensional Hilbert space.
We take each set of the $d$ basis vectors, each set that has $n$ elements and contains no repeats, and antisymmeterize those basis vectors. We'll actually end up with $\begin{pmatrix} d \\ n \end{pmatrix}$ basis vectors. We can form a rectangular matrix out of them, which will let us translate between a state of $\begin{pmatrix} d \\ n \end{pmatrix}$ dimensions, and a state of $n$ antisymmeterized particle states.
import qutip as qt
from qutip.qip.operations import swap
import itertools
def perm_parity(lst):
parity = 1
for i in range(0,len(lst)-1):
if lst[i] != i:
parity *= -1
mn = min(range(i,len(lst)), key=lst.__getitem__)
lst[i],lst[mn] = lst[mn],lst[i]
return parity
def antisymmetrize(*states):
a = sum([perm_parity(list(perm))*qt.tensor(*[states[p] for p in perm])\
for perm in itertools.permutations(list(range(len(states))))])
if a.norm() != 0:
a = a.unit()
return a
def antisymmetric_basis(n, d=2):
states = []
for perm in itertools.combinations_with_replacement(list(range(d)), r=n):
a = antisymmetrize(*[qt.basis(d, p) for p in perm])
if a.norm() != 0:
states.append(a)
Q = qt.Qobj(np.array([state.full().T[0].tolist() for state in states]))
Q.dims[1] = [d]*n
return Q
a = qt.rand_ket(2)
print(a)
b = qt.rand_ket(2)
print(b)
c = antisymmetrize(a,b)
print(c)
print(swap()*c)
antisymmetric_basis(3, d=3)
Now we put the two together.
We tensor, say, $3$ fermionic oscillators, and construct the total number operator. Then we construct a permutation matrix that rearranges the number operator so that its eigenvalues are in order. We can also use this permutation matrix $P$ to rearrange a Fock state so that it's of the form: $H_{1} \oplus H_{3} \oplus H_{3} \oplus H_{1}$. The dimensionality of these Hilbert spaces corresponds to the dimensionality of the antisymmetric subspace of the tensor product of $0, 1, 2$, and $3$ particles. So once we've picked apart out state, we can apply the linear map constructed above which sends:
$H_{1} \oplus H_{3} \oplus H_{3} \oplus H_{1} \rightarrow 0 \oplus A(H_{3}) \oplus A(H_{3}, H_{3}) \oplus A(H_{3}, H_{3}, H_{3})$, where $A()$ denotes the antisymmetric subspace of the tensor product of the particles.
import qutip as qt
import numpy as np
from itertools import product
def construct_permutation(n):
tensor_basis_labels = list(product([0,1], repeat=n))
total_n_basis_labels = []
for i in range(n+1):
total_n_basis_labels.extend(\
list(filter(lambda x: sum(x) == i, list(product([0,1], repeat=n))))[::-1])
P = np.zeros((2**n, 2**n))
for i, label in enumerate(tensor_basis_labels):
P[total_n_basis_labels.index(label)][i] = 1
P = qt.Qobj(P)
P.dims = [[2]*n, [2]*n]
sums = [sum(label) for label in total_n_basis_labels]
unique_sums = set(sums)
dims = [sums.count(us) for us in unique_sums]
return P, dims
def extract_states(q, dims):
v = q.full().T[0]
running = 0
blocks = []
for d in dims:
blocks.append(qt.Qobj(v[running:running+d]))
running += d
return blocks
def osc_antisymmetrics(state):
n = len(state.dims[0])
P, dim = construct_permutation(n)
extracted = extract_states(P*state, dim)
finished = [extracted[0]]
for i in range(1, len(extracted)):
finished.append(antisymmetric_basis(i, d=n).dag()*extracted[i])
return finished
#########################################################################################
n = 3
IDn = qt.identity(2**n)
IDn.dims = [[2]*n, [2]*n]
f = fermion_operators(n)
N = sum([a.dag()*a for a in f])
print(N)
P, dim = construct_permutation(n)
N_ = P*N*P.dag()
print(N_)
state = qt.rand_ket(2**n)
state.dims = [[2]*n, [1]*n]
print("\n Extracted sectors:")
print(extract_states(P*state, dim))
print("\nTower of antisymmetric states:\n")
antisymmetric_tower = osc_antisymmetrics(state)
print(antisymmetric_tower)
So there you have it: fermionic quantization.
It is often said that bosons are permutation symmetric particles with integer spin and fermions are permutation antisymmetric particles with half integer spin. This doesn't actually follow from anything we've discussed so far, even though we know quite a lot about spin, and also how to quantize things bosonically and fermionically. What it has to do with is the so-called "spin-statistics theorem."
First, a word about integer vs. half integer spin. This is the first clue. After all, what is the effect of representing rotations with $SU(2)$ acting on spinors instead of $SO(3)$ acting on $(x, y, z)$ vectors? We actually get something extra. $SU(2)$ isn't just identical to the 3D rotation group: in fact, it's the double cover of the 3D rotation group: for each element of the 3D rotation group, there are two elements of $SU(2)$. Take a spin-$\frac{1}{2}$ state and rotate it a full turn by exponentiating some Pauli matrix. You'll find that the state actually comes back to itself times $-1$. It'll only come back to itself exactly if you give it another full turn. So you actually have to turn a spin-$\frac{1}{2}$ state $720^{\circ}$ to make a full "rotation": it takes twice has much. Hence, the double cover: there are effectively two copies of the rotation group sitting inside $SU(2)$ so they can keep track of this negative sign. This is true for all half integer spin states: a full turn takes them to the negative of themselves, but the phase will wind around a bunch of times as it gets there. In contrast, integer spin states come back to themselves under a $360^{\circ}$ rotation (although the phase might wind around a lot for higher states.
import qutip as qt
import numpy as np
import vpython as vp
from examples.magic import *
def get_phase(q):
c = sum(q.full().T[0][::-1])
return np.exp(1j*np.angle(c))
scene = vp.canvas(background=vp.color.white)
j = 3/2
n = int(2*j + 1)
dt = 0.001
state = qt.rand_ket(n)
H = qt.jmat(j, 'y')
U = (-1j*H*dt).expm()
vp.sphere(color=vp.color.blue, opacity=0.5)
initial_vstars = [vp.sphere(pos=vp.vector(*xyz),\
radius=0.2, emissive=True,\
color=vp.color.yellow)
for xyz in spin_XYZ(state)]
vstars = [vp.sphere(pos=vp.vector(*xyz),\
radius=0.2, emissive=True)
for xyz in spin_XYZ(state)]
phase = get_phase(state)
initial_vphase = vp.arrow(pos=vp.vector(0,2,0),color=vp.color.yellow,\
axis=vp.vector(phase.real, phase.imag, 0))
vphase = vp.arrow(pos=vp.vector(0,2,0),\
axis=vp.vector(phase.real, phase.imag, 0))
while True:
state = U*state
for i, xyz in enumerate(spin_XYZ(state)):
vstars[i].pos = vp.vector(*xyz)
phase = get_phase(state)
vphase.axis = vp.vector(phase.real, phase.imag, 0)
vp.rate(1000)
So already, it sort of might make intuitive sense that half integer particles might have to do with fermions: there's that $-1$!
But that assumes some kind of connection between rotations and permutations in the tensor product.
After all, we've quantized spin-$\frac{1}{2}$ states as bosons and everything turned out fine. We can form permutation symmetric states, and swap the particles in them without picking up a $-1$. That only comes in if we were to rotate one of the particles individually, but of course, that would break the permutation symmetry, which suggests we can't act on just one boson, but have to act on them all symmetrically. Either way, the full state will pick up a global negative $-1$ due to the rotation.
In fact, the classic proofs of the spin statistics theorem involve invoking some high powered machinery from quantum field theory and special relativity in 3+1 dimensions, which we can't go into here.
Instead, however, let's look at an interesting construction thanks to Berry and Robbins, who approach the proof from a different angle. It involves the oscillator representation of spin, so it'll be a good example of that too!
So we have three concepts we want to combine: that of a swap in the tensor product, a rotation of spin, and an actual exchange of position. Take two objects $A$ and $B$, put them side to side, and swap them. Imagine you were $A$: what would you see? As $B$ circles around you, you start to see $B$'s underside, and then eventually $B$'s far side. In other words, $A$ would see a swap as: half a rotation of $B$, a rotation of $\pi$ radians. So let's ignore their actual position, and describe things just from $A$ and $B$'s point of view on each other. That leaves us with two things: we want to combine a $180^{\circ}$ rotation with a swap in the tensor product.
In order to do this, we're going to use our oscillator construction of spin, where each spin gets two oscillators. By the way, you can think of this like taking the "square root" of the spin. Being able to split a spin into two oscillators turns out to be very useful here.
Suppose we have two spins. So we have four oscillators in total, with annihilators: $a_{\uparrow}, a_{\downarrow}, b_{\uparrow}, b_{\downarrow}$. As we well know, we can upgrade our spin-$\frac{1}{2}$ $X, Y, Z$ operators to act on the pair of $a$ oscillators.
$X = \begin{pmatrix} a_{\uparrow}^{\dagger} & a_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} a_{\uparrow} \\ a_{\downarrow}\end{pmatrix} = \begin{pmatrix} a_{\uparrow}^{\dagger} & a_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} a_{\downarrow} \\ a_{\uparrow}\end{pmatrix} = a_{\uparrow}^{\dagger}a_{\downarrow} + a_{\downarrow}^{\dagger}a_{\uparrow}$
$Y = \begin{pmatrix} a_{\uparrow}^{\dagger} & a_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} \begin{pmatrix} a_{\uparrow} \\ a_{\downarrow}\end{pmatrix} = \begin{pmatrix} a_{\uparrow}^{\dagger} & a_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} -ia_{\downarrow} \\ ia_{\uparrow}\end{pmatrix} = -ia_{\uparrow}^{\dagger}a_{\downarrow} + ia_{\downarrow}^{\dagger}a_{\uparrow}$
$Z = \begin{pmatrix} a_{\uparrow}^{\dagger} & a_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \begin{pmatrix} a_{\uparrow} \\ a_{\downarrow}\end{pmatrix} = \begin{pmatrix} a_{\uparrow}^{\dagger} & a_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} a_{\uparrow} \\ -a_{\downarrow}\end{pmatrix} = a_{\uparrow}^{\dagger}a_{\uparrow} - a_{\downarrow}^{\dagger}a_{\downarrow}$
We can do the same with the $b$ oscillators:
$X = \begin{pmatrix} b_{\uparrow}^{\dagger} & b_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} b_{\uparrow} \\ b_{\downarrow}\end{pmatrix} = \begin{pmatrix} b_{\uparrow}^{\dagger} & b_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} b_{\downarrow} \\ b_{\uparrow}\end{pmatrix} = b_{\uparrow}^{\dagger}b_{\downarrow} + b_{\downarrow}^{\dagger}b_{\uparrow}$
$Y = \begin{pmatrix} b_{\uparrow}^{\dagger} & b_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} \begin{pmatrix} b_{\uparrow} \\ b_{\downarrow}\end{pmatrix} = \begin{pmatrix} b_{\uparrow}^{\dagger} & b_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} -ib_{\downarrow} \\ ib_{\uparrow}\end{pmatrix} = -ib_{\uparrow}^{\dagger}b_{\downarrow} + ib_{\downarrow}^{\dagger}b_{\uparrow}$
$Z = \begin{pmatrix} b_{\uparrow}^{\dagger} & b_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \begin{pmatrix} b_{\uparrow} \\ b_{\downarrow}\end{pmatrix} = \begin{pmatrix} b_{\uparrow}^{\dagger} & b_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} b_{\uparrow} \\ -b_{\downarrow}\end{pmatrix} = b_{\uparrow}^{\dagger}b_{\uparrow} - b_{\downarrow}^{\dagger}b_{\downarrow}$
But why stop there? We have four oscillators after all. We could consider spins built out of cross-oscillators, taking as pairs: the two $\uparrow$ oscillators (of $a$ and $b$ each):
$X = \begin{pmatrix} a_{\uparrow}^{\dagger} & b_{\uparrow}^{\dagger} \end{pmatrix} \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} a_{\uparrow} \\ b_{\uparrow}\end{pmatrix} = \begin{pmatrix} a_{\uparrow}^{\dagger} & b_{\uparrow}^{\dagger} \end{pmatrix} \begin{pmatrix} b_{\uparrow} \\ a_{\uparrow}\end{pmatrix} = a_{\uparrow}^{\dagger}b_{\uparrow} + b_{\uparrow}^{\dagger}a_{\uparrow}$
$Y = \begin{pmatrix} a_{\uparrow}^{\dagger} & b_{\uparrow}^{\dagger} \end{pmatrix} \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} \begin{pmatrix} a_{\uparrow} \\ b_{\uparrow}\end{pmatrix} = \begin{pmatrix} a_{\uparrow}^{\dagger} & b_{\uparrow}^{\dagger} \end{pmatrix} \begin{pmatrix} -ib_{\uparrow} \\ ia_{\uparrow}\end{pmatrix} = -ia_{\uparrow}^{\dagger}b_{\uparrow} + ib_{\uparrow}^{\dagger}a_{\uparrow}$
$Z = \begin{pmatrix} a_{\uparrow}^{\dagger} & b_{\uparrow}^{\dagger} \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \begin{pmatrix} a_{\uparrow} \\ b_{\uparrow}\end{pmatrix} = \begin{pmatrix} a_{\uparrow}^{\dagger} & b_{\uparrow}^{\dagger} \end{pmatrix} \begin{pmatrix} a_{\uparrow} \\ -b_{\uparrow}\end{pmatrix} = a_{\uparrow}^{\dagger}a_{\uparrow} - b_{\uparrow}^{\dagger}b_{\uparrow}$
As well as the two $\downarrow$ oscillators:
$X = \begin{pmatrix} a_{\downarrow}^{\dagger} & b_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} a_{\downarrow} \\ b_{\downarrow}\end{pmatrix} = \begin{pmatrix} a_{\downarrow}^{\dagger} & b_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} b_{\downarrow} \\ a_{\downarrow}\end{pmatrix} = a_{\downarrow}^{\dagger}b_{\downarrow} + b_{\downarrow}^{\dagger}a_{\downarrow}$
$Y = \begin{pmatrix} a_{\downarrow}^{\dagger} & b_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} \begin{pmatrix} a_{\downarrow} \\ b_{\downarrow}\end{pmatrix} = \begin{pmatrix} a_{\downarrow}^{\dagger} & b_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} -ib_{\downarrow} \\ ia_{\downarrow}\end{pmatrix} = -ia_{\downarrow}^{\dagger}b_{\downarrow} + ib_{\downarrow}^{\dagger}a_{\downarrow}$
$Z = \begin{pmatrix} a_{\downarrow}^{\dagger} & b_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \begin{pmatrix} a_{\downarrow} \\ b_{\downarrow}\end{pmatrix} = \begin{pmatrix} a_{\downarrow}^{\dagger} & b_{\downarrow}^{\dagger} \end{pmatrix} \begin{pmatrix} a_{\downarrow} \\ -b_{\downarrow}\end{pmatrix} = a_{\downarrow}^{\dagger}a_{\downarrow} - b_{\downarrow}^{\dagger}b_{\downarrow}$
Now suppose we evolve with $e^{-i\theta (Y_{a_{\uparrow}b_{\uparrow}} + Y_{a_{\downarrow}b_{\downarrow}})}$, in other words, perform a simultaneous $Y$ rotation of $\theta$ degrees to the two cross spins. If we watch what happens to the original spins during this $Y$ rotation, after a turn of $\pi$, the two spins will have swapped places in the tensor product. A turn of $2\pi$ returns them back to where they began. And so, we've unified rotations and permutations.
Moreover, if we consider doing this permutation/rotation to states with different $j$ values, trying out $\frac{1}{2}, 1, \frac{3}{2}, 2, \dots$, we find that: for half-integer spins, the state is multiplied by $-1$ after a swap, but not for integer spins. So that if we do a rotation/swap, half integer spins act like fermions, picking up a sign; but integer spins act like bosons, permuting without picking up a sign. That's hardly a complete proof, but it's suggestive!
import qutip as qt
import numpy as np
import math
import vpython as vp
##########################################################################################
def get_phase(q):
c = sum(q.full().T[0][::-1])
return np.exp(1j*np.angle(c))
def second_quantize_operator(O, a):
O = O.full()
terms = []
for i in range(len(a)):
for j in range(len(a)):
terms.append(a[i].dag()*O[i][j]*a[j])
return sum(terms)
def second_quantize_state(state, a):
v = state.full().T[0]
return sum([v[i]*a[i] for i in range(len(a))]).dag()
def second_quantize_spin_state(spin, a):
n = spin.shape[0]-1
j = (spin.shape[0]-1)/2.
v = spin.full().T[0]
terms = []
z, w = [a_.dag() for a_ in a]
for m in np.arange(-j, j+1, 1):
i = int(m+j)
terms.append(v[i]*(z**(n-i))*(w**i)*\
(math.sqrt(math.factorial(2*j)/\
(math.factorial(j-m)*math.factorial(j+m)))))
return sum(terms)
def upgrade(O, i, n):
ID = qt.identity(O.shape[0])
ID.dims = O.dims
return qt.tensor(*[O if i==j else ID for j in range(n)])
##########################################################################################
max_ex = 5
a = [upgrade(qt.destroy(max_ex), 0, 4), upgrade(qt.destroy(max_ex), 1, 4)]
b = [upgrade(qt.destroy(max_ex), 2, 4), upgrade(qt.destroy(max_ex), 3, 4)]
XYZ = [0.5*qt.sigmax(), 0.5*qt.sigmay(), 0.5*qt.sigmaz()]
XYZa = [second_quantize_operator(O, a) for O in XYZ]
XYZb = [second_quantize_operator(O, b) for O in XYZ]
XYZab = [second_quantize_operator(O, [a[0], b[0]]) for O in XYZ]
XYZba = [second_quantize_operator(O, [a[1], b[1]]) for O in XYZ]
Na = sum([a_.dag()*a_ for a_ in a])/2
Nb = sum([b_.dag()*b_ for b_ in b])/2
Nab = (a[0].dag()*a[0] + b[0].dag()*b[0])/2
Nba = (a[1].dag()*a[1] + b[1].dag()*b[1])/2
##########################################################################################
vac = qt.basis(a[0].shape[0])
vac.dims = [a[0].dims[0], [1]*len(a[0].dims[0])]
a_state = qt.rand_ket(2)
b_state = qt.rand_ket(2)
# try spin-1/2
state = second_quantize_spin_state(a_state, a)*second_quantize_spin_state(b_state, b)*vac
# try spin-1
#state = second_quantize_spin_state(a_state, a)*second_quantize_spin_state(b_state, b)*second_quantize_spin_state(a_state, a)*second_quantize_spin_state(b_state, b)*vac
# try tensor state
#ab_state = qt.bell_state("00")
#ab_tensor = [a[0]*b[0], a[0]*b[1], a[1]*b[0], a[1]*b[1]]
#state = second_quantize_state(ab_state, ab_tensor)*second_quantize_state(ab_state, ab_tensor)*vac
state = state.unit()
##########################################################################################
scene = vp.canvas(background=vp.color.white, width=1000, height=800)
vA = vp.sphere(pos=vp.vector(-1.5, 0, 0), radius=qt.expect(Na, state), color=vp.color.blue, opacity=0.3)
vB = vp.sphere(pos=vp.vector(1.5, 0, 0), radius=qt.expect(Nb, state), color=vp.color.red, opacity=0.3)
vAB = vp.sphere(pos=vp.vector(0, 1.5, 0), radius=qt.expect(Nab, state), color=vp.color.yellow, opacity=0.3)
vBA = vp.sphere(pos=vp.vector(0, -1.5, 0), radius=qt.expect(Nba, state), color=vp.color.green, opacity=0.3)
vA_arrow = vp.arrow(pos=vA.pos,\
axis=vp.vector(*[qt.expect(O, state) for O in XYZa]))
vB_arrow = vp.arrow(pos=vB.pos,\
axis=vp.vector(*[qt.expect(O, state) for O in XYZb]))
vAB_arrow = vp.arrow(pos=vAB.pos,\
axis=vp.vector(*[qt.expect(O, state) for O in XYZab]))
vBA_arrow = vp.arrow(pos=vBA.pos,\
axis=vp.vector(*[qt.expect(O, state) for O in XYZba]))
vA_origA = vp.sphere(pos=vA.pos+vA_arrow.axis, color=vp.color.blue, radius=0.15)
vA_origB = vp.sphere(pos=vB.pos+vA_arrow.axis, color=vp.color.blue, radius=0.15)
vB_origA = vp.sphere(pos=vA.pos+vB_arrow.axis, color=vp.color.red, radius=0.15)
vB_origB = vp.sphere(pos=vB.pos+vB_arrow.axis, color=vp.color.red, radius=0.15)
phase = get_phase(state)
initial_vphase = vp.arrow(pos=vp.vector(0,0,0),color=vp.color.yellow,\
axis=0.5*vp.vector(phase.real, phase.imag, 0))
vphase = vp.arrow(pos=vp.vector(0,0,0),\
axis=0.5*vp.vector(phase.real, phase.imag, 0))
##########################################################################################
dt = 0.0005
U = (-1j*(XYZba[1]+XYZab[1])*dt).expm()
T, t = 0, 0
#while T < np.pi:
while True:
state = U*state
vA_arrow.axis = vp.vector(*[qt.expect(O, state) for O in XYZa])
vB_arrow.axis = vp.vector(*[qt.expect(O, state) for O in XYZb])
vAB_arrow.axis = vp.vector(*[qt.expect(O, state) for O in XYZab])
vBA_arrow.axis = vp.vector(*[qt.expect(O, state) for O in XYZba])
vA.radius = qt.expect(Na, state)
vB.radius = qt.expect(Nb, state)
vAB.radius = qt.expect(Nab, state)
vBA.radius = qt.expect(Nba, state)
#print(qt.entropy_vn(state.ptrace(0,1)), qt.entropy_vn(state.ptrace(2,3)))
T += dt
t += 1
if t % 100 == 0:
phase = get_phase(state)
vphase.axis = 0.5*vp.vector(phase.real, phase.imag, 0)
vp.rate(1000)
Some references:
Back and forth from Fock space to Hilbert space: a guide for commuters
Indistinguishability for Quantum Particles: Spin, Statistics and the Geometric Phase
Group Representations and the Quantum Statistics of Spin
The Geometry of Classical Particles
Atiyah's Edinburgh Lectures on Geometry, Analysis and Physics