Solutions
to
the
Dirac
Equation
and
Their
Connection
to
Quantum
Cellular
Automata
by
Michael
Q.
May
Professor
Frederick
Strauch,
Advisor
A
thesis
submitted
in
partial
fulllment
of
the
requirements
for
the
Degree
of
Bachelor
of
Arts
with
Honors
in
Physics
WILLIAMS
COLLEGE
Williamstown,
Massachusetts
May
19,
2017
2
Abstract
We
explore
the
relationship
between
quantum
cellular
automata
(QCAs)
and
the
free
Dirac
and
Weyl
equations
in
two
and
three
dimensions.
First,
we
derive
then
characterize
localized,
positive
energy
solutions
of
the
Dirac
equation
in
two
and
three
dimensions.
We
nd
their
moments
and
conditions
on
their
spinor
components
for
symmetric
propagation.
Next
we
derive
cellular
automata
on
square
and
hexagonal
lattices,
thus
countering
D'Ariano
and
Perinotti's
New
NoGo
Lemma,
supposedly
restricting
QCAs
to
only
the
square
lattice
in
twodimensions.
Finally,
we
consider
the
recovery
of
the
Weyl
equation
from
our
numerical
simulations
of
twodimensional
automata.
We
nd
that
in
the
long
wavelength
limit,
using
the
same
spinor
conditions
found
for
the
Dirac
and
Weyl
equations
for
symmetric
propagation,
that
the
hexagonal
lattice
is
isotropic,
and
its
automaton
yields
a
dispersion
relation
identical
to
that
of
the
Weyl
equation.
i
ABSTRACT
Executive
Summary
It
is
the
goal
of
this
thesis
to
develop
methods
to
analyze
the
question
of
which
lattices
can
support
quantum
cellular
automata
and
simulate
the
Dirac
and
Weyl
equations.
0.1
Quantum
Cellular
Automata
Quantum
cellular
automata,
or
QCAs,
like
cellular
automata,
CAs,
act
on
a
lattice.
They
apply
local,
temporally
and
spatially
homogeneous
update
rules
to
the
values
of
lattice
points
at
each
timestep.
Lattice
sites
operated
on
by
QCAs,
unlike
CAs,
may
take
nondiscrete,
complex
values,
and
the
operation
of
a
QCA
must
also
be
linear
and
unitary
such
that
the
norm
of
all
lattice
sites
is
preserved.
It
has
been
proved
by
Meyer
[1]
that
any
scalar
QCA
in
one
dimension
is
identical
to
a
translation
operator
multiplied
by
a
complex
phase.
As
a
result
of
this
so
called
\NoGo
Lemma,”
any
interesting
1D
QCAs
are
further
distinguished
from
CAs
through
their
use
of
twocomponent
spinors
to
describe
each
lattice
point.
The
NoGo
Lemma
generalizes
to
higher
dimensions,
with
twoand
threedimensional
lattices
requiring
two
or
more
components
as
well.
Figure
1:
Meyer's
onedimensional
QCA
with
nonzero
mass
and
periodic
boundary
conditions.
Time
runs
from
left
to
right.
More
recently,
D'Ariano
and
Perinotti
[2]
have
proposed
stringent
conditions
on
the
lattices
on
which
QCAs
may
operate.
Their
\New”
NoGo
Lemma
states
that
there
do
not
exist
temporally
and
spatially
homogeneous
QCAs
on
lattices
in
two
and
three
dimensions
besides
the
square
and
bodycentered
cubic,
respectively.
This
New
NoGo
Lemma
has
played
an
important
part
in
directing
their
continuing
research
as
well
as
the
research
of
others
in
the
area
of
QCAs.
Though
there
continues
to
be
work
on
nonsquare
cellular
automata
[3],
much
iii
EXECUTIVE
SUMMARY
(a)
(b)
Figure
2:
Translation
steps
on
a
hexagonal
lattice
in
the
multistep
(a)
and
reduced
(b)
schemes.
recent
research
on
twodimensional
quantum
cellular
automata
has
focused
exclusively
on
square
automata
[4],
even
in
the
related
area
of
quantum
walks
[5,
6,
7].
0.1.1
Countering
the
New
NoGo
Lemma
Following
the
work
of
BialynickiBirula
[8],
we
can
construct
a
hexagonal
cellular
automaton
using
a
temporally
inhomogeneous
(multistep)
process.
We
unitarily
update
along
one
of
the
three
hexagonal
translation
axes
at
each
time
step,
with
the
full
action
of
the
automaton
shown
in
Fig.
2a.
After
three
time
steps,
the
initial
lattice
point
has
propagated
to
seven
points,
to
six
of
its
nextnearest
neighbors
and
itself.
Because
each
of
the
points
of
propagation
lies
along
an
original
translation
axis,
we
are
motivated
to
rewrite
the
action
of
the
QCA
as
a
singlestep
operation,
shown
in
Fig.
2b.
This
now
temporally
homogeneous
QCA
moves
to
each
of
its
nearest
neighbors
at
each
timestep,
with
some
probability
of
not
propagating
at
all.
Our
construction
of
this
QCA
acts
as
a
counterexample
to
D'Ariano
and
Perinotti's
New
NoGo
Lemma,
opening
up
a
relatively
unexplored
class
of
QCAs
for
consideration.
0.2
Connecting
QCAs
to
the
Dirac
and
Weyl
Equation
s
Having
shown
that
in
addition
to
the
square
automaton,
there
exists
a
hexagonal
QCA,
we
may
ask
how
numerical
simulations
of
these
automata
compare
with
solutions
to
the
Dirac
and
Weyl
equations.
To
that
end,
we
derive
localized,
positive
energy
solutions
to
the
Dirac
equation
in
two
and
three
dimensions
to
determine
conditions
on
their
initial
spinor
0.2.
CONNECTING
QCAS
TO
THE
DIRAC
AND
WEYL
EQUATIONS
v
3
(a)
Ca
=
0
(b)
Ca
=
(c)
Ca
=1
2
Figure
3:
Density
plot
of
3D
wavepacket
projected
onto
the
xy
plane
at
t
=
1,
localization
parameter
a
=
:2,
for
various
initial
spinor
values.
Figure
4:
Hexagonal
QCA
after
50
time
steps
with
an
initial
rudimentary
wavepacket.
Each
of
the
seven
initial
nonzero
cells
has
only
a
single
nonzero
spinor
component.
components
for
symmetric
propagation.
As
can
be
seen
in
Fig.
3,
when
the
spinor
parameter
Ca
=
1,
there
will
be
symmetric
propagation.
It
happens
that
a
solution
to
the
Dirac
equation
and
a
QCA
require
the
same
number
of
spinor
components
per
dimension
as
one
another.
Thus,
allowing
the
analogous
Ca
parameter
in
our
hexagonal
automaton
to
be
1
also
results
in
symmetric
propagation.
We
also
nd
that
in
accordance
with
our
Dirac
equation
solutions,
a
wavepacket
on
the
hexagonal
lattice
spreads
with
a
welldened
wavefront
at
nearly
the
speed
of
light
on
the
lattice.
Indeed,
in
the
long
wavelength
limit
our
hexagonal
QCA,
show
in
Fig.
4,
and
our
Weyl
equation
solution
give
identical
dispersion
relations.
Future
work
will
apply
the
methods
used
to
derive
the
hexagonal
automaton
to
the
derivation
of
other
QCAs
on
nonsquare
lattices.
Additionally,
the
relationship
between
the
hexagonal
QCA
and
the
Dirac
equation
will
be
more
fully
explored
by
introducing
inertial
or
chiral
mass
terms
to
the
automaton.
EXECUTIVE
SUMMARY
Acknowledgments
First,
I
would
like
to
thank
my
advisor,
Professor
Frederick
Strauch,
for
all
of
his
guidance
and
mentorship
throughout
the
year.
I
began
work
with
him
uncertainly,
having
previously
sworn
off
Bessel
functions
after
his
Math
Methods
course
my
freshman
year.
His
passion
for
physics,
however,
and
his
kind
encouragement
proved
to
change
my
mind
on
the
matter.
This
thesis
would
not
have
been
possible
without
the
direction
he
provided.
Second,
I
want
to
thank
Professor
TuckerSmith
for
his
helpful
feedback
as
the
second
reader.
Finally,
I
have
appreciated
the
ongoing
support
of
my
friends
and
family
throughout
my
college
career.
Thank
you.
vii
ACKNOWLEDGMENTS
Contents
Abstract
i
Executive
Summary
iii
0.1
QuantumCellularAutomata...........................
iii
0.1.1
CounteringtheNewNoGoLemma
...................
iv
0.2
ConnectingQCAstotheDiracandWeylEquations
.
.
.
.
.
.
.
.
.
.
.
.
.
.
iv
Acknowledgments
vii
1
Introduction
1
1.1
DerivationoftheDiracEquation
........................
1
1.2
IntroductiontoCellularAutomata
.......................
3
1.2.1
The
OneDimensional
Quantum
Cellular
Automaton
.
.
.
.
.
.
.
.
.
4
1.2.2
HigherDimensional
Quantum
Cellular
Automata
.
.
.
.
.
.
.
.
.
.
.
5
1.3
OutlineforPaper
.................................
7
2
Solving
the
Dirac
Equation
9
2.1
TheGeneralSolution
...............................
9
2.1.1
Solving
for
I
(r;t)
............................
10
2.2
WavepacketswithSpinorComponents
.....................
13
2.2.1
TwoDimensions
.............................
13
2.2.2
ThreeDimensions.............................
14
3
Moments
of
Dirac
Equation
Solutions
19
3.1
2DMoments....................................
19
3.1.1
FirstMoments
..............................
19
3.1.2
SecondMoments
.............................
20
3.2
3DMoments....................................
21
3.2.1
FirstMoments
..............................
21
3.2.2
SecondMoments
.............................
22
ix
CONTENTS
4
Quantum
Cellular
Automata
25
4.1
SquareCellularAutomata
............................
25
4.1.1
Numerical
Simulation
of
Square
Automaton
.
.
.
.
.
.
.
.
.
.
.
.
.
.
26
4.2
HexagonalQuantumCellularAutomata
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
28
4.2.1
NewNoGoLemma............................
30
5
Discussion
33
5.1
TheNewNoGoLemma
.............................
33
5.2
FutureWork....................................
33
5.2.1
ExtendingtheDiracEquationSolution
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
33
5.2.2
AddingMass
...............................
34
5.2.3
OtherLattices...............................
34
5.2.4
Analytic
Connections
of
QCAs
to
the
Dirac
Equation
.
.
.
.
.
.
.
.
.
35
A
Appendix
37
A.1
MomentCalculations
...............................
37
A.1.1
TwoDimensions
.............................
37
A.1.2
ThreeDimensions.............................
38
A.2
Hexagonal
Automaton
Unitarity
Conditions
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
38
List
of
Figures
1
OneDimensionalQCAEvolution
........................
iii
2
HexagonalQCATranslationSteps........................
iv
3
ThreeDimensional
Wave
Packet
Density
Plots
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
v
4
HexagonalQCAEvolutionofaWavepacket
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
v
1.1
1DCellularAutomata,Rule110.........................
3
1.2
OneDimensionalQCAEvolution
........................
4
1.3
TwoDimensionalLattices
............................
7
2.1
I
(r,
t)atFourTimes
..............................
13
2.2
TwoDimensionalWavePackets
.........................
14
2.3
TwoDimensional
Wave
Packet
Localization
Plots
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
15
2.4
ThreeDimensional
Wave
Packet
Density
Plots
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
17
3.1
Time
Independent
Elements
of
3D
Second
Moments
.
.
.
.
.
.
.
.
.
.
.
.
.
.
23
3.2
TimeDependentElementsof3DSecondMoments.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
23
4.1
TranslationStepsforSquareQCA........................
26
4.2
ReducedTranslationStepsforSquareQCA
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
27
4.3
SquareQCAEvolutionfromaSinglePoint.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
27
4.4
HexagonalTranslationSteps
...........................
29
4.5
HexagonalReducedTranslationSteps......................
30
4.6
HexagonalQCAEvolutionfromaSinglePoint.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
31
4.7
HexagonalQCAEvolutionofaWavepacket
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
32
5.1
GrapheneLattice
.................................
34
xi
LIST
OF
FIGURES
Chapter
1
Introduction
Quantum
cellular
automata
(QCAs)
are
intimately
tied
to
fundamental
equations
in
physics.
At
the
least,
this
means
that
areas
such
as
quantum
information
theory
can
benet
from
the
insights
gained
by
studying
the
discretization
of
normally
continuous
unitary
operators
in
QCAs
[1].
The
Weyl
and
Dirac
equations,
governing
relativistic
wave
and
particle
behavior,
could
potentially
be
fully
simulated
using
QCAs
[8].
On
a
grander
scale,
QCAs
may
be
used
to
probe
questions
of
discrete
space
and
time
[9].
In
this
paper,
we
explore
the
relationship
between
QCAs
and
the
Dirac
equation
in
two
and
three
dimensions,
providing
new
insights
into
the
various
lattices
on
which
a
QCA
yielding
the
Dirac
propagator
may
act.
We
begin
with
a
derivation
of
the
Dirac
equation
and
explain
its
meaning
compared
to
the
Schrodinger
equation.
We
then
introduce
cellular
automata
and
quatum
cellular
automata
before
exploring
their
connection
to
the
Dirac
equation.
The
New
NoGo
Lemma,
which
places
strong
restrictions
on
which
lattices
QCAs
may
operate,
is
introduced
and
explained.
1.1
Derivation
of
the
Dirac
Equation
A
generalized
wave
equation
may
be
written,
^
E =
H^
,
(1.1)
^
where
we
have
E
=
i@t
and
the
Hamiltonian
is
a
function
of
spatial
derivatives.
The
Schrodinger
equation
describes
the
time
evolution
of
nonrelativistic
probability
distributions
2
where
H^
=
2p^
m
and
^p
=
..i'
(we
have
set
h
=
1).
It
has
rst
order
time
derivatives
and
second
order
spatial
derivatives
which
mirrors
the
classical
energymomentum
relation.
2
p
E
=
(1.2)
2m
If
one
wants
to
take
special
relativity
into
account,
however,
energy
and
momentum
should
be
of
the
same
order,
i.e.
the
relativistic
energymomentum
equation
should
be
used.
22
24
E2
=
cp
+
mc
(1.3)
1
CHAPTER
1.
INTRODUCTION
Simply
making
the
substitution
for
the
relativistic
energy
in
Eq.
1.1
does
not
give
a
satisfactory
solution
to
the
problem,
though.
It
cannot
take
into
account
particle
spin,
and
it
lacks
rstorder
time
derivatives
just
as
the
Schrodinger
equation
does.
In
1928,
Paul
Dirac
realized
that
by
allowing
the
wavefunction
to
be
multicomponent,
01
1
BC
2
BC
=
B.
C,
(1.4)
.
@.
A
n
in
addition
to
using
the
relativistic
energymomentum
equation,
he
could
construct
a
relativistic
wave
equation
which
was
rst
order
in
both
time
and
space.
Letting
H^
pc+mc2
=
^,
where
a
and
ß
are
matrices,
and
requiring
Eq.
1.3
be
satised,
one
nds
that
in
one
dimension
the
wavefunction
must
have
two
components.
Additionally,
the
matrices
a
and
ß
must
satisfy
2
=
I,
(1.5)
2
=
I,
(1.6)
ß
+
a
=0.
(1.7)
We
may
write
the
a
or
ß
matrices
in
terms
of
the
the
Pauli
matrices.
01
0
..i
10
x
=;y
=;z
=(1.8)
10i
00
..1
Our
full
relativistic
wave
equation
in
one
dimension1
can
thus
be
written,
10
0
..i
12 1
i@t=pc
+mc
.
(1.9)
20
..1i
0 2
In
two
dimensions,
the
a
matrix
becomes
two
a
matrices,
one
corresponding
to
each
momentum
component.
In
this
case
we
may
again
use
the
Pauli
matrices,
giving
.
two
components.
In
general,
for
n
dimensions,
one
will
need
na
matrices
which
satisfy
jß
+
j
=
0
(1.10)
jk
+
kj
=2jkI.
(1.11)
Since
there
do
not
exist
any
set
of
threebythree
matrices
satisfying
these
conditions,
threedimensional
solutions
to
the
Dirac
equation
must
have
four
by
four
a
matrices.
These
may
be
written
in
a
number
of
ways,
such
as
the
tensor
products
of
the
Pauli
matrices.
We
will
use
a
variation
on
the
Weyl
basis,
such
that
0
=
1
.
I2;
1
=
3
.
1;
2
=
3
.
2;
3
=
3
.
3.
(1.12)
1This
derivation
of
the
Dirac
equation
follows
closely
that
found
in
Thaller
[10]
1.2.
INTRODUCTION
TO
CELLULAR
AUTOMATA
Setting
s
=(x;y;z),
we
can
write
these
Dirac
matrices
as
.
.
a
=
.
=
s
0
0
..s
,
(1.13)
where
ß
=
0.
For
the
massless
case,
the
Dirac
equation
reduces
to
the
Weyl
equation,
and
we
can
again
use
Pauli
matrices.
1.2
Introduction
to
Cellular
Automata
Cellular
Automata
(CAs)
act
on
a
discretevalued
lattice
in
discrete
time
steps
with
a
local,
spatially
homogeneous
update
rule.
That
is,
every
site
can
only
\see”
the
sites
around
it,
and
each
site
must
update
using
the
same
rule.
The
simplest
cases
of
CAs
are
those
rst
exhaustively
described
by
Wolfram
wherein
at
each
time
step
a
lattice
point
considers
the
values
of
itself
and
its
two
nearest
neighbors
to
decide
its
new
binary
value
[11].
There
are
256
possible
1D
automata,
although
very
few
are
interesting.
Some,
however,
show
surprising
complexity
for
their
simple
rules.
Shown
in
Fig.
1.1
is
rule
110
which,
given
the
appropriate
initial
condition,
can
simulate
any
calculation.
A
Quantum
Cellular
Automaton
(QCA)
acts
similarly
to
a
CA,
except
lattice
sites
may
take
nondiscrete,
complex
values
and
the
operation
of
the
QCA
must
further
be
linear
and
unitary
such
that
the
norm
of
all
lattice
sites
is
preserved.
Further,
for
a
QCA
to
yeild
the
Weyl
or
Dirac
propagators,
it
must
be
symmetric
on
a
1D
lattice
and,
(a)Sixiterations.
more
generally,
isotropic
on
higherdimensional
lattices.
It
has
been
proved
by
Meyer
[1]
that
any
scalar
QCA
in
one
dimension
is
identical
to
a
translation
operator
multiplied
by
a
complex
phase.
As
a
result
of
this
so
called
\NoGo
Lemma,”
any
interesting
1D
QCAs
are
further
distinguished
from
CAs
through
their
use
of
two
component
spinors
to
describe
each
lattice
point.
This
NoGo
Lemma
generalizes
to
higher
dimensions,
with
twoand
threedimensional
lattices
requiring
two
or
more
components.
The
necessity
for
multicomponent
spinors
is
analogous
to
that
of
the
Dirac
Equation,
with
each
dimension
requiring
the
same
order
of
spinor.
(b)
One
hundred
iterations.
Figure
1.1:
Onedimensional
cellular
automata
updating
using
Wolfram's
rule
110.
Time
proceeds
from
top
to
bottom.
CHAPTER
1.
INTRODUCTION
1.2.1
The
OneDimensional
Quantum
Cellular
Automaton
The
simplest
1D
QCAs
have
twocomponent
spinors,
whose
components
can
be
thought
of
as
leftmoving
and
rightmoving.
With
each
time
step,
there
is
some
probability
that
each
component
will
move
in
the
appropriate
direction
and
some
probability
that
each
component
will
ip
directions,
accruing
a
complex
phase
in
the
process.
We
can
write
this
QCA
in
its
most
general
form
as
R(n,
t
+
1)
=
cos()
(cos() R(n

1;t)+
i
sin() L(n
+1;t))
+
sin()(..i
cos() R(n,
t)
+
sin() L(n,
t))
L(n,
t
+
1)
=
cos()
(cos() L(n

1;t)+
i
sin() R(n
+1;t))
+
sin()
(sin() R(n,
t)

i
cos() L(n,
t))
,
(1.14)
where
R
and
L
are
the
right
moving
and
left
moving
components
of
.
This
cellular
automaton
has
been
been
exhaustively
characterized
by
Meyer
[1,
12]
in
some
of
the
seminal
papers
on
QCAs.
Of
note
in
this
setup
are
the
two
elements
.
and
,
which
both
might
be
called
mass
elements.
The
sin()
terms
cause
(n,
t
+
1)
to
depend
on
itself
in
addition
to
its
neighbors;
we
will
call
such
terms
the
inertial
mass.
The
sin()
terms
cause
the
left
and
right
components
of
.
to
switch
back
and
forth;
such
terms
are
analogous
to
the
mass
terms
in
the
Dirac
equation.
Using
this
analogy,
we
will
call
.
the
chiral
mass.
Without
an
inertial
mass
term,
i.e.
setting
.
=
0,
this
QCA
acts
on
two
noninteracting
lattices
which
alternate
positions
in
time.
Adding
a
mass
term
couples
the
two
QCAs.
In
terms
of
what
that
would
physically
look
like
in
the
context
of
Fig.
1.2,
the
propagation
cone
would
separate
from
the
light
cone.
Instead
of
propagating
at
an
angle
of
p
2
,
the
wavefunction
would
spread
at
a
shallower
angle.
Figure
1.2:
Meyer's
onedimensional
QCA
with
.
=
p
2
,
nonzero
mass,
and
periodic
boundary
conditions.
Time
runs
from
left
to
right.
Another
way
to
write
the
action
of
Eq.
1.14
is
to
split
its
operations
into
a
coin
component
and
a
translation
component
acting
on
two
dierent
Hilbert
spaces.
We
can
thus
rewrite
the
action
as
U
=
A0
.
I+
A+
.
T
+
A
.
T
..1
,
(1.15)
where
the
rst
components,
An,
are
2
×
2
matrices
acting
on
.
and
the
second
components,
1.2.
INTRODUCTION
TO
CELLULAR
AUTOMATA
T
1,
are
translation
operators
acting
on
the
lattice.
In
the
case
above,
A+
=
cos()
0
0
i
sin()
cos()
,
(1.16)
A
=
cos()
cos()
i
sin()
0
0,
(1.17)
A0
=
sin()
sin()
..i
cos()
..i
cos()
sin()
.
(1.18)
2
It
is
easy
to
check
that
U
is
in
fact
unitary.
I
=
UyU
=
A†
A0
+
A†
A+
+
Ay..A
(1.19)
0+
Ay
+..A0
+
A†
0A+T
+
A+†
A0
+
A†
0A..T
..1
2
+Ay..A+T
2
+
A†
+A..T
..2
(1.20)
We
only
need
to
ensure
that
Eq.
1.19
sums
to
the
identity
and
that
each
of
the
elements
in
1.20
are
zero.
Previous
work
on
this
1D
QCA
has
shown
it
to
be
linked
to
the
Dirac
Equation.
Indeed,
in
the
long
wavelength
limit,
it
has
been
proved
that
the
1D
QCA
and
the
Dirac
equation
give
equivalent
propagators
[1].
1.2.2
HigherDimensional
Quantum
Cellular
Automata
In
dimensions
higher
than
one,
the
general
construction
of
the
QCA
is
not
at
all
obvious
even
when
the
mass
is
zero.
In
two
dimensions
one
can
choose
between
the
square
lattice,
a
natural
extension
of
the
onedimensional
lattice,
where
each
point
has
four
neighbors,
or
more
exotic
graphs
such
as
the
hexagonal,
graphene,
or
kagome
graphs
with
each
point
having
six,
three,
and
four
neighbors,
respectively.
One
of
the
rst
3D
QCAs
was
introduced
by
BialynickiBirula
[8].
Starting
from
the
condition
of
unitarity,
he
derives
the
most
general
QCA
on
a
bodycentered
cubic
lattice
with
8
neighbors.
He
does
not
go
into
detail
about
the
the
primitive
cubic
and
the
facecentered
cubic
lattices
with
6
and
12
neighbors,
though.
Having
chosen
primitive
translations,
further
complications
arise.
For
instance,
the
exact
order
of
those
translations
can
lead
to
translation
steps
on
dierent
eective
lattices.
BialynickiBirula
[8]
and
Strauch
[13],
using
the
same
primitive
translations,
order
them
such
that
their
QCAs
operate
on
a
bodycentered
cubic
lattice
and
something
more
complicated,
respectively.
Using
a
Fourier
representation
for
the
translation
operators,
we
can
write
the
unitary
operators
..ikxx
e..ikyy
e..ikzz
BialynickiBirula
=
e,
(1.21)
ky
ky
..i
kx
x
..iy
..ikzz
e..iy
..i
kx
x
Strauch
=
(1.22)
e
e
2
e
2
e
.
CHAPTER
1.
INTRODUCTION
The
eigenvalues
of
these
matrices,
ei!,
correspond
to
the
dispersion
relations:
cos(!BB)
=
cos(kx)
cos(ky)
cos(kz)
+
sin(kx)
sin(ky)
sin(kz),
(1.23)
cos(!S)
=
cos(kx)
cos(ky)
cos(kz),
(1.24)
which
can
be
expanded
in
the
long
wavelength
limit,
where
kx,
ky,
kz
«
1.
Keeping
only
terms
up
to
the
third
order,
we
have,
!2
BB
=
kx
2
+
ky
2
+
kz
2

