CITE THIS NOTEBOOK: QCPITT: quantum computer simulation packages by Frank Tabakin. Wolfram Community OCT 5 2022.
The quantum computer simulation codes QDENSITY, QCWAVE, BTSYSTEM, and CIRCUITS are presented as a combined version: QCPITT. These Mathematica codes provide a pedagogic tool for understanding the basic ideas of a quantum computer and in addition provide a tool for studying the advantages and vulnerabilities of a quantum computer. In QDENSITY, a density matrix description is used to describe quantum systems; whereas, QCWAVE uses a quantum state description. BTSYSTEM is used not only for qubits (Binary), but also for qutrits (Triplet) states and can be extended to qudrits (more than 3 levels). CIRCUITS provides a graphical display of quantum circuits. These mathematical codes were developed at the University of Pittsburgh during the 1990s and have been published and posted, accompanied by many detailed tutorials.
​
These codes are designed to be simple, open and approachable. Overly sophisticated definitions are avoided, in favor of readily understood, revisable and expandable versions. A modular, instructive approach is used to generate the basic elements that make up a quantum circuit. The QCPITT package includes: Qdensity.m, QCWave.m, Circuits.m, and BTSystem.m, which contain the tools needed in quantum circuits, e.g. multi-qubit kets, projectors, gates, etc. Selected examples of the basic commands are presented here and in tutorial notebooks that serve as a full guide and templates. Application is made to a variety of relevant cases, including Teleportation, Quantum Fourier Transform, Grover’s search and Shor’s algorithm, in separate notebooks: QFT.nb, Teleportation.nb, Grover.nb and Shor.nb, where each algorithm is explained in detail. A Mathematica palette containing some commands is also included: QCPITTpalette.nb. Recent developments include: (1) application to qutrits and qudits; (2) extension to closed and open system dynamics using a Lindblad/Beretta formulation, (3) studies of one and two qubit gate dynamics. and (4) multi-dimensional tunnelling (in preparation). Future goals are delineated.

Introduction

QCPITT is a collection of Mathematica packages that simulate a quantum computer on a classical one. This simulation closely follows the procedures that would occur in an actual quantum computer, wherein quantum interference, superposition and entanglement play key roles. Although emphasis is on using a density matrix, QCPITT also handles state vectors. The package comes with numerous tutorials that enable the user to build any quantum circuit out of conventional quantum gates. These notebooks illustrate the most relevant quantum algorithms: teleportation, searching and factoring algorithms in a self contained way. Extensions to both binary (qubit ) and triplet (qutrit) systems are included.

Installation

Note: To download the QCPITT package go to:
​https://sites.pitt.edu/~tabakin/QDENSITY/​
Once installed on your computer, the file QCPITTessay will be interactive.
Version: 6.0-April 2022--For Mathematica 13
Click QCPITT to download and then unzip QCPITT.zip. All the quantum computing simulation codes will start by loading the Qdensity, QCWave , BTSystem and QCircuits packages. To make the packages fully available put the QDENSITY folder directory for example to the location: “/Users/youraccount/Library/Mathematica/Applications/” on your user account, as deduced for a MacOS below:
In[]:=
FileNameJoin[{$UserBaseDirectory,"Applications","QDENSITY"}]
Out[]=
/Users/franktabakin/Library/Mathematica/Applications/QDENSITY
To check proper installation type in:
SystemOpen[“/Users/franktabakin/Library/Mathematica/Applications/QDENSITY”]
To access the tutorials on a range of examples click here Click
TUTORIALS
. These and a lots more are also included directly in the full download.
Once the Context is set with Mathematica running , load Qdensity, QCWave and, if required, Circuits and BTSystem using:
Off[General::shdw];​​Needs["QDENSITY`Qdensity`"];​​Needs["QDENSITY`QCWave`"];​​Needs["QDENSITY`BTSystem`"];​​Needs["QDENSITY`Circuits`"];
In[]:=
?introBT
Out[]=
Symbol
The quantum states and density matix of quantum systems which consist ofmixtures of qubits and qutrits are provided in this notebook.Partial trace, and othe commands are also provided and are illustrated.
To list all commands type in: ?QDENSITY`Qdensity`* , ?QDENSITY`QCWave`*, ?QDENSITY`BTSysten`*, ?QDENSITY`Circuits`*. To list commands starting with a particular letter such as "C" in each package type:
In[]:=
?QDENSITY`Qdensity`C*
Out[]=
QDENSITY`Qdensity`
CNOT
Comm
CompPTr
ControlledY
CPHASE
Coeff
CompMeasureKet
ControlledX
ControlledZ
CPHASEV
In[]:=
?QDENSITY`QCWave`D*
Out[]=
QDENSITY`QCWave`
DForm
DFormA
In[]:=
?QDENSITY`Circuits`C*
Out[]=
QDENSITY`Circuits`
C2
C3X
CircuitW
CNPHI
CNTRR3
CNTRRnm2
CNTRX
CX
C2X
Circuit
CNOTC
CNTRG
CNTRRn
CNTRSM
CNTRY
CY
C3
CircuitG
CNOTS
CNTRR2
CNTRRnm1
CNTRSP
CNTRZ
CZ
In[]:=
?QDENSITY`BTSystem`B*
Out[]=
QDENSITY`BTSystem`
B
BL
BraT
BTbit
BTop
BTOp1
BTop2
BTOp2

