In this article we are going to simulate different kinds of indistinguishable particles and introduce that concept briefly in case if you are not yet familiar with them. We will also write few words about the second quantization formalism in 11D to be able to formally represent quantum states involving those particles. For fermions we will use Jordan-Wigner transformation, for bosons we will use native operators already available in QuTip. In the end you should be able define operators that act on fermions or bosons in QuTip, we will test their properties using pytest testing framework and eventually as we conclude, we will briefly discuss the idea of test driven science.

For a starter, brief introduction to indistinguishable particles. Let ∣ψαψβ>\left\vert \psi_\alpha \psi_\beta \right> be a quantum state representing a system holding two indistinguishable particles α\alpha and β\beta. Their indistinguishability implies that

∣<ψβψα∣ψαψβ>∣2=1.\begin{aligned} \left \vert \left< \psi_\beta \psi_\alpha \vert \psi_\alpha \psi_\beta \right> \right \vert^2 &= 1. \end{aligned}

Note that particles α\alpha and β\beta are indistinguishable from measurement perspective and what vanishes during quantum measurement is the phase factor. In such case exchanging particles should allow a phase factor without violating their indistinguishability. To express that let UU be an operator that exchanges two particles

U∣ψαψβ>=eiϕ∣ψβψα>.\begin{aligned} U \left \vert \psi_\alpha \psi_\beta \right> &= e^{i\phi} \left \vert \psi_\beta \psi_\alpha \right>. \end{aligned}

The value of ϕ\phi characterizes groups of indistinguishable particles. For bosons phase factor is ϕ=0\phi = 0 and for fermions it is ϕ=π\phi = \pi. Particles that have any different value of ϕ\phi are commonly referred to as anyons. Let us change the notation a bit to something more workable computationally, for this let us use the second quantization formulation which allows us to describe quantum many-body states from perspective of occupation number. For this we consider the Fock state basis, consisting of single particle states occupancies. For example, a Fock state for a quantum 44-body system in a vacuum would be ∣0000>\left \vert 0000 \right> and if last site is occupied by a single particle we would write it as ∣0001>\left \vert 0001 \right> and so on.

Such Fock states are manipulated using creation and annihilation operators. For fermions we write them as aj†a^\dagger_j and aja_j. For bosons bj†b^\dagger_j and bjb_j. The subscript jj represents the lattice site index. For example we could place fermions on the edges using a1†a4†∣0000>=∣1001>a^\dagger_1 a^\dagger_4\left \vert 0000 \right> = \left \vert 1001 \right>. Values in the ket represent occupancy number. Note that this notation does not really distinguish what occupies the site so let us just assume those were fermionic sites and from now on let us try to just define if site is meant for bosons or fermions. We should try to avoid mixing them up. The properties of indistinguishable particles (i.e. what phase factor their exchange pulls out) is entirely described by their operators, thus for fermions we have

{an,an′}=1{an,an′†}=δnn′an2=0\begin{aligned} \{a_n, a_{n^\prime}\} &= 1 \\ \{a_n, a_{n^\prime}^\dagger\} &= \delta_{nn^\prime} \\ a_n^2 &= 0 \end{aligned}

and for bosons

[bn,bn′]=0[bn,bn′†]=δnn′.\begin{aligned} [b_n, b_{n^\prime}] &= 0 \\ [b_n, b_{n^\prime}^\dagger] &= \delta_{nn^\prime}. \end{aligned}

the above properties are true in the infinite-dimensional Hilbert space. Here we consider numerical simulation so our Hilbert space is fixed to hold maximum NpN_p bosons per site and commutation relation becomes

[bn,bn′]=0[bn,bn′†]=δnn′[1−Np+1Np!(bn)pN(bn†)pN]\begin{aligned} [b_n, b_{n^\prime}] &= 0 \\ [b_n, b_{n^\prime}^\dagger] &= \delta_{nn^\prime} [1 - \frac{N_p + 1}{N_p!} (b_n)^N_p (b_n^\dagger)^N_p] \end{aligned}

from (Somma, Ortiz, Knill, & Gubernatis, 2003).