2kxkykz,
(1.25)
!S
2
=
kx
2
+
ky
2
+
kz
2
.
(1.26)
Eq.
1.26
is
identical
to
the
dispersion
relation
of
the
Weyl
equation,
but
Eq.
1.25
includes
a
higherorder
term.
This
shows
that
dierent
QCAs
can
better
simulate
fundamental
equations
than
others,
and
there
are
multiple
possibilities
for
valid
update
rules
in
higher
dimensions.
New
NoGo
Lemma
Despite
the
seeming
breadth
of
possibilities
of
QCAs
in
dimensions
greater
than
one,
there
are
in
fact
many
restrictions
which
can
be
placed
on
them.
In
two
dimensions,
for
example,
there
are
only
two
topologically
distinct
Bravais
lattices
on
which
a
cellular
automata
may
act,
namely
the
hexagonal
and
square
lattices.
These
lattices
may
be
constructed,
assuming
a
lattice
spacing
of
1,
as:
rsquare
=
n1x^+
n2y^(1.27)
v
rhexagonal
=
n1..1
2
x^+
2
3
y^+
n2x^(1.28)
8n1;n2
.
Z.
All
other
twodimensional
Bravais
lattices
simply
rotate
the
square
lattice,
lengthen
the
step
size
in
the
x
or
y
directions
(i.e.,
changing
the
numerical
factor
in
front
of
n1x^or
n2y^to
a
nonunity
value),
or
skew
the
the
x
and
y
axes.
Other
regular
graphs,
such
as
that
of
graphene,
do
not
satisfy
the
condition
of
spatial
homogeneity
(see
Fig.
1.3).
Specically,
despite
there
being
six
dierent
translation
vectors
in
graphene,
as
in
the
hexagonal
lattice,
not
all
translation
vectors
are
available
at
each
position
on
the
graph.
Indeed,
at
any
point,
only
three
translation
vectors
may
be
dened.
While
arguments
may
be
made
that
there
exists
a
clear
global
sense
of
homogeneity,
since
the
translation
rules
for
each
point
cannot
be
mapped
on
every
other
point
on
the
graph,
graphene
does
not
allow
for
a
primitive
QCA.
In
this
same
spirit
of
spatial
homogeneity,
one
could
further
require
QCAs
to
satisfy
strict
temporal
homogeneity.
That
is
to
say
that
in
addition
to
every
point
on
a
lattice
being
subject
to
the
same
translation
rules,
at
every
discrete
time
step
the
QCA
must
implement
the
same
action.
This
would
rule
out
a
multistep
rule
for
graphene
in
which
the
two
sublattices
updated
on
alternating
timesteps.
1.3.
OUTLINE
FOR
PAPER
(c)
Graphene
Lattice
(a)
Square
Lattice
(b)
Hexagonal
Lattice
Figure
1.3:
Twodimensional
lattices.
Note
that
the
red
and
green
sites
in
graphene
are
geometrically
distinct,
while
all
sites
on
the
square
and
hexagonal
lattices
are
identical.
In
light
of
the
above,
one
might
ask
how
strongly
these
conditions
of
spatial
and
temporal
homogeneity
restrict
the
possible
lattices
for
QCAs
in
two
and
three
dimensions.
Bialynicki
Birula
[8]
hinted
that
the
bodycentered
cubic
lattice
may
be
the
only
one
to
support
a
3D
QCA.
More
recently,
D'Ariano
and
Perinotti
[2]
have
apparently
proven
that
the
temporal
homogeneity
condition
restricts
QCAs
to
only
the
square
lattice
in
two
dimensions
and
the
bodycentered
cubic
lattice
in
three.
This
\New”
NoGo
Lemma
has
played
an
important
part
in
directing
their
continuing
research
as
well
as
the
research
of
others
in
the
area
of
QCAs.
Though
there
has
continued
to
be
work
on
nonsquare
cellular
automata
[3],
much
recent
research
in
twodimensional
cellular
automata
has
focused
exclusively
on
square
automata
[4],
even
in
the
related
area
of
quantum
walks
[5,
6,
7].
It
is
the
goal
of
this
thesis
to
develop
methods
and
analyze
the
question
of
what
lattices
can
support
QCAs
and
simulate
the
Dirac
and
Weyl
equations.
1.3
Outline
for
Paper
To
help
determine
a
connection
between
QCAs
and
the
Dirac
and
Weyl
equations,
in
Chapter
2
we
nd
closedform
solutions
to
the
Dirac
equation
in
two
and
three
dimensions.
These
positive
energy,
localized,
wavepackets
are
analogous
to
that
of
a
spreading
wavepacket
on
a
QCA.
Such
wavepacket
solutions
are
very
useful
to
interpret
the
evolution
of
the
QCA.
We
determine
the
necessary
conditions
on
their
initial
spinor
components
for
symmetric
propagation,
and
conrm
the
use
of
the
localization
parameter
in
our
solution.
We
can
use
these
solutions
to
compare
with
spreading
wavepackets
constructed
on
QCAs,
especially
those
which
have
not
yet
been
analytically
shown
to
limit
to
the
Dirac
equation.
To
more
fully
describe
the
wavepackets
found
in
Chapter
2,
moments
of
the
wavepackets
are
calculated
in
Chapter
3.
We
conrm
the
properties
of
our
packets
match
with
the
properties
we
should
expect
from
a
highly
localized,
relativistic
wavepacket
and
determine
conditions
for
symmetric
propagation
of
the
packet.
We
also
check
that
the
wavepackets
propagate
causally.
CHAPTER
1.
INTRODUCTION
In
Chapter
4
we
present
numerical
simulations
of
2D
QCAs
on
square
and
hexagonal
lattices
and
discuss
their
signicance
in
the
long
wavelength
limit
when
compared
to
our
solutions
of
the
Weyl
equation.
We
also
determine
limits
of
the
New
NoGo
Lemma,
showing
that
in
certain
circumstances
the
hexagonal
QCA
can
support
a
unitary
operation
subject
to
the
constraint
of
temporal
homogeneity.
In
Chapter
5,
we
discuss
our
conclusions
with
regard
to
the
New
NoGo
Lemma
and
our
derived
hexagonal
QCA.
We
discuss
extensions
and
future
work
in
connecting
higherdimensional
QCAs
and
the
Dirac
equation.
Chapter
2
Solving
the
Dirac
Equation
To
guide
our
work,
we
calculate
localized,
positive
energy
solutions
to
the
Dirac
Equation
in
two
and
three
dimensions
following
the
method
prescribed
in
Strauch
[14].
This
method
reduces
the
problem
of
nding
a
wavepacket
solution
to
a
fundamental
function
and
an
initial
spinor.
In
general,
the
wavepacket
spreads
asymmetrically,
with
a
direction
determined
by
initial
spinor
components.
Conditions
on
symmetric
propagation
of
the
solutions
are
found
and
compared
with
the
onedimensional
solution.
2.1
The
General
Solution
It
has
been
shown
in
Strauch
[14]
that
given
an
equation
of
the
form
of
Eq.
1.1,
one
can
nd
a
large
set
of
solutions
of
the
form
Y
j(r,
t)=!k(..ir)