Single-Qubits

Single-qubit states

As an introduction to using the QCPITT package, let us start by setting up the single qubit states, first in the form of Dirac state vectors and then as a density matrix. A qubit is a simple quantum system with two discrete quantum states. For example, this can be a spin 1/2 system that could point in the ± z-directions. Or it could be the two circular polarization states of a photon. We denote these two possibilities as a |0〉 or a |1〉 state. In a classical computer these are distinct, separate and measurable states, which are physical realizations for binary counting and are called bits. For a quantum system, the act of measurement disturbs the system and therefore the best description is that the quantum system is simultaneously in the |0〉 and |1〉 state with associated probability amplitudes. This is called a qubit and is described as a linear combination or superposition. There are two complex probability amplitudes,
a
and
b
, which yield the probability for the state |0〉 as
*
a
a
and for |1〉 as
*
b
b
. Using various forms along with the DForm to display a Dirac ket we have
In[]:=
Clear[a,b]​​|ψ〉=a|0〉+b|1〉​​ψ=aKet[0]+bKet[1]​​DForm[%%]​​DForm[%%]​​Ket[0]==|0〉&&Ket[1]==|1〉
Out[]=

a
b

Out[]=

a
b

Out[]=
+ (a) |0 > + (b) |1 >
Out[]=
+ (a) |0 > + (b) |1 >
Out[]=
True
A single qubit state is also written above using our equivalent Ket and Bra notation. Set the cell option “convert to” -> standard.
We normalize the state and set
{a
*
a
+b
*
b
}1
, which is a statement that we have the system and all of the system in hand. The equation below asserts that the two states are distinct
In[]:=
Bra[1].Ket[0]=={{0}}
Out[]=
True
The above Ket and Bra state vectors can also be displayed in Dirac notation using our DForm and DFormA “Dirac Form” commands. DForm is contained in QCWAVE. It also prints out results in Dirac Form:
In[]:=
aKet[0]+bKet[1]​​a|0〉+b|1〉==%​​DForm[aKet[0]+bKet[1]]
Out[]=

a
b

Out[]=
True
Out[]=
+ (a) |0 > + (b) |1 >
In[]:=
aBra[0]+bBra[1]​​a〈0|+b〈1|​​DFormA[aBra[0]+bBra[1]]
Out[]=
(
a
b
)
Out[]=
(
a
b
)
Out[]=
+ ( a ) < 0 | + ( b ) < 1 |
Here we introduce the row (bra form or adjoint/dual space) and the column (ket form) to
complete Dirac's humorous label of the overlap 〈0 | 1〉 as a "bracket." The advantage of this brief notation is that it handles superposition, normalization, orthogonality and completeness with grace, with every step having a carefully designed physical quantum purpose
In[]:=
Ket[0]==
1
0
&&Ket[1]==
0
1
​​Bra[0]==(
1
0
)&&Bra[1]==(
0
1
)
Out[]=
True
Out[]=
True
For a qubit pointing in the direction specified by the Euler angles α  (θ,ϕ), we have:
In[]:=
ψ
α
=Ket1[θ,ϕ]​​DForm[
ψ
α
]
Out[]=
-
ϕ
2