Where [A,B]=AB−BA[A, B] = AB - BA is the commutator and {A,B}=AB+BA\{A, B\} = AB + BA is the anticommutator. Also note how bosonic operators do not square to zero!

To use those operators in our simulations we need to translate them into spin-12\frac{1}{2}, we do not worry about bosonic operators because QuTip already provides creation and annihilation operators for NN-level Fock basis for bosons. For fermions (unfortunately) we cannot be so lazy and we need to make them ourselves and for this we need the famous Jordan-Wigner transformation

aj=(∏k=1j−1σkα)(σjβ+iσjγ)aj†=(∏k=1j−1σkα)(σjβ−iσjγ)\begin{aligned} a_j &= (\prod_{k=1}^{j-1} \sigma^\alpha_k)(\sigma^\beta_j + i \sigma^\gamma_j) \\ a_j^\dagger &= (\prod_{k=1}^{j-1} \sigma^\alpha_k)(\sigma^\beta_j - i \sigma^\gamma_j) \end{aligned}

I do not like to be attached to any particular basis so I generalized the operators for Jordan-Wigner transformation to Pauli operators labeled σα\sigma^\alpha, σβ\sigma^\beta and σγ\sigma^\gamma. The choice of those operators would affect the spin-12\frac{1}{2} representation of the vacuum state and could potential mix-up the creation and annihilation operators. The choice of σα\sigma^\alpha determines how vacuum is represented, vacuum state would be the state that has +1+1 eigenvalue of the σα\sigma^\alpha-operator. The occupied state would be one that has −1-1 eigenvalue. Regarding the choice of remaining operators, if σβσγ=iσα\sigma^\beta \sigma^\gamma = i\sigma^\alpha then aa and a†a^\dagger will represent annihilation and creation respectively (as expected). On the other hand choosing σβσγ=−iσα\sigma^\beta \sigma^\gamma = -i\sigma^\alpha would swap their interpretations.

Let us implement the Jordan-Wigner transformation that works for an arbitrary choice of those operators for a NN-body system. First some spin-12\frac{1}{2} operators.

import numpy as np
from qutip import qeye, sigmax, sigmay, sigmaz, tensor

def Is(i): return [qeye(2) for j in range(0, i)]
def Sx(N, i): return tensor(Is(i) + [sigmax()] + Is(N - i - 1))
def Sy(N, i): return tensor(Is(i) + [sigmay()] + Is(N - i - 1))
def Sz(N, i): return tensor(Is(i) + [sigmaz()] + Is(N - i - 1))

Also, as usual, product, sum and power operators that act on generic Object-type which is very useful for defining quantum mechanical operators using QuTip.

def osum(lst): return np.sum(np.array(lst, dtype=object))

def oprd(lst, d=None):
    if len(lst) == 0:
        return d
    p = lst[0]
    for U in lst[1:]:
        p = p*U
    return p

def opow(op, N): return oprd([op for i in range(N)])

Using these we can defined Jordan-Wigner operators.

def a_(N, n, Opers=None):
    Sa, Sb, Sc = Sz, Sx, Sy
    if Opers is not None:
        Sa, Sb, Sc = Opers
    return oprd([Sa(N, j) for j in range(n)], d=I(N))*(Sb(N, n) + 1j*Sc(N, n))/2.

def ad(N, n, Opers=None):
    Sa, Sb, Sc = Sz, Sx, Sy
    if Opers is not None:
        Sa, Sb, Sc = Opers
    return oprd([Sa(N, j) for j in range(n)], d=I(N))*(Sb(N, n) - 1j*Sc(N, n))/2.

We could have just use a_.dag() for the definition of ad as they only differ by a sign but I find it beneficial to sometimes follow the definitions strictly even if they differ only by a sign. Benefit of that is it helps to correlate things for those who read this article and see those things for the first time. We also need identity operator, but not just any identity… We need identity that matches the QuTip’s quantum object type of Hilbert space of NN 22-level particles.

def I(N): return osum([Sz(N, i)*Sz(N, i) for i in range(N)])/N

Lets write a little test script that checks basic properties of fermions and bosons. It must cover all possible choices of Jordan-Wigner operators (work in any basis) for fermions and for bosons we expect it to cover cases of different number of possible levels. Lets start by making the necessary imports.