H^
(..ir)I
(j)(r;t),
(2.1)
k6
=j
where
I
(j)(r;t)
satises
the
equation,
i@tI
(j)(r;t)=
!j(..ir)I
(j)(r;t),
(2.2)
and
where
!j(..ir)
are
the
partial
eigenvalues
of
the
diagonalized
nnH^
matrix
and
(..ir)
is
an
arbitrary
dierential
operator
with
n
components.
We
may
restrict
our
solutions
to
those
which
satisfy
the
canonical
dispersion
relation
with
h
=
c
=
1,
X
!22
=
m
+ki
2
=
m
2
..r2
,
(2.3)
i
and
further
let
f
be
simply
an
arbitray
ndimensional
spinor.
Now,
we
may
write
positive
energy
solutions
to
the
Dirac
equation
in
the
form,
(r,
t)=
N..i@t

H^
I
(r;t),
(2.4)
2
9
CHAPTER
2.
SOLVING
THE
DIRAC
EQUATION
where
N
is
a
normalization
constant,
the
second
term
corresponds
to
taking
the
positive
energy
projection
of
the
wave
equation,
and
f
is
an
arbitrary
multidimensional
spinor.Recall
that
f
must
have
two
components
in
one
and
two
dimensions,
and
four
components
in
three
dimensions.
A
set
of
solutions
to
Eq.
2.2
can
be
written
as
Z
..!(k)
I
(r;t)=
1
ee
i(kr..!(k)t)dkv,
(2.5)
v
g(k)
where
g(k)
is
an
arbitrary
function
of
k.
We
will
set
g(k)
.
!(k)
for
the
rest
of
the
paper,
following
the
onedimensional
derivation
of
I
(r;t)
described
in
Strauch
[14].
With
analogy
to
the
localization
parameter
of
the
Schrodinger
Gaussian
wavepacket,
a
localization
parameter
a
has
been
introduced.
2.1.1
Solving
for
I
(r;t)
In
one
dimension
I
(r;t)
may
be
identied
with
the
zerothorder
modied
Bessel
function
of
second
kind,
whereas
in
two
dimensions
the
integral
must
be
solved
numerically.
The
solution
for
I
(r;t)
in
three
dimensions,
similar
to
that
in
one
dimension,
is
also
a
modied
Bessel
function
of
the
second
kind,
but
the
rstorder
rather
than
zeroth.
Two
Dimensions
In
two
dimensions
we
use
the
canonical
dispersion
relation
!2
=
kx
2
+
ky
2
+
m2
.
For
simplicity,
we
choose
to
evaluate
the
integral
along
along
xaxis
to
write
k
·
r
=(kx;ky)
·
(x,
y)=
kxr
(2.6)
Rewriting
I
(r;t)
in
polar
coordinates
thus
yields,
Z
..!(k)
1
ee
i(kxr..!(k)t)dkxdky
(2.7)
v
!(k)
Z2p
Z8
ke..(+it)!(k)
=
1
e
ikr
cos()dkd.
(2.8)
00
!(k)
The
angular
part
of
the
integrand
may
be
identied
with
the
integral
representation
of
the
Bessel
function
of
the
rst
kind
[15],
Z21
i(nt..x
sin(t))dt;
Jn(x)=
e
(2.9)
2
0
where
n
=
0.
Our
nal
solution
is
thus,
Z8
I
(r,
t)=2
1
ke..(+it)!(k)J0(kr)dk.
(2.10)
0
!(k)
2.1.
THE
GENERAL
SOLUTION
Eq.
2.10
does
not
have
a
closed
form
solution
and,
from
this
point,
must
be
evaluated
numerically.
We
can
simplify
the
expression,
though,
in
the
massless
case.
Then,
.
=
v
m2
+
k2
=
k,
and
our
expression
for
I
(r,
t)
simplies
to
Z8
I
(r,
t)=2e
..(+it)kJ0(kr)dk,
(2.11)
0
which
must
also
be
evaluated
numerically.
Three
Dimensions
An
analytical
solution
to
Eqn.
2.5
may
be
found
in
three
dimensions
using
a
similar
setup
as
in
two
dimensions.
Using
again
the
expected
dispersion
relation,
!2
=
k2
+
k2
+
k2
+
m2
,
xyz
as
well
as
evaluating
the
integral
along
the
zaxis,
we
nd
Z
..!(k)
I
(r;t)=
1
ee
i(kzr..!(k)t)dkxdkydkz.
(2.12)
v
!(k)
The
rst
step
in
solving
this
integral
will
be
to
convert
to
spherical
coordinates.
Z8
Zp
Z2p
..(+it)!(k)
I
(r,
t)=
1
ee
ikr
cos()k2
sin()
df
d.
dk.
(2.13)
00!(k)
0
Next,
we
make
the
substitution
q
=
cos()
to
nd,
Z8
Z1
..(+it)!(k)k2
I
(r,
t)=2
1
eeikrqdq
dk
Z8
0
!(k)
..1
11
..
..(+it)!(k)k2
ikr