Cos
θ
2

ϕ
2

Sin
θ
2

Out[]=
+ (
-
ϕ
2

Cos
θ
2
) |0 > + (
ϕ
2

Sin
θ
2
) |1 >
This state can also be obtained by a rotation.
In[]:=
RotZ[ϕ].RotY[θ].Ket[0]​​Ket1[θ,ϕ]==%
Out[]=
-
ϕ
2

Cos
θ
2

ϕ
2

Sin
θ
2

Out[]=
True
A general Wigner rotation with Euler angles θ, ϕ, and γ is
In[]:=
WR[ϕv_,θv_,γv_]:=RotZ[ϕv].RotY[θv].RotZ[γv]​​WRA[ϕv_,θv_,γv_]:=RotZ[-γv].RotY[-θv].RotZ[-ϕv]​​WR[ϕ,θ,γ].WRA[ϕ,θ,γ]//FullSimplify;​​%==s[0]​​WR[ϕ,θ,γ]
So we recover Ket1. This is called a spinor rotation--a spin 1/2 rotation, where the Pauli operators are the generators of the unitary Wigner rotations.
This can be used to rotate states, operators and/or density matrices to expose any useful symmetries.
​
If we need to find the probability amplitude of the |0〉 or |1〉 state in Ket1, we form the "bracket" or overlap or projection as
This example shows how to pick out the desired probability amplitude, calculate the probability and after repeated "measurements" confirm that the full system has yielded a net probability of 1. Soon we will deal with an extension of this idea where we get that selected part of a state vector in hand with its associated amplitude. This is called a projection operator.

Pauli Spin Operators

Other properties follow:
We define raising and lowering operators that jump up from an assumed lower level
|0> to |1> (sp) or fall down from |1> to |0> (sm). Indeed, one could start from this operator, define all of its desired attributes and then find you have again developed the Pauli operators
The Pauli matrices form a basis for any trace 1 Hermitian matrix. Later we will see that a density matrix for one qubit is a trace-1, Hermitian matrix and can be expressed in terms of a real polarization vector {Px, Py, Pz} of magnitude P. An Hermitian matrix has real eigenvalues.
The Wigner rotation WR[ϕ, θ ,γ] can be used to rotate a spin operator to a new direction:

Single Qubit Projection Operators

For the |1〉 state, we have
These Ket and Bra state vectors are directly displayed in Dirac notation (set the cell option “convert to” -> standard)
These projection operators satisfy
Additional projection operators are also available in the package:

Density Matrix

And the ensemble average energy is

Density Matrix Examples

or by using the ρ[ψ] command (seen below to be Hermitian)
And we can evaluate the trace and purity of the above density matrix.
For a normalized state, these are both equal to 1 . The entropy is zero in line with this being a highly ordered ensemble with all spins aligned. ( We use EnTropy to avoid the MM Entropy command.)
The above is an example of a pure state; one that satisfies ρ . ρ = ρ
In papers associated with these codes, the general definition of a density matrix is presented, along with a discussion of a pure state, a mixed state, and for two or more qubits an entangled state.
The package also contains definitions for | + 〉 and | - 〉 states, i.e. for states in the plus or minus x-direction :
Scalar products are simple to compute; for instance, if we have a general state in the computational basis, we can compute its scalar product with the | + 〉 or | - 〉 operators (KetX[0] and KetX[1]). We can also get the Dirac form for the adjoint or Bra space using DFormA :
Now compute scalar product <state1 | + 〉 and <state1 | - 〉

The  and  unitary operators (* as used in qiskit *)