import pytest
import math
import itertools

from qm import Sx, Sy, Sz
from qm import a_, ad, b_, bd, I
from qm import opow
from qm import commutator, anticommutator

Now parameters, we want to cover cases of N=2…4N = 2 \dots 4 particles and in case of bosons Np=2…4N_p = 2 \dots 4 levels. Using the itertools library we prepare all possible permutations of spin operators for Jordan-Wigner and provide them with string type label for test naming.

Ns = [2, 3, 4]
Nps = [2, 3, 4]
jws = itertools.permutations([(Sx, 'X'), (Sy, 'Y'), (Sz, 'Z')])

Prepare products of those parameters to cover all cases for fermions and bosons as well as pre-generate the test names. Fermion tests will be dynamically labeled using fermionTestName function that accepts parameters and returns test case name as a string.

def fermionTestName(param):
    N, jw = param
    (_, la), (_, lb), (_, lc) = jw
    return 'N={0},JW={1}'.format(str(N), la + lb + lc)

fermion_params = list(itertools.product(Ns, jws))
fermion_names = [fermionTestName(param) for param in fermion_params]

boson_params = list(itertools.product(Ns, Nps))
boson_names = ['N={0},Ns={1}'.format(str(N), str(Ns)) for N, Ns in boson_params]

Finally, let us implement test for fermion operators. It is a literal implementation of fermion operator properties stated above. Note how zero operator is prepared by multiplying identity by constant zero

@pytest.mark.parametrize('N,jw', fermion_params, ids=fermion_names)
def testFermions(N, jw):
    (Sa, _), (Sb, _), (Sc, _) = jw
    Opers = Sa, Sb, Sc
    zero = 0.*I(N)
    # test all the pairs
    for n in range(N):
        a_n = a_(N, n, Opers=Opers)
        adn = ad(N, n, Opers=Opers)
        for np in range(N):
            a_np = a_(N, np, Opers=Opers)
            adnp = ad(N, np, Opers=Opers)
            assert anticommutator(a_n, a_np) == zero
            if n == np:
                assert anticommutator(a_n, adnp) == I(N)
            else:
                assert anticommutator(a_n, adnp) == zero
        assert a_n*a_n == zero
        assert adn*adn == zero

we proceed in the same way to design test for bosons, following the finite Hilbert space boson commutation relations stated above

@pytest.mark.parametrize('N,Np', boson_params, ids=boson_names)
def testBosons(N, Np):
    zero = 0.*b_(N, Np, 0)*bd(N, Np, 0)
    # test all the pairs
    for n in range(N):
        b_n = b_(N, Np, n)
        bdn = bd(N, Np, n)
        for np in range(N):
            b_np = b_(N, Np, np)
            bdnp = bd(N, Np, np)
            # test anticommutation properties
            assert commutator(b_n, b_np) == zero
            LHS = commutator(b_n, bdnp)
            RHS = zero
            if n == np:
                NpF = math.factorial(Np)
                RHS = (1. - ((Np+1)/NpF)*opow(bdn, Np)*opow(b_n, Np))
            assert LHS == RHS

$ python -m pytest tests.py -v
============================= test session starts ==============================
platform darwin -- Python 3.5.3, pytest-4.2.0, py-1.8.0, pluggy-0.12.0 -- /Users/marek/.pyenv/versions/3.5.3/bin/python
cachedir: .pytest_cache
metadata: {'Platform': 'Darwin-18.6.0-x86_64-i386-64bit', 'Python': '3.5.3', 'Packages': {'pluggy': '0.12.0', 'py': '1.8.0', 'pytest': '4.2.0'}, 'Plugins': {'steps': '1.6.4', 'timeout': '1.3.3', 'html': '1.22.1', 'metadata': '1.8.0', 'dynamodb': '1.2.0', 'flaky': '3.6.1'}}
rootdir: /Users/marek/Development/bloggists/jordan-wigner, inifile:
plugins: steps-1.6.4, timeout-1.3.3, dynamodb-1.2.0, flaky-3.6.1, html-1.22.1, metadata-1.8.0
collected 27 items                                                             