e
..ikr
=2e
edk
(2.14)
0
!(k)
ikr
After
identifying
the
exponentials
as
a
sine
function,
we
may
then
rewrite
the
sine
in
terms
of
a
dervative
of
a
cosine,
at
the
same
time
making
the
substitution,
+
=
a
+
it.
Because
the
.
and
cosine
are
even
functions,
the
integral
may
also
be
extended
to
..8
to
1.
Z8
I
(r,
t)=
..4p
11
e
..(+it)!(k)k
sin(kr)dk
r0
!(k)
Z1
11
d
..+!(k)
=
..2e
cos(kr)dk.
(2.15)
r!(k)
dr
..8
We
can
remove
the
derivative
with
respect
to
r
from
the
integral
since
only
the
cosine
function
has
r
dependence
and
rewrite
the
cosine
as
an
exponential
function,
remembering
that
we
are
taking
the
integral
over
symmetric
bounds.
We
can
then
rewrite
as,
CHAPTER
2.
SOLVING
THE
DIRAC
EQUATION
Z1
1
d
1
..+!(k)
I
(r,
t)=
..2p
ee
ikrdk.
(2.16)
rdr!(k)
..8
Next
we
make
the
substitution
k2
=
m2
sinh2(y),
such
that
!2
=
k2
+
m
2
=
m
2(sinh2(y)+1)
=
m
2
cosh2(y)
(2.17)
and
thus
Z1
1
d
imr
sinh(y)..m+
cosh(y)dy.
I
(r,
t)=
..2e
(2.18)
r
dr
..8
In
order
to
reduce
the
components
in
the
exponential
to
a
single
term,
we
use
the
hyperbolic
trigonometric
identity,
cosh(x
+
y)
=
cosh(x)
cosh(y)
+
sinh(x)
sinh(y).
(2.19)
and
make
the
following
substitutions:
p
s
=r2
+
2
+
z
=
..is
cosh()
r
=
..is
sinh().
(2.20)
I
(r,
t)
may
now
be
written
in
the
form,
1
d
Z1..ms
cosh(y..)dy.
I
(r,
t)=
..2e
(2.21)
r
dr
..8
0
Finally,
we
can
make
the
substitution
y=
y

f
to
yield,
1
d
Z1..ms
cosh(y)dy.
I
(r,
t)=
..2e
(2.22)
r
dr
..8
This
is
of
the
same
form
as
the
integral
representation
of
the
modied
Bessel
function
of
the
second
kind
[15],
Z8
Kn(z)=
1
e
..z
cosh(x)
cosh(nx)dx,
(2.23)
2
..8
where
n
=
0.
Substituting
in
for
K0
we
now
have
1
d
I
(r,
t)=
..4K0(ms).
(2.24)
r
dr
Taking
the
derivative
using
the
identity,
dK0(z)
=
..K1(z),
(2.25)
dz
2.2.
WAVEPACKETS
WITH
SPINOR
COMPONENTS
Figure
2.1:
I
(r;t)
along
the
xaxis
with
scaling
parameter
a
=
:5
for
ve
times,
t
=
0,
2:5,
5,
7:5,
and
10
at
two
dierent
ranges.
gives
our
nal
answer
for
I
(r,
t)
in
three
dimensions:
I
(r,
t)=
4p
1
s
K1(ms)
(2.26)
p
=4p
v
1
K1(mx2
+
y2
+
z2
+(a
+
it)2).
(2.27)
x2+y2+z2+(+it)2
As
can
be
seen
in
Fig.
2.1,
I
(r,
t)
spreads
with
time
almost
identically
to
the
1D
Dirac
solution/quantum
walk
discussed
in
Chapter
1.
At
time
t
=
0,
the
wavepacket
resembles
a
narrow
Gaussian.
For
t>
0
the
function
develops
a
welldened
wavefront
which
propogates
at
approximately
the
speed
of
light,
one
unit
length
per
time.
2.2
Wavepackets
with
Spinor
Components
2.2.1
Two
Dimensions
In
two
dimensions
Eq.
2.4
may
be
written
with
f
as
an
arbitrary
two
component
vector
and
H
dened
as:
H^
=
1px
+
2py
+
m
=
xpx
+
ypy
+
zm
m
..i@x

@y
=.
(2.28)
@y

i@x
..m
After
substituting
^
H
and
our
solution
for
I
(r,
t),
we
can
write
(r,
t):
..i@t

m@y
+
i@x
C0
(r,
t)=
NI
(r,
t).
(2.29)
..@y
+
i@x
..i@t
+
mC1
The
resultant
probability
distribution
of
Eqn.
2.29
at
some
nonzero
time
may
be
seen
in
Fig.
2.2.
We
can
rewrite
the
spinor
components,
which
have
three
degrees
of
freedom
CHAPTER
2.
SOLVING
THE
DIRAC
EQUATION
(a)
Cx
=
Cy
=0;Cz
=
1
(b)
Cx
=1;Cy
=0;Cz
=0
Figure
2.2:
Twodimensional
wavepackets
with
localization
parameter
a
=
:2
at
time
t
=
1.
(two
real
parts,
two
imaginary
parts,
minus
the
condition
that
the
spinor
be
normalized),
in
terms
of
three
new
components
Cx,
Cy,
and
Cz.
Cx
=
C0
*
C1
+
C1
*
C0
Cy
=
i
(C0
*
C1

C0C1
*
)
Cz
=
jC0j2
..jC1j2
(2.30)
These
new
spinor
components
have
the
advantage
of
being
more
physically
meaningful,
as
Cx;Cy,
and
Cz
correspond
to
the
spinor
's
direction
in
threedimensional
space.
As
can
be
seen
in
Fig.
2.2(a),
allowing
Cz
=
1,
i.e.
setting
either
of
the
original
spinor
components
to
zero,
results
in
a
radially
symmetric
wave
packet.
Conversely,
setting
Cz
=
0,
resulting
in
the
original
spinor
pointing
somewhere
on
the
xy
plane
of
the
Bloch
sphere,
creates
a
maximally
asymmetric
wavepacket
as
can
be
seen
in
2.2(b).
Thus
the
magnitude
of
Cz
corresponds
to
the
degree
of
asymmetry
of
the
wave
packet,
and
the
magnitudes
of
Cx
and
Cy
determine
the
direction
of
the
asymmetry
about
the
zaxis.
Fig.2.3
conrms
the
use
of
a
as
a
localization
parameter.
For
smaller
values
of
a
the
wavepacket
is
increasingly
initially
localized.
For
values
of
a
greater
than
one,
the
wavepacket
spreads
similarly
to
a
Gaussian,
nonrelativistic
wavepacket,
but
when
a
is
much
less
than
one,
the
wavepacket
develops
in
a
highly
relativistic
fashion.
As
can
be
seen
in
Fig.2.2,
where
the
intial
value
of
a
is
much
smaller
than
one,
the
packet
has
a
characteristic
bowl
like
pattern
with
most
of
the
wavepacket
found
near
the
light
cone.
This
is
equivalent
to
a
revolution
about
the
zaxis
of
the
onedimensional,
twohorned
wavepacket
described
in
Chapter
1.
2.2.2
Three
Dimensions
Beginning
again
with
Eq.
2.4
we
can
write
the
components
of
the
3D
wavepacket.
First,
the
matrix
corresponding
to
the
the
positive
energy
projection
of
the
Hamiltonian
may
be
2.2.
WAVEPACKETS
WITH
SPINOR
COMPONENTS
(a)
a
=0:1
(b)
a
=1:1
(c)
a
=2:1
Figure
2.3:
Localization
of
the
2D
wave
packet
for
three
values
of
the
localization
parameter
a
at
time
t
=
0.
written:
..i@t

H^
=
..i@t
I4

(..i(@x
+
@y
+
@z)
·
a
+
m)
1
0
=
BB
.
..i@t
+
i@z
@y
+
i@x
..m
0
..@y
+
i@x
..i@t

i@z
0
..m
..m
0
..i@t

i@z
..@y

i@x
0
..m@y