The following unitary operators  and  are used extensively in Qiskit. They are basically rotations about the z-axis, with a phase change on just the Ket[1] part. Note that  is the Matrix square root of s[3]. The properties are:
The phase changes produced are as follows:
The density matrix changes by rotating the x-y plane polarizations :

General 0ne-qubit Density Matrix

From the above discussion, we see that a one-qubit density matrix is a Hermitian operator with unit trace and as mentioned earlier can be described in terms of the Pauli matrices. It has positive definite eigenvalues and its observables stay within Bloch spheres.

Bell states

The Bell (two qubit) and the GHZ (three qubit) states are important for understanding entanglement.
They can all be constructed by building a Bell operator from a Hadamard and then a CNOT gate.
Then full 4 Bell basis states are then:
The following steps illustrate how to convert MM to python code, as part of a future plan to build qiskit code for
submission to a real quantum computer.
These Bell states are of great importance in quantum circuits. Check, for instance, the Teleportation circuit in the tutorials.
Then the density matrices ρBell for these Bell states are :
Bell states have zero vector polarizations and diagonal tensor polarization.
Bell states are the prime example of an entangled state and are often used to gauge a states entanglement. Entanglement is a correlation between spins that can not be produced by a product, a mixed or any classical procedure. It is a correlation that arises only in the quantum theory. It is a key resource for quantum computing. See the tutorial on Bell states.

GHZ states

Another type of entangled states are the GHZ states which can be prepared by having the following operator act on a state with three qubits :
Then the basic GHZ state is:
The density matrix for this GHZ state is:
Which is mathematically equivalent to building the state and then its density matrix as:
Then GHZ observables for this case are then:
The 8 distinct GHZ basis states are then:

Subsystems, Partial Traces

The density matrix of a specific subsystem can be computed easily by taking partial traces of the density matrix of the full system, e.g. we will build the 4 qubit state |0111〉 and then recover the state of either the subsystem (1) or the subsystems (24). For example, starting with a four qubit state |0111〉, its density matrix is constructed and then partial traces taken:
We can check basic properties of these density matrices, such as their length, trace and Purity:
Similarly we can produce a more complicated density matrix that is not pure,
It is interesting to show the density matrix for each qubit taking partial traces. These subspace density matrices can be thought of as the density matrices of one qubit in the environment of the remaining 2 qubits:
Pauli spin matrices can be used to build the more general 2 and n-qubit density matrices, e.g. for a general 2 qubit state the density matrix is:
Sometimes, a random one-qubit density matrix or a random 2 qubit density matrix is needed. Being Hermitian, it has real eigenvalues, which sum to one and are positive definite. For a density matrix, the system described is always pure in contrast to random state vectors. (The procedure used to build the random 2 qubit density matrix can be easily generalized to the n-qubit case.)
​
Consider a random 2 qubit density matrix.
Now the subsystems are also not pure states, for instance, the subsystem corresponding to qubit number 2 of the two qubits density matrix is obtained by performing the partial trace of qubit 1,
Use RandomQubitN[n] for n qubits
Analogous functions for one-qubit Ket and Bra exist. These are however pure states.

Qubit Dynamics

where ωAL and ωBL are the Larmor frequencies for qubit A and B; Planck's constant is set to 1. For example, for ωAL= 100 there is a lower level at - 50 and an upper level at + 50. These are the two level of qubit A. How does a qubit move in this case? Clearly the polarization vector and the associated density matric are time - dependent and the density matrices are of the general form below. Each density matrix has unit trace and is Hermitian for a real polarization vector.
The above coupled matrix equation can be solved using NDSolve with an initial condition of
ρB[0] or ρB[0]. Or one can deduce specific dynamic equations for the polarization vectors using
We adopt the second approach because many other attributes, such as energy, power and entropy, can be evaluated once ρA[t] is known. A sample case is shown for a pure Larmor precession case (ω2=0), followed by a Rabi case at resonance:
Top plots show the time-dependence of the polarization. Next as above but with rotated uniform field:
Next a Rabi driving field at resonance w=wL.
Next a Rabi at resonance w=wL with a rotated Hamiltonian and a rotated initial qubit:
The above Rabi resonance case simulates the action of a NOT gate and flips the spin.
Now an off-resonance w=1.4 wL Rabi oscillation case:
That off-resonance Rabi oscillation shows the result is much closer to the non Rabi result, i,e, it is off-resonance and does not spin flip.
Once we have a dynamic density matrix, it is simple to examine other dynamic aspects, such as the energy and power evolution and changes in entropy
The on-resonance single qubit Rabi oscillation case is a simple dynamic example of a spin flip or NOT gate. This can be simple extended to a phase change, or an Hadamard or any other spin rotation gate. In addition a Lindblad/Beretta term can be added to study thermal bath, entropy, dissipation, noise and measurement processes for a single qubit. An example of such a study is in Annals of Physics Volume 383, 2017, 33 link. Two qubit gates can also be implemented following the same procedures with QCPITT; a publication is in preparation.