tests.py::testFermions[N=2,JW=XYZ] PASSED                                [  3%]
tests.py::testFermions[N=2,JW=XZY] PASSED                                [  7%]
tests.py::testFermions[N=2,JW=YXZ] PASSED                                [ 11%]
tests.py::testFermions[N=2,JW=YZX] PASSED                                [ 14%]
tests.py::testFermions[N=2,JW=ZXY] PASSED                                [ 18%]
tests.py::testFermions[N=2,JW=ZYX] PASSED                                [ 22%]
tests.py::testFermions[N=3,JW=XYZ] PASSED                                [ 25%]
tests.py::testFermions[N=3,JW=XZY] PASSED                                [ 29%]
tests.py::testFermions[N=3,JW=YXZ] PASSED                                [ 33%]
tests.py::testFermions[N=3,JW=YZX] PASSED                                [ 37%]
tests.py::testFermions[N=3,JW=ZXY] PASSED                                [ 40%]
tests.py::testFermions[N=3,JW=ZYX] PASSED                                [ 44%]
tests.py::testFermions[N=4,JW=XYZ] PASSED                                [ 48%]
tests.py::testFermions[N=4,JW=XZY] PASSED                                [ 51%]
tests.py::testFermions[N=4,JW=YXZ] PASSED                                [ 55%]
tests.py::testFermions[N=4,JW=YZX] PASSED                                [ 59%]
tests.py::testFermions[N=4,JW=ZXY] PASSED                                [ 62%]
tests.py::testFermions[N=4,JW=ZYX] PASSED                                [ 66%]
tests.py::testBosons[N=2,Ns=2] PASSED                                    [ 70%]
tests.py::testBosons[N=2,Ns=3] PASSED                                    [ 74%]
tests.py::testBosons[N=2,Ns=4] PASSED                                    [ 77%]
tests.py::testBosons[N=3,Ns=2] PASSED                                    [ 81%]
tests.py::testBosons[N=3,Ns=3] PASSED                                    [ 85%]
tests.py::testBosons[N=3,Ns=4] PASSED                                    [ 88%]
tests.py::testBosons[N=4,Ns=2] PASSED                                    [ 92%]
tests.py::testBosons[N=4,Ns=3] PASSED                                    [ 96%]
tests.py::testBosons[N=4,Ns=4] PASSED                                    [100%]

========================== 27 passed in 4.72 seconds ===========================

We tested commutation / anticommutation properties of different kinds of indistinguishable particles. We defined QuTip operators for fermions using Jordan-Wigner operators. We have not simulated any physical system although using those operators we can start working with different physical systems in 1D geometry. Jordan-Wigner transformation although primarily designed for 1D systems, it can be extended to 2D (Azzouz, 1993).

I am strongly attached to the idea of testing operator properties. Writing pytest tests takes huge amount of my daily research work but it is a good time investment. I am a huge advocate of test-driven development (TDD) and I brought those habits with me to research work as I career switched from software engineering. Research is primarily exploration of the unknown. Doing research requires exploring ideas which have never been tried before. In case of numerical simulations it remains true. I do not believe it would be possible to get right results if at least one operator is not working properly. I observe people working with tools like Mathematica notebooks or Jupyter, often looping around being blocked on same bug resulting from some old state being preserved as one of the cells was not re-run after changes were made. I do large part of my numerics with pytest, whenever I test a hypothesis literally everything is being tested so that if I temporarily change something to try an idea test and if that breaks anything then pytest will catch it immediately.

For the complete source code please consult the gist repository.

  1. Somma, R. D., Ortiz, G., Knill, E. H., & Gubernatis, J. (2003). Quantum simulations of physics problems. (E. Donkor, A. R. Pirich, & H. E. E. Brandt, Eds.)Quantum Information and Computation. 10.1117/12.487249
  2. Azzouz, M. (1993). Interchain-coupling effect on the one-dimensional spin-1/2 antiferromagnetic Heisenberg model. Physical Review B, 48(9), 6136–6140. 10.1103/physrevb.48.6136

[Back to Top]