i@x
..i@t
+
i@z
CC
.
.
(2.31)
We
also
allow
the
the
intial
spinor
to
be
in
an
arbitrary
state,
.
.
f
=
BB
.
C0
C1
C2
C3
CC
.
.
(2.32)
Combining
Eqs.
2.31
and
2.32
with
our
result
for
I
(r,
t),
and
taking
the
prescribed
derivatives
results
in
our
solution
for
the
3D
wavepacket.
The
resultant
probability
distribution
is
found
to
be
2
..
1
.
*
=N
2
m
jK+j2
r
2
+
ja+j2
t
(x
Re(Cx)

y
Im(Cy)+
z
Re(Cz))
4
4
s
..
*
K
*
1
K+
(+

x
Im(Cx
)+
y
Re(Cy)

z
Im(Cz
+
s
(2.33)
))

x
Im(Cx)+
y
Re(Cy)

z
Im(Cz)
+
sK1K+*
a
*
+
+
jsj2jK1j2,
CHAPTER
2.
SOLVING
THE
DIRAC
EQUATION
where
we
have
let
Cx
=
C0
C3
+
C1
C2
(2.34)
Cy
=
C0
C3

C1
C2
(2.35)
Cz
=
C0
C2

C1
C3
(2.36)
Kn
=
Kn(ms)
(2.37)
K+
=
K0
+
K2
+2K1
.
(2.38)
ms
As
in
the
twodimensional
case,
we
can
choose
to
rewrite
the
spinor
components
in
a
more
suggestive
basis.
We
can
also
dene
a
parameter
analogous
to
Cz
in
the
2D
case,
Ca
=
jC0j2+jC1j2..jC2j2..jC3j2
.
(2.39)
A
relationship
between
the
Cx;y;z
spinor
components
and
their
respective
axes
can
be
noted
in
Eq.
2.38.
Each
instance
of
Cx;y;z
is
accompanied
by
its
respective
variable,
x,
y,
or
z,
in
the
probability
distribution.
Thus,
altering
the
values
of
Cx,
Cy,
and
Cz
changes
the
asymmetry
axis
of
the
packet.
Increasing
Cx
(Cy;Cz)
points
the
asymmetry
closer
towards
the
^x
(^y,
z^)
direction.
When
either
the
rst
or
second
two
spinor
components
is
zero,
Ca
=
1,
and
Cx;y;z
=
0.
In
the
case
that
C0
=
1,
C1;2;3
=
0
or
Ca
=
1,
Cx;y;z
=
0
the
probability
distribution,
Eq.
2.33,
reduces
to
j j2
=
N
2
2
+
r
2
mK+
2s2
2
mK1
K+
+
m(a
+
it)
.
(2.40)
2s2
s
All
nonsymmetric
components
are
zero,
so
the
value
of
the
wavepacket
at
a
particular
point
depends
only
on
the
size
of
the
localization
parameter,
the
time,
and
the
distance
of
that
point
from
the
origin.
Fig.
2.4
(c)
shows
such
a
packet.
Ca
determines
the
degree
of
asymmetry
of
propagation.
As
in
two
dimensions
with
Cz
=
1,
when
Ca
=
1,
the
wavepacket
propagates
symmetrically,
with
a
uniform
wavefront.
Allowing
Ca
to
be
zero
results
in
a
maximally
asymmetric
wavepacket.
2.2.
WAVEPACKETS
WITH
SPINOR
COMPONENTS
(a)
Cx
=
1,
Cy
=
=
0,
(b)
Cx
=
:5,
Cy
=
Cz
=
0,
(c)
Cx
=
0,
Cy
=
=0,
Cz
v
Cz
3
Ca
=0
Ca
=
Ca
=1
2
Figure
2.4:
Density
plot
of
3D
wavepacket
projected
onto
the
xy
plane
at
t
=
1,
a
=
:2
for
various
initial
spinor
values.
CHAPTER
2.
SOLVING
THE
DIRAC
EQUATION
Chapter
3
Moments
of
Dirac
Equation
Solutions
In
this
Chapter
we
derive
moments
for
our
constructed
twoand
threedimensional
wavepackets.
We
discuss
the
results
in
the
context
of
the
packets’
evolution
and
compare
the
results
to
the
probability
distributions
calculated
in
Chapter
2.
3.1
2D
Moments
3.1.1
First
Moments
The
expectation
value
of
x,
written
as
R1
hx)
=
N
2..8
(x,
y,
t)x (x,
y,
t)dxdy,
(3.1)
where
we
have
previously
dened
.
(in
Eqs.
2.4
and
2.5)
as
Z
..+!(k)
(r;t)=
N..i@t

H^
1
ee
ikrdkv,
(3.2)
2p
v
!(k)
may
be
most
easily
calculated
in
momentum
space,
where
the
position
operators
x
=
i@=@px
and
y
=
i@=@py.
We
can
then
rewrite
3.1
as
Z8
hx)
=
N
2F
*
(px;py;t)i
.
(px;py;t)dpxdpy.
(3.3)
@px
..8
The
Fourier
tranform
of
,
,
may
be
readily
identied
with
part
of
the
integrand
of
Eq.
3.2,
where
Z8
(r)=(k)e
ikrdkv.
(3.4)
..8
Plugging
F
into
Eq.
3.3
yields
Z8
^^
H†
+
.
dH
+
!
..!yi
hx)
=
N
2ee
..!f
dpxdpy.
(3.5)
.
dpx
!
..8
19
CHAPTER
3.
MOMENTS
OF
DIRAC
EQUATION
SOLUTIONS
Taking
the
derivative
in
Equation
3.5
and
rewriting
in
terms
of
the
previously
dened
Cx
and
Cy
we
have
(with
assistance
from
Mathematica),
Z8
..
..2!2
hx)
=
N
2!..2
eCym

2tCxpxdpxdpy.
(3.6)
..8
As
shown
in
the
Appendix,
this
integral
may
be
evaluated
in
terms
of
upper
incomplete
gamma
functions,
Z8
tn..1
..n
=e
..tdt,
(3.7)
2ma
resulting
in
x
and
y
moments:
hx)
=
N
2
mCy..2

tmCx
1
..0

..1,
(3.8)
22
4m
1
hy)
=
N
2
mCx..2

tmCy22
..0

..1.
(3.9)
4m
While
the
time
dependent
part
of
the
x
moment
does
depend
on
Cx,
the
stationary
part
depends
on
Cy
pointing
to
the
inverse
relationship
between
initial
packet
size
and
speed
of
spreading
noted
in
Chapter
2.
In
the
case
that
Cx
=
Cy
=
0
we
can
clearly
see
that
the
rst
moments
are
zero
too,
resulting
in
symmetric
propagation
of
the
packet
with
respect
to
the
x
and
y
axes.
Looking
at
the
hxy)
moment,
which
depends
on
Cx
and
Cy
in
much
the
same
way
at
hxi,
conrms
that
the
wavepacket
will
develop
fully
symmetrically
if
Cx
=
Cy
=
0.
Though
neither
of
the
moments
depend
explicitly
on
Cz,
since
Cx
2
+
Cy
2
+
Cz
2
=
1
the
asymmetry
parameter
shows
up
implicitly.
Additionally,
besides
the
rst
moments’
dependence
on
the
scaling
parameter
a
through
..n,
there
is
also
an
explicit
a
dependence.
Smaller
values
of
a
increase
the
apparent
asymmetry
of
the
wavepacket
with
time,
as
a
piece
of
the
time
dependent
component
of
the
moments
depends
on
the
inverse
of
2
.
3.1.2
Second
Moments
We
use
the
same
procedure
for
deriving
the
second
moments
as
we
used
for
the
rst,
changing
only
x
x2
.
R8
R1
hx
2)
=
N
2..1..8
(x,
y,
t)x2 (x,
y,
t)dxdy
R8
R8
=
R
N
2(px;py;t)
@2
x
(px;py;t)dpxdpy
..1..8
@p2
e2
22
=pv
..!2
5
!a
((!2A

pxB)(.
+
Czm)

(A
+
1)pxCzm
+2t2px!2(.
+
Czm))
dpv,
(3.10)
where
we
have
dened
the
unitless
parameters:
A
=1+2!,
B
=
A
+22!2
.
(3.11)
3.2.
3D
MOMENTS
Upon
exchanging
px
for
py,
Eq.
3.10
is
identical
to
that
of
hy2i.
Since
Cz
does
not
distinguish
between
the
x
and
y
directions,
appearing
in
the
x2
and
y2
moments
in
the
same
way,
it
has
no
eect
on
the
direction
of
the
propagation
of
the
wavepacket.
Furthermore,
the
second
moments
do
not
depend
on
either
Cx
or
Cy.
This
conrms
our
previous
observation
that
Cx
and
Cy
do
not
appear
to
aect
the
amount
the
distance
the
wavepacket
travels,
only
its
direction.
Upon
taking
the
integral
in
Eq.
3.10
we
arrive
at
our
nal
expression
for
the
second
moments:
(2)
X(2)
2
(2)
V
(2)
hx
2)
=
X+
Cz+
tV
+
Cz.
(3.12)
0
z
0
z
The
second
moments
may
be
broken
up
in
two
ways,
rst
according
to
whether
or
not
the
component
depends
on
time
and
second
based
on
whether
or
not
the
component
depends
on
Cz.
..
X0(2)
=
1
2
..0
+..1
+..2
+22m2
(....2
+
2..3
+
2..4)(3.13)
..
X(2)
z
=2m..2
1
..1

..3
+22m2
(..3
+
4..4
+
6..5)(3.14)
(2)
V0
=
22
(..0

42m2..2)
(3.15)
V
(2)
m
z
=
a
(..1

42m2..3)
(3.16)
Interestingly,
although
the
2D
wavepacket
has
no
closed
form
solution,
both
the
rst
and
second
moments
may
be
written
entirely
in
terms
of
gamma
functions.
The
second
moments
are
entirely
dependent
on
the
localization
parameter
a
and
the
third
spinor
component
Cz.
This
is
similar
to
the
1D
case
explored
in
Strauch
[14]
wherein
two
of
the
spinor
components
may
be
found
in
the
rst
moment
and
the
third
spinor
component
is
only
found
in
the
second
moment.
3.2
3D
Moments
The
procedure
for
nding
moments
in
3D
is
similar
to
that
of
the
2D
case
with
only
a
change
in
the
variable
of
integration.
This
change
yeilds
modied
Bessel
functions
and
their
integrals
as
being
the
general
solutions
to
components
of
the
moments
rather
than
the
upper
incomplete
gamma
functions
in
two
dimensions.
3.2.1
First
Moments
We
begin
with
Eq.
3.1,
and
after
applying
the
derivative
we
nd
the
integral
to
be
similar
to
that
of
Eq.
3.6.
e..2!a
..
hx)
=
N
2
!2..2mIm
(Cx)+4tRe
(Cx)
p
2
xdpv
(3.17)
pv
CHAPTER
3.
MOMENTS
OF
DIRAC
EQUATION
SOLUTIONS
As
can
be
seen
in
the
Appendix,
taking
the
integral
of
Eq.
3.17
and
analogous
integrals
for
y
and
z
results
in
rst
moments
which
are
combinations
of
Bessel
functions:
....
hx)
=
N
2m2Im(Cx)(K1

Ki2)+4t
Re(Cx)m3Ki1

4
5
K1
+
1
4
K3,
(3.18)
....
N
231
hy)
=
..m2Re(Cy)(K1

Ki2)+4t
Im(Cy)m..Ki1

4
5
K1
+
4
K3,
(3.19)
..
hz)
=
N
2m2Im(Cz)(K1

Ki2)+4t
Re(Cz)m3Ki1