Rotations

The general rotation of a qubit state is given by the Wigner rotation WR[ϕ, θ ,γ ].

Quantum Gates

In the density matrix language any operator Ω acting on the density matrix of the system ρ , produces a new density matrix ρ’

a) Single qubit operations: Hadamard and rotations

or the script form:
If we have 3 qubits and want to Hadamard the second one,
which is equivalent to
If we want to apply Hadamards as above to an initial system of |000〉, we construct the initial density matrix and transform to the new density matrix, also we check the subsystems 2 and 1 to verify that it is only qubit 2 that has been subjected to a Hadamard:
Sometimes we need a Hadamard on every one of the qubits; in that case you may use the command HALL:
The action of a Hadamard on state vectors is also simple to evaluate,
A rotation around the y axis by an angle θ applied to a state |0〉 can be generated by the user as:
or the user can use the QCPITT command RotY[ θ ]

b) Two Qubit operations: More on the CNOT gate

If we apply the CNOT to an initial density matrix |10〉〈10|, we write,
which trivially is seen to be equal to the density matrix for |11〉〈11|
A more general case would be to have a CNOT gate with control in qubit 3 and target in qubit 5 in a circuit with 6 qubits:

c) Two Qubit operations (CPHASE)

The CPHASE operation, used mainly in the cluster model approach, is defined by:
CPHASE | a b 〉 = (-1)^(ab) | a b 〉
Its syntax is is similar to a CNOT gate, but in this case the order of the control and target qubits is not relevant.
Therefore, CPHASE[L,q1,q2], performs a CPHASE gate on qubits q1 and q2 on a system of L qubits.
In a multi-qubit environment: of say three qubits, a CPHASE between qubits 1 and 2 is:
A different version of the same gate is defined by CPHASEV[ i, j, a ] and is used to get the CPHASE action on qubits i and j of a state vector a.
For instance, if we apply CPHASE to the two qubit vector:
CPHASE[ 2, 1, 2] [ a | 0 0 〉+ b | 0 1 〉 + c | 1 0 〉 + d | 1 1 〉 ]
= a | 0 0〉+ b | 0 1〉+ c | 1 0〉 - d | 1 1〉
Which corresponds to the following command:
Thus CPHASEV[ i , j , a] gives the result of a CPHASE acting on qubits i and j for a state vector with expansion coefficients a: {a,b,c.....}. The length of the array gives the total number of qubits.

d) Two Qubit operations, TwoOp

General Two Qubit operations are build using TwoOp. Its syntax is very simple and allows one to build any two qubit operation embedded in an L qubit system.

Circuits to Operators