5
K1
+
1
,
(3.20)
44
K3
where
the
normalization
constant,
easily
calculable
in
momentum
space,
is
found
to
be
N
2
a
..1
=(K2(2m)+
CK1(2m)).
(3.21)
4m2
As
in
the
2D
case,
the
moments
have
components
which
either
do
not
depend
on
or
depend
linearly
on
time.
These
same
combinations
of
spinor
components
appeared
in
the
2.33,
the
expectation
value
of
the
threedimensional
wavepacket.
Predictably
Cx,
which
in
the
wavepacket
only
appears
in
terms
which
also
depend
on
x,
only
appears
in
the
x
moment.
The
same
is
true
of
the
other
variable
y
and
z.
The
mixed
moments,
hxyi,hyzi,
and
hzx)
also
depend
on
unique
combinations
of
of
the
spinor
components.
Cxy
=
jC0j2..jC1j2+jC2j2..jC3j2
(3.22)
Cyz
=
C0
C1
+
C1
C0
+
C2
C3
+
C3
C2
(3.23)
Czx
=
CC1

CC0
+
CC3

CC2
(3.24)
0123
3.2.2
Second
Moments
The
integral
representation
of
the
second
moments
may
be
written,
Z
e..2!a
....
2
222
hx
2)
=
!5!2A

pxB(.
+
Cm)

(A
+
1)pxCm
+2tpx!2(.
+
Cm)dpv.
pv
(3.25)
Ca
is
maximized
when
either
the
rst
two
or
second
two
spinor
components
is
zero.
This
is
identical
to
requiring
that
Cx
=
Cy
=
Cz
=
0
or
that
the
wavepacket
propagate
fully
symmetrically.
Also
note
that
Eqs.
3.10
and
3.25
are
identical
save
for
dierences
in
their
spinor
components
and
the
area
over
which
the
integral
is
being
taken.
We
can
now
write
the
second
moment
as
2)
=
hy
2)
=
hz
2)
=
N
2(2)
+
CX(2)
2(2)
+
CV
(2)
hxX0
a
+
tV0
,
(3.26)
where
we
have
dened
the
following
components:
41
X0(2)
=
m2K1

Ki1

Ki3
+2m
(K2

Ki2)+
2
m
2
(5K1

K3

4Ki1),
(3.27)
32
4
..
X(2)
a
=
m3(Ki2

Ki4)+2m
(K1
+
Ki1

2Ki3)+
2
m
2
(3K0

K2

Ki2),
(3.28)
3
V0(2)
=
2
m3
(..5K1
+
K3
+4Ki1)
,
(3.29)
3
V
(2)
a
=
4
m3
(..3K0
+
K2
+2Ki2)
:(3.30)
3
3.2.
3D
MOMENTS
These
moments
are
plotted
as
functions
of
a
in
Figs.
3.1a
and
3.1b.
The
value
of
X0(2)
scales
with
a
2
for
large
values
of
a
conrming
its
use
as
a
localization
parameter
in
three
dimensions.
For
small
values
of
,
X0(2)
approaches
zero
much
faster
than
.
The
value
(2)
(2)
(2)
of
Xa
is
much
smaller
than
X0
for
a
given
value
of
.
In
the
large
a
limit,
Xa
only
(2)
(2)
(2)
approaches
X0
=2,
while
when
a
is
small,
X0
=Xa
approaches
innity.
(b)
(a)
Figure
3.1:
Time
independent
elements
of
hx2i,
rescaled
with
the
localization
parameter
.
The
xaxes
are
on
a
log
scale.
The
velocity
component
of
second
moment
approaches
1
3
for
small
values
of
,
increasingly
initially
localized
packets.
Adding
the
the
x,
y,
and
z
components
of
the
second
moment
velocity
term
gives
the
second
moment
velocity
term
for
r.
Thus,
for
increasingly
localized
packets,
the
velocity
expectation
value
of
r
approaches
one,
the
speed
of
light.
(2)
(2)
Va
is
at
a
maximum
when
a
is
one.
Like
the
Xa
term,
it
is
much
less
than
the
term
(2)
which
does
not
depend
on
.
For
values
of
a
much
less
than
one,
Va
approaches
zero
with
(2)
V0
.
This
conrms
the
observation
that
initially
wide
packets
spread
much
slower
than
localized
ones.
Figure
3.2:
The
time
dependent
elements
of
hx2i.
The
xaxes
are
on
a
log
scale.
CHAPTER
3.
MOMENTS
OF
DIRAC
EQUATION
SOLUTIONS
Chapter
4
Quantum
Cellular
Automata
In
this
Chapter
we
examine
numerical
simulations
of
QCAs
on
square
and
hexagonal
lattices.
We
begin
with
a
derivation
of
the
square
QCA
using
a
multistep
scheme
before
rewriting
it
as
a
temporally
homogeneous
QCA.
Using
the
square
QCA
as
a
guide
we
then
create
a
hexagonal
QCA
using
a
similar
scheme
before
also
transitioning
to
a
temporally
homogeneous
case.
We
discuss
the
status
of
the
New
NoGo
lemma
in
the
context
of
our
derivation
of
a
hexagonal
QCA.
4.1
Square
Cellular
Automata
We
want
to
construct
a
square
QCA
with
an
update
rule
like
(x,
y,
t
+1)
=
A++ (x

1;y

1;t)+
A+.. (x

1;y
+1;t)
+A..+ (x
+1;y

1;t)+
A.... (x
+1;y
+1;t),
(4.1)
where
the
2
×
2
A
matricesrecall
that
in
two
dimensions,
at
each
lattice
point
.
is
a
two
component
spinorpreserve
the
normalization
of
the
lattice.
Following
a
multistep
method
discussed
in
Mackay
[6],
we
can
dene
our
unitary
operator,
UT
otal
=
UyUx,
(4.2)
where
the
Ux
and
Uy
components
are
unitary
operators
governing
1D
transitions
along
the
x
and
y
axes,
respectively.
In
the
rst
time
step,
Ux
is
applied
to
lattice;
in
the
second
Uy
is
applied,
Fig.
4.1.
Drawing
inspiration
from
the
results
of
BialynickiBirula
[8],
we
can
write
these
steps
as
Un
=
1(I+
n)
.
Tn
+
1(I
n)
.
Tn
..1
;n
=
x,
y.
(4.3)
22
Un
is
composed
of
coin
operations
on
the
lattice
points
and
shift
operators.
Tn
translates
one
lattice
point
in
the
positive
n
direction,
while
Tn
..1
translates
in
the
negative
n
direction.
25
CHAPTER
4.
QUANTUM
CELLULAR
AUTOMATA
Figure
4.1:
The
two
translation
steps
for
a
square
cellular
automata.
The
lattice
updates
rst
along
the
xaxis,
then
the
yaxis,
alternating
back
in
forth
in
subsequent
steps.
The
matrices
n
can
be
any
two
of
the
Pauli
matrices,
but
we
have
chosen
x
to
act
along
the
xaxis
and
y
to
act
along
the
yaxis.
We
have
now
fully
described
the
square
multistep
quantum
walk.
This
method
does
not
obviously
yield
a
temporally
homogeneous
QCA,
though,
because
the
action
of
U
alternates
acting
along
the
xand
yaxes.
Fortunately
one
needs
only
to
multiply
the
Un
components
to
nd
a
temporally
homogeneous
QCA
of
the
form
in
Eq.
4.1.
We
nd
1+
i
1

i
1
4
1
4
A++
(4.4)
=
1+
i
1

i
1

i
1+
i
A+
(4.5)
=
1

i
1+
i
1

i
..1

i
1
4
1
4
A..+
(4.6)
=
..1+
i
1+
i
1+
i
..1+
i
A..
(4.7)
=
:
..1

i
1

i
As
shown
in
Fig.
4.2,
we
can
reduce
each
of
the
four
pairs
of
translation
operators
into
single
diagonal
translations.
That
is,
the
new
translation
steps
do
not
lie
along
the
primary
directions
as
the
rst,
but
on
a
square
lattice
rotated
by
45.
.
We
can
nally
check
if
the
unitarity
condition
of
this
QCA
is
satised.
UyU
=
UxyUy
yUyUx
=
I
(4.8)
Indeed,
upon
plugging
in
our
expressions
for
Un,
we
nd
the
condition
holds.
4.1.1
Numerical
Simulation
of
Square
Automaton
The
anisotropies
of
the
square
cellular
automata
are
evident
when
there
is
only
a
single
nonzero
cell
on
the
initial
lattice.
In
Fig.
4.3,
one
can
note
two
preferential
axes
along
the
4.1.
SQUARE
CELLULAR
AUTOMATA
Figure
4.2:
The
result
of
reducing
the
two
translation
steps
of
4.1
into
a
single
translation
step
on
two
axes
simultaneously.
The
intermediate
sites
have
been
eliminated.
Figure
4.3:
Square
QCA
after
50
time
steps.
At
time
t
=
0
only
the
centermost
cell
has
a
nonzero
value
where
C0
=
1,
C1
=
0.
CHAPTER
4.
QUANTUM
CELLULAR
AUTOMATA
directions
of
the
reduced
translation
vectors.
Along
the
y
=
..x
line
in
particular
are
clear
lines
of
increased
probability
which
fan
out
as
they
approach
the
edge
of
of
the
wavefront.
Additionally,
with
only
a
single
initially
nonzero
cell,
the
spreading
of
the
wavefront
does
not
look
relativistic.
The
probability
density
inside
the
wavefront,
even
at
the
initial
nonzero
cell
location,
is
of
the
same
order
of
magnitude
as
that
on
the
edge
of
the
wavefront.
However,
for
dierent
initial
states
one
can
produce
wavepackets
with
dierent
properties.
This
will
be
shown
for
the
hexagonal
QCA
below.
4.2
Hexagonal
Quantum
Cellular
Automata
As
discussed
in
Chapter
1,
D'Ariano
and
Perinotti
[2]
argued
that
there
are
no
strictly
spatially
and
temporally
homogeneous
QCAs
in
two
dimensions
besides
the
square
one.
One
of
the
simplest
conditions
which
we
may
relax
in
an
eort
to
create
a
hexagonal
QCA
is
that
of
temporal
homogeneity.
We
have
already
seen
how
a
multistep
QCA
on
a
square
lattice
is
equivalent
to
that
of
the
temporally
homogeneous
one.
Changing
from
the
translation
vectors
to
the
reduced
translation
vectors
eectively
converts
the
temporally
inhomogeneous
square
QCA
into
a
temporally
homogeneous
one.
This
property
of
square
QCAs
does
raise
questions
about
the
temporal
homogeneity
requirement
though.
Are
there
other
lattices
that
share
this
property?
Since
the
only
other
topologically
distinct
lattice
in
two
dimensions
is
the
hexagonal
lattice,
it
is
the
natural
choice
for
considering
this
question.
We
begin
by
dening
an
alternate
QCA
in
much
the
same
way
that
we
did
on
the
square
lattice.
First,
the
entire
unitary
operator
is
composed
of
three
other
operators
which
act
on
the
lattice
sequentially:
UT
otal
=
U3U2U1.
(4.9)
We
then
dene
translation
vectors
and
transformation
matrices,
taking
into
account
that
we
are
no
long
translating
along
the
x
or
y
axes.
To
that
end
we
can
construct
the
spin
matrices
as
linear
combinations
of
the
rst
and
second
Pauli
matrices,
and
the
translation
vectors
similarly.
v
v
13
1
=
1
2
x
+
2
3
y;2
=

x
+
y;3
=
x
(4.10)
22
v
v
1
3
1
3
..ia
kx+
ky
..ia..
kx+
ky..iakx
T1
=
e
2
2
,
T2
=
e
2
2
T3
(4.11)
=
e
,
Putting
everything
together,
we
have
the
general
form
of
an
operation
for
each
of
the
three
translation
steps:
Un
=
1(I+
n)
.
Tn
+
1(I
n)
.
Tn
..1
;n
=1,
2,
3.
(4.12)
22
The
hexagonal
Un
is
dened
identically
to
the
square
Un,
save
for
components
being
modied
to
handle
the
translations
not
being
along
the
primary
axes.
4.2.
HEXAGONAL
QUANTUM
CELLULAR
AUTOMATA
Figure
4.4:
Translation
steps
along
the
three
axes
of
a
hexagonal
lattice.
Also,
in
color
are
the
four
points
on
the
hexagonal
primitive
cell
which
support
four
noninteracting
QCAs.
As
can
be
seen
in
Fig.
4.4,
after
applying
the
three
translation
steps,
a
single
lattice
point
has
in
uenced
seven
lattice
points,
six
of
which
are
two
lattice
points
away
from
itself.
As
a
consequence
of
it
not
in
uencing
most
of
its
neighbors,
as
with
the
square
automaton
splitting
into
two
disconnected
lattices,
the
hexagonal
automata
splits
into
four
disconnected
lattices,
with
each
point
on
the
primitive
cell
on
a
dierent
eective
lattice.
We
can
now
dene
new
translation
vectors
T~1,
T~2,
and
T~3:
~
T1
=
T3T2T1
(4.13)
~T
..1
=
(4.14)
T23
T2T1
T~3
=
T3T2
..1T1
(4.15)
T~..1
T
..1T
..1T
..1
=
(4.16)
1
321
T~..1
T3T
..1T
..1
2
=
21
(4.17)
T~..1
T
..1T2T
..1
=
.
(4.18)
3
31
In
addition
to
those
reduced
translation
vectors,
we
must
also
dene
a
zero
translation
vector,
T
..1
T0
=
T1T2
..1
3
=
T1
..1T2T3
=
I.
(4.19)
We
have
thus
temporally
homogenized
the
action
of
the
QCA.
As
can
be
seen
in
Fig.
4.5,
the
nonzero
probability
of
not
moving
in
a
time
step
is
explicit
in
this
new
scheme.
Additionally,
unlike
in
the
case
of
the
square
automaton,
the
new
translation
vectors
of
the
hexagonal
automaton
lie
along
the
same
directions
as
the
original
translation
vectors,
no
rotation
necessary.
We
can
then
dene
the
full
action
of
the
unitary
operator
as,
rewriting
T~n
as
Tn,
U
=
A0
+
A1T1
+
A..1T
..1
+
A2T2
+
A..2T
..1
+
A3T3
+
A..3T
..1
.
(4.20)
123
We
expect,
for
this
to
be
considered
a
true
QCA,
for
U
to
be
unitary.
We
nd
that,
as
noted
in
the
Appendix,
U
is
unitary
for
the
rule
given
above
(this
can
also
be
seen
from
Eqns.
4.9
and
4.12).
CHAPTER
4.
QUANTUM
CELLULAR
AUTOMATA
Figure
4.5:
Reduced
translation
steps
along
the
three
axes
of
a
hexagonal
lattice.
4.2.1
New
NoGo
Lemma
We
have
found
a
counterexample
to
D'Ariano
and
Perinotti's
New
NoGo
Lemma,
in
the
process
developing
a
method
which
might
be
used
to
explore
QCAs
on
other,
previously
forbidden
lattices.
Where
did
their
argument
go
wrong?
In
D'Ariano
and
Perinotti
[2]
the
hexagonal
lattice
is
dismissed
as
a
possible
QCA
because
they
argue
that
the
unitary
conditions
cannot
be
satised
in
the
massless
case,
that
is,
one
with
the
inertial
A0
matrix
equal
to
0.
It
is
easy
to
see,
however,
that
allowing
a
nonzero
inertial
matrix
provides
more
free
parameters
for
satisfying
the
unitary
conditions.
As
shown
in
the
Appendix,
adding
an
inertial
term
also
does
not
add
any
more
equations
to
Eq.
4.20.
Thus,
adding
an
inertial
term
gives
more
options
for
satisfying
the
unitary
conditions
than
it
takes
aways.
This
indicates
that
even
if
the
argument
presented
in
[2]
is
correct,
the
New
NoGo
Lemma
could
only
ever
apply
to
Weyl
automata.
The
Dispersion
Relation
The
eigenvalues
of
U
give
the
dispersion
relation
for
the
hexagonal
QCA:
......v
..v
cos(!)=
8
1
..1
+
3cos(2kx)
+
coskx
+3ky+
coskx