Another feature to be included in a future QCPITT package is the ability to specify a circuit and then automatically generate the explicit operator or matrix corresponding to that circuit. This is done by a few additional commands. For example, we start by specifying a circuit "mat." Then use command CircuitG[matr_] to produce a cleaned-up circuit. Follow that by opB=CircuitToOp[mat, to obtain the associated operator. To look at the operator at each stage use the command CircuitToOp2[Bell][[1]] and CircuitToOp2[Bell][[2]],etc. See tutorial CircuitToOp.nb for details. A few examples are given below:
Hopefully, this will be extended to a drag and drop method for circuit assembly and also for direct input to Qiskit codes.

Conclusion

An introduction to the package QCPITT has been presented, followed by samples of its capabilities and facilities. It taps into the basic features of Mathematica and hopefully provides a viable and practical tool for learning, handling, enjoying and understanding quantum computing. In addition to the selected examples presented here there is much more available and much more is possible. Publications of earlier versions, along with tutorials and applications, are available. Some of the methods used here were adapted from another related supercomputer version called QCMPI written in Fortran. Links to these resources are provided below. Recent developments include: (1) application to qutrits and qudits; (2) extension to closed and open system dynamics using a Lindblad/Beretta formulation, (3) studies of one and two qubit gate dynamics. and (4) multi-dimensional tunnelling (in preparation). Future goals are: (1) to link to the qiskit/Jupyter setup to allow direct submission of QC jobs. (2) To provide a drag and drop method to assemble executable circuits. (3) Error correction studies (4) examine tunneling, scattering, and bound states in one to three dimensions.

Resource Links

Publications
QDENSITY—A Mathematica Quantum Computer simulation. Comput.Phys.Commun.174 June 2006,914. Bruno Juliá-Díaz, Joseph M. Burdis and Frank Tabakin. link​
​
QDENSITY-A Mathematica quantum computer simulation.Comput.Phys.Commun.180(3):474 (2009)
Bruno Juliá-Díaz ,Joseph M.Burdis,Frank Tabakin. link​
​
QCMPI:A parallel environment for quantum computing. Comput.Phys.Commun.180(6):948-964 (2009)
Frank Tabakin, Bruno Juliá-Díaz. link​
​
QCWAVE-A Mathematica quantum computer simulation update. Comput.Phys.Commun.182(8):1693-1707 (2011)
Frank Tabakin and Bruno Juliá-Díaz. link​

QDENSITY/QCWAVE:A Mathematica quantum computer simulation update. Comput.Phys.Commun.201:171-172(2016)
Frank Tabakin.
​
Model dynamics for quantum computing. Annals of Physics Volume 383, August 2017, Pages 33-78 link
Tutorials Download Tutorials

References

Michael A. Nielsen and Isaac I. Chuang,
``Quantum Computation and Quantum Information ’’, Cambridge University Press (2000).
​
Teleportation C. H. Bennett, G. Brassard, C. Crepeau, R. Jozsa, A. Peres, and W. K. Wootters,
Phys. Rev. Lett. 70, 1895-1899 (1993).
​
L. K. Grover, Phys. Rev. Lett. 79, 325-328 (1997).

Peter W. Shor, SIAM J. COMPUT 26 (5) 1484 (1997).

J. S. Bell, Rev. Mod. Phys. 38 447 (1966).

D. M. Greenberger , M. A. Horne, and A. Zeilinger, in
Bell’s Theorem, Quantum Theory, and Conceptions of the Universe,
edited by M. Kafatos (Kluwer Academics, Dordrecht, The Netherlands, 1989), pp. 73.
​
R. F. Werner,Phys. Rev. A 40, 4277 (1989).
​
G. Lindblad, Commun. Math. Phys. 4 ,119 (1976).

R. Raussendorf and H. J. Briegel, Phys. Rev. Lett. 86 (2001) 5188;
R. Raussendorf, D. Brownie and H.J. Briegel, Phys. Rev. A 68, 022312 (2003);
Michael I. Nielsen, quant - ph/0504097
​
John von Neumann, ``Mathematical Foundations of Quantum Mechanics’’,
Princeton University Press (1955).

Acknowledgements

The QCPITT project originally included Prof. Bruno Juliá-Díaz, (Universitat de Barcelona, Barcelona (Spain)), and Dr. Joseph M. Burdis, to whom I am very grateful. Questions provided by Dr. Kapil K. Sharma are also appreciated; he stimulated the extensions to hybrid systems, partial transposition, and to quantum discord. The author’s interest in QC was enhanced by communications with Dr. Victor Volkov. Earlier this project was supported by the National Science Foundation.