3ky(4.21)
Ignoring
terms
on
the
order
of
k3,
in
the
long
wavelength
limit,
this
reduces
to:
9
..
!2
k2
+
k2
=
xy.
(4.22)
4
This
is
of
the
same
form
as
the
dispersion
relation
for
the
Weyl
equation,
but
with
a
dierence
in
the
numerical
factor
preceding
the
wavenumbers.
This
dierence
may
be
partly
attributed
to
the
length
in
lattice
units
of
each
reduced
translation
step2,
and
partly
to
the
nonzero
A0
matrix.
We
could
redene
either
the
lattice
units
or
the
time
units
to
give
exact
agreement
with
the
Dirac
equation.
4.2.
HEXAGONAL
QUANTUM
CELLULAR
AUTOMATA
Figure
4.6:
Hexagonal
QCA
after
50
time
steps
with
a
single
nonzero
cell,
with
spinor
component
Cz
=
1
at
time
t
=
0.
Plots
With
only
a
single
nonzero
cell
at
time
t
=
0,
the
preferred
directions
of
the
lattice
become
apparent,
just
as
in
the
case
of
the
square
automaton.
In
Fig.
4.6,
the
probability
density
is
highest
along
the
three
translation
axes,
at
the
wavefront,
and
near
the
initially
nonzero
cell.
Using
spinor
components
C0
=
1,
C1
=
0
yields
symmetric
propagation,
as
in
the
case
of
the
square
automaton.
A
rudimentary
wave
packet
can
be
constructed
by
allowing
the
some
of
the
cells
surrounding
the
origin
cell
to
also
have
nonzero,
though
small,
initial
magnitudes.
Giving
each
of
the
cells
the
same
spinor
components,
we
nd
that
the
wavepacket
propagates
remarkably
uniformly,
as
shown
in
Fig.
4.7.
The
anisotropies
of
the
lattice
have
been
washed
out,
and
there
is
no
longer
signicant
probability
in
the
center
of
propagation.
The
wavepacket
additionally
has
a
clear
and
uniform
wavefront,
though
it
is
double
edged.
Constructing
a
wavepacket
in
momentum
space,
as
we
did
with
the
Dirac
equation,
would
no
doubt
increase
the
similarities
between
the
packet
and
our
solution
for
the
Dirac
equation.
CHAPTER
4.
QUANTUM
CELLULAR
AUTOMATA
Figure
4.7:
Hexagonal
QCA
after
50
time
steps
with
an
initial
rudimentary
wavepacket.
Each
of
the
seven
initial
nonzero
cells
has
only
a
single
nonzero
spinor
component.
Chapter
5
Discussion
In
this
Chapter
we
discuss
the
implications
of
our
results
in
the
context
of
D'Ariano
and
Perrinotti's
New
NoGo
Lemma.
We
also
describe
future
work
to
be
done
with
QCAs
on
nonsquare
or
BCC
lattices
and
their
relationship
to
the
Dirac
equation.
5.1
The
New
NoGo
Lemma
We
have
found
the
QCAs
D'Ariano
and
Perrinotti
allow
in
two
dimensions
to
be
limited
by
their
denition
of
QCAs.
Besides
the
square
QCA
in
two
dimensions,
a
multistep,
massless
hexagonal
QCA
and
a
massive,
temporally
homogeneous
hexagonal
QCA
may
also
be
found.
By
not
considering
the
possibility
of
an
inertial
term
(the
A0
matrix
element)
in
their
derivation
of
the
possible
2D
QCAs,
D'Ariano
and
Perrinotti
fail
to
consider
the
hexagonal
QCA.
Thus,
we
have
provided
an
extended
framework
for
classifying
QCAs
on
higherdimensional
lattices,
and
for
which
the
New
NoGo
Lemma
does
not
apply.
5.2
Future
Work
We
discuss
extensions
of
our
work
with
the
Dirac
equation
and
its
connection
to
QCAs.
5.2.1
Extending
the
Dirac
Equation
Solution
Our
solution
for
the
Dirac
equation
took
the
form
of
Eq.
5.1,
where
g(k)
was
chosen
to
be
1=!(k).
Other
distributions
in
momentum
space
are
possible.
..!(k)
I
(r;t)=
g(k)ee
i(kzr..!(k)t)dkv
(5.1)
v
Indeed,
certain
values
of
g(k)
may
yield
closed
form
solutions
to
the
twodimensional
Dirac
equation.
33
CHAPTER
5.
DISCUSSION
Figure
5.1:
Lattice
structure
of
graphene.
Each
lattice
point
has
three
nearest
nieghbors.
5.2.2
Adding
Mass
We
have
considered
several
ways
to
incorporate
mass
in
the
hexagonal
and
square
QCAs
which
would
take
us
from
Weyl
to
Dirac
automata.
But,
the
problem
of
what
kind
of
mass
term,
inertial
or
chiral,
should
be
added
to
a
Weyl
QCA
to
give
the
Dirac
propagator
does
not
have
an
obvious
answer.
In
the
square
case,
besides
the
issue
of
disconnected
lattices,
the
two
types
of
mass
terms
seem
to
have
similar
functions.
For
the
hexagonal
automaton,
the
inertial
mass
term
is
integral
to
its
unitarity,
but
whether
another,
variable,
inertial
mass
term
may
be
added
to
the
automaton,
remains
open.
5.2.3
Other
Lattices
Another
open
question
is
how
to
extend
primitive
QCAs
to
other
lattices
in
two
and
three
dimensions,
such
as
the
kagome
and
graphene
lattices.
Although
they
are
not
true
lattices
in
the
mathematical
sense
that
a
set
of
translation
vectors
and
their
inverses
can
be
applied
from
a
single
starting
point
to
tile
the
plane,
they
are
natural
extensions
of
the
QCA
model.
In
the
square
and
hexagonal
lattices,
each
point
has
four
and
six
nearest
neighbors,
respectively.
However,
there
are
interesting
nonBravais
lattices
in
two
dimensions
with
dierent
numbers
of
nearest
neighbors,
such
as
the
graphene
lattice
in
Fig.
5.1.
These
lattices
may
yield
new
and
interesting
behavior.
Three
Dimensions
Most
of
the
numerical
work
for
this
project
is
in
two
dimensions
because
1)
twodimensional
simulations
are
easier
to
visualize,
especially
in
the
context
of
a
paper,
than
their
three
dimensional
analogues
and
2)
numerical
simulations
on
threedimensional
lattices
take
a
long
time
to
compute.
In
three
dimensions,
the
most
interesting
lattice
overlooked
by
D'Ariano
and
Perinotti
is
the
rhombohedral,
which
looks
like
stacked
hexagonal
lattices
connected
vertically,
but
the
simple
cubic
and
facecentered
cubic
lattices
were
also
ruled
out
by
the
New
NoGo
lemma.
Further
work
should
be
done
to
see
if
inertial
terms
can
be
introduced
5.2.
FUTURE
WORK
to
construct
primitive
or
other
natural
QCAs
for
higherdimensional
lattices.
5.2.4
Analytic
Connections
of
QCAs
to
the
Dirac
Equation
Further
work
should
also
be
done
to
show
analytically
that
the
hexagonal
QCA
and
our
solution
to
the
Dirac
equation
are
equivalent,
or
under
what
circumstances
they
may
be
considered
equivalent.
Numerically,
QCA
wavepackets
should
be
constructed
in
momentum
space
as
for
our
2D
Dirac
wavepacket
in
Chapter
2.
The
propagation
of
these
QCA
wavepackets
could
then
be
directly
compared
to
that
given
by
our
solution
to
the
Dirac
equation.
Using
either
of
these
methods
could
provide
better
information,
besides
simply
conditions
for
symmetric
propagation,
on
the
relationship
between
the
initial
spinor
in
our
Dirac
solutions
and
that
of
our
described
automata.
CHAPTER
5.
DISCUSSION
Appendix
A
Appendix
A.1
Moment
Calculations
From
chapter
3,
we
detail
the
calculations
for
the
2D
and
3D
moments
of
our
solution
to
the
Dirac
Equation.
A.1.1
Two
Dimensions
To
solve
an
expression
of
the
form
ZR8
!..nk
..2!dp2
p0e,
(A.1)
..8
with
n
an
integer,
k
either
0
or
2,
and
p0
representing
either
px
or
py,
we
rst
change
into
cylindrical
coordinates.
This
yields
an
expression
independent
of
our
choice
of
p0,
(
R8
2k
=0
1+k!..n..2atdr
0
re(A.2)
k
=2
(
R8
2k
=0
k!1..n..2atd.
=re.
(A.3)
m
k
=2
The
variable
of
integration
was
changed
from
r
to
.
in
the
second
expression.
Upon
making
the
substitution
.
=
t=2a
and
recalling
(in
the
k
=
2
case)
that
r2
=
!2

m2,
this
integral
may
readily
be
identied
with
that
of
the
incomplete
gamma
function
[15].
ZR8
ta..1
..(a,
x)=
e
..tdt
(A.4)
x
In
this
case,
the
x
in
A.4
is
2ma
while
a
is
an
integer
ranging
from
3
to
2.
37
APPENDIX
A.
APPENDIX
In
the
k
=
0
case
of
A.3,
we
have
Z8
Z8
..n+1
t
dt
!..n+1
..t
..2!d.
=
ee
(A.5)
2
m
2m2a
..n+2Z8
..t
=
1
et(2..n)..1dt
=
(2)n..2
..(2

n,
2m)
(A.6)
2
2ma
In
the
k
=
2
case
or
A.3,
Z8
Z8
..
2!..n+1
!22!..n+1
..2!d!
..2!d.
=
re
+
me
(A.7)
m
m!Z8
..n+3
..n+1
tt
=
m
2e
..tdt
(A.8)
22
2ma
=
(2)n..4
..(4

n,
2m)

m
2(2)n..2
..(2

n,
2m)
(A.9)
We
have
derived
expressions
for
the
solutions
of
each
component
of
the
2D
moment
integrals.
A.1.2
Three
Dimensions
We
want
to
solve
an
expression
identical
to
that
of
A.1,
except
in
three
dimensions.
First,
we
change
the
expression
to
spherical
coordinates
to
nd,
R8
!..nk
..8
pxe..2!dp3
(A.10)
Rv
km2..r22
=V
sin()r2
(r
sin()
cos())e..2a
(m2
+
r2)
n
drddf
(A.11)
Rp
R2p
R1v

n
2+k..2m2+r2
+
r2)2
=sink+1()dcosk()dre2
(mdr
(A.12)
0
00
We
can
now
solve
the
angular
integrals.
If
k
=
0,
the
angular
part
of
the
integral
is
4,
and
when
k
=2,
it
is
4=3.
Making
the
substitution
r
=
m
sinh(y),
we
nd
Z8
..2m
cosh(y)(m
=(m
sinh(y))2+k
e
cosh(y))..n+1dy.
(A.13)
0
If
k
=
2,
we
next
rewrite
the
sinh
component
in
terms
of
cosh.
This
set
of
integrals
may
then
be
identied
with
that
of
the
modied
Bessel
function
of
the
second
kind
[15].
In
table
A.1,
we
tabulate
our
nal
solutions
in
terms
of
the
power
of
k
and
n
in
Eq.
A.10.
A.2
Hexagonal
Automaton
Unitarity
Conditions
We
discuss
the
unitary
condition
on
the
hexagonal
automata,
mentioned
in
chapter
4,
in
two
dierent
cases.
The
hexagonal
automaton
must
satisfy
the
unitarity
condition,
UyU
=1.
(A.14)
A.2.
HEXAGONAL
AUTOMATON
UNITARITY
CONDITIONS
k
=
0
k
=
2
n
=
1
n
=
2
n
=
3
n
=
4
n
=
5
m2
2
(K0

K2)
m
(K1

Ki1
)
(K0

Ki2
)
m3....5
4
K1
+
1
4
K3
+
Ki1m2....3
2
K0
+
1
2
K2
+
Ki2m
(K1

2Ki1

Ki3
)
K0

2Ki2

Ki4
Table
A.1:
Solutions
to
equations
of
the
form
A.10
For
the
hexagonal
automaton
in
the
reduced
translation
scheme,
the
most
general
denition
of
a
unitary
operator
on
the
lattice
is,
U
=
A0
+
A1T1
+
A..1T1
..1
+
A2T2
+
A..2T2
..1
+
A3T3
+
A..3T3
..1
,
(A.15)
which
we
have
also
dened
in
Eq.
4.20.
Multiplying
this
by
its
complex
conjugate
yeilds:
UyU
=
jA0j2
+
jA1j2
+
jA..1j2
+
jA2j2
+
jA..2j2
+
jA3j2
+
jA..3j2
(A.16)
+A†
0A1
+
A..
†
1A0T1
+A†
0A..1
+
A1†
A0T1
..1
+AyA2
+
A†
A0T2
+AyA..2
+
A†
A0T
..1
0..2022
A†
A†
T
..1
+0A3
+
A†
T3
+0A..3
+
A†
(A.17)
..3A03A03
+Ay..3A2
+
Ay..2A3T2T3
+A†
3A..2
+
A2yA3T2
..1T
..1
3
+Ay..1A..3
+
A3yA1T1T
..1
+A†
1A3
+
A†
T
..1T3
3
..3A..11
+AyA2
+
A†
T
..1
A†
A..2
+
A†
T
..1
(A.18)
1..2A..11
T2
+..12A1T12
+A..
†
1A3
+
A†
T1T3
+A†
3A..1T1
..1T3
..1
..3A11A..3
+
A†
y†
y†
T
..1T
..1
+AT1T2
+A
..1A2
+
A..2A11A..2
+
A2A..112
+Ay..2A..3
+
A†
3A2T2T3
..1
+A†
2A3
+
Ay..3A..2T2
..1T3
(A.19)
yy
T
2
T
..2
+A..1A11
+A1A..11
A†
T
2
A†
T
..2
+..2A22
+2A..22
+A..
†
3A3T3
2
+A3yA..3T3
..2
(A.20)
APPENDIX
A.
APPENDIX
For
the
hexagonal
lattice,
the
elements
in
A.18
can
be
rewritten
using
the
following
identities:
T1
=
T2T3
=
T3T2
(A.21)
T2
=
T1T3
..1
=
T3
..1T1
(A.22)
T3
=
T1T2
..1
=
T2
..1T1.
(A.23)
Thus
combining
the
appropriate
elements
of
A.17
and
A.18,
we
have
our
equations
of
contraint
on
the
A
elements.
Those
of
A.16
must
equal
the
identity
and
each
of
the
combinations
of
A
elements
multiplied
by
translation
vectors
must
be
individually
zero.
UyU
=
jA0j2
+
jA1j2
+
jA..1j2
+
jA2j2
+
jA..2j2
+
jA3j2
+
jA..3j2
(A.24)
)
+
0A1
+
Ay..3A2
+
Ay
Ay..1A0
+
Ay..2A3
T1
)
†
y†
†
T
..1
+
A0A..1
+
A1A0
+
A3A..2
+
A2A31
)
+
A0†
A2
+
Ay..1A..3
+
A†
3A1
..2A0
+
A†
T2
)
†
yy†
T
..1
+
A0A..2
+
A2A0
+
A1A3
+
A..3A..12
)
+
A†
A3
+
A†
A0
+
A†
A..2
+
AyA1
T3
0..3..12
)
†
yy†
T
..1
+
A0A..3
+
A3A0
+
A1A2
+
A..2A..13
(A.25)
(
(
)
+
Ay..1A3
+
Ay..3A1
T1T3
+
A†
1A..3
+
A†
3A..1
T1
..1T3
..1
T
..1
+(::.
)T1T2
+(::.
)T
..1
+(::.
)T2T
..1
+(::.
)T
..2T3
(A.26)
12
32
(
)
Ay
+
A†
T
2
+
T
..2
..1A11
1A..11
(
(
)
A†
T
2
A†
T
..2
+
A2
+
A..2
..2222
)
T
..2
+
A†
T
2
+
A†
(A.27)
..3A33
3A..33
The
unitarity
requirement
gives
19
constraints
on
the
56
real
matrix
elements
of
the
seven
A
matrices.
Without
the
inertial
mass
term
A0,
there
are
the
same
number
of
constraint
equations,
19,
but
eight
fewer
matrix
elements,
48.
Bibliography
[1]
D.
A.
Meyer,
J.
Stat.
Phys.
85,
551
(1996).
[2]
G.
M.
D'Ariano
and
P.
Perinotti,
Phys.
Rev.
A
90,
062106
(2014).
[3]
C.
M.
Chandrashekar,
Scientic
Reports
3,
2829
(2013).
[4]
G.
M.
D'Ariano,
N.
Mosco,
and
P.
Perinotti,
Europhys.
Lett.
109,
40012
(2015).
[5]
G.
J.
de
Valcarcel,
M.
Hinarejos,
E.
Roldan,
A.
Perez,
and
A.
Romanelli,
New
J.
Phys
15,
073041
(2013).
[6]
T.
D.
Mackay,
S.
D.
Bartlett,
L.
T.
Stephenson,
and
B.
C.
Sanders,
J.
of
Phys.
A
35
(2002).
[7]
E.
Roldan,
C.
D.
Franco,
F.
Silva,
and
G.
J.
de
Valcarcel,
Phys.
Rev.
A
87,
022336
(2013).
[8]
I.
BialynickiBirula,
Phys.
Rev.
D
49
(1994).
[9]
S.
R.
Beane,
Z.
Davoudi,
and
M.
J.
Savage,
Eur.
Phys.
J.
A50,
148
(2014).
[10]
B.
Thaller,
The
Dirac
Equation
(SpringerVerlag,
New
York,
1992).
[11]
S.
Wolfram,
A
New
Kind
of
Science
(Wolfram
Media,
Champaign,
IL,
2002).
[12]
D.
A.
Meyer,
Phys.
Rev.
E
55,
5261
(1997).
[13]
F.
W.
Strauch,
Phys.
Rev.
A
73,
054302
(2006).
[14]
F.
W.
Strauch,
J.
Math.
Phys.
48,
082102
(2007).
[15]
Modied
bessel
functions
(2016),
URL
http://dlmf.nist.gov/10.32.
[16]
F.
W.
Strauch,
Phys.
Rev.
A
74,
030301
(2006).
[17]
G.
M.
D'Ariano,
N.
Mosco,
P.
Perinotti,
and
A.
Tosini,
Entropy
18,
228
(2016).
[18]
R.
Feynman
and
A.
R.
Hibbs,
The
Dirac
Equation
(McGrawHill,
New
York,
1965).
41