2006
, published 8 October
, doi: 10.1098/rspa.2006.1701
462 2006 Proc. R. Soc. A
M.W Janowicz, J.M.A Ashbourn, Arkadiusz Orlowski and Jan Mostowski
wave propagation in dispersive media
Cellular automaton approach to electromagnetic
Supplementary data
62.2074.2927.DC1.html
http://rspa.royalsocietypublishing.org/content/suppl/2009/03/09/4
"Data Supplement"
References
html#reflist1
http://rspa.royalsocietypublishing.org/content/462/2074/2927.full.
This article cites 26 articles
Email alerting service
herethe box at the top righthand corner of the article or click
Receive free email alerts when new articles cite this article  sign up in
http://rspa.royalsocietypublishing.org/subscriptions go to: Proc. R. Soc. ATo subscribe to
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
Cellular automaton approach to
electromagnetic wave propagation
in dispersive media
B
Y
M.W.J
ANOWICZ
1,3,
*
,J.M.A.A
SHBOURN
2
,A
RKADIUSZ
O
RŁOWSKI
3
AND
J
AN
M
OSTOWSKI
3
1
Fakulta
¨
t 5,Institut fu
¨
r Physik,CarlvonOssietzky Universita
¨
t,
26111 Oldenburg,Germany
2
Department of Engineering Science,University of Oxford,Parks Road,
Oxford OX1 3PJ,UK
3
Instytut Fizyki PAN,Al.Lotniko´w 32/46,02668 Warsaw,Poland
Extensions of BiałynickiBirula’s cellular automaton are proposed for studies of the one
dimensional propagation of electromagnetic ﬁelds in Drude metals,as well as in both
transparent,dispersive and lossy dielectrics.These extensions are obtained by
representing the dielectrics with appropriate matter ﬁelds,such as polarization together
with associated velocity ﬁelds.To obtain the different schemes for the integration of the
resulting systems of linear partial differential equations,splitoperator ideas are
employed.Possible further extensions to twodimensional propagation and for the
study of lefthanded materials are discussed.The stability properties of the cellular
automaton treated as a difference scheme are analysed.
Keywords:cellular automata;electromagnetic waves;propagation;dispersive media
1.Introduction
This paper is concerned with the simulation of light propagation in dispersive
materials with the help of suitably chosen cellular automata (CA).CA in a broad
sense form an idealization of a physical system in which space and time are
discrete.A multitude of applications has already been found for these,including
dynamical systems theory,formal languages,statistical mechanics,theoretical
biology and neural networks (cf.Ilachinski (2001) and references therein).
Techniques of the socalled latticegas CA,as developed in Hardy et al.(1976)
and Frisch et al.(1986),have been successfully applied to several systems of
partial differential equations (cf.Rothman & Zaleski (1994) and Chopard &Droz
(1998)).The fact that there have only been a few attempts to adapt the CA
related methods in order to model the electromagnetic wave propagation could
Proc.R.Soc.A (2006) 462,2927–2948
doi:10.1098/rspa.2006.1701
Published online 18 April 2006
The electronic supplementary material is available at http://dx.doi.org/10.1098/rspa.2006.1701 or
via http://www.journals.royalsoc.ac.uk.
*Author and address for correspondence:Instytut Fizyki PAN,Al.Lotniko
´
w 32/46,02668
Warsaw,Poland (mjanow@ifpan.edu.pl).
Received 2 December 2005
Accepted 21 February 2006
2927
q 2006 The Royal Society
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
probably be attributed to the great efﬁciency of more traditional methods,such
as ﬁnite differences or ﬁnite elements in problems of this kind.Among the most
widespread methods for the integration of the Maxwell equations in dispersive
and lossy media (which are of most interest to us),in addition to ﬁnitedifference
timedomain methods,there are approaches based on fast Fourier transformation
(e.g.Dvorak & Dudley 1995,1996),various numerical integration methods,
asymptotic techniques (e.g.Oughstun & Sherman 1989a,b,1990;Oughstun &
Balictsis 1996) and hybrid analytical–numerical methods,as described,for
example,in Dvorak et al.(1998).Nevertheless,the standard Boltzmann lattice
gas methods have been applied to Maxwell’s equations in three spatial
dimensions in Simons et al.(1999).
In two other approaches to CAbased integrators of the wave equations,the
authors decided to relax one of the most important deﬁning requirements of CA,
that of the discreteness of the dependent variables.In Chopard et al.(1997) and
Chopard & Droz (1998),the Boltzmann method was still retained.On the other
hand,in Sornette et al.(1993),the Smatrix formulation of the wave propagation
was used,which turned out to be equivalent to a ﬁnitedifference discretized
version of various wave equations.In addition,Vanneste &Sebbah (2001) as well
as Sebbah & Vanneste (2002) have discussed numerical propagation in random
dispersive media,which is very much in the spirit of our present work.In this
paper,we develop the approach initiated by BiałynickiBirula (1994),who has
been mostly interested in problems of a theoretical nature.One of our aims is to
demonstrate that in the case of onedimensional propagation at least,the
BiałynickiBirula algorithm (BBA) is also practically very useful.We show this
by extending it to the problem of the dynamics of signals in dispersive materials.
2.BiałynickiBirula’s cellular automaton in (1D1) dimensions
Let us start with the ﬁrst timedependent pair of Maxwell’s equations in
a homogeneous dispersionless dielectric in the absence of charges and currents
(in SI units):
V!E ZK
vB
vt
;V!BZm
0
e
0
me
vE
vt
;ð2:1Þ
where e and m are the dielectric constant and the magnetic permeability of
the dielectric,respectively.In one spatial dimension,it is sufﬁcient to represent
the electric ﬁeld vector as EZð0;0;EðxÞÞ and the magnetic ﬁeld vector as
BZð0;BðxÞ;0Þ (the other polarization may be considered in an analogous way).
It is convenient to introduce a new variable F,such that EZcF.
Then Maxwell’s equations take the form
vF
vt
Z
c
n
2
vB
vx
;
vB
vt
Zc
vF
vx
;ð2:2Þ
where n
2
Zem.If we introduce the vector GZðFðxÞ;BðxÞÞ
T
,(2.2) takes the form
vG
vt
Z c s
K
C
1
n
2
s
C
vG
vx
;ð2:3Þ
M.W.Janowicz and others
2928
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
where
s
K
Z
0 0
1 0
!
;s
C
Z
0 1
0 0
!
:
Together with s
K
and s
C
,we shall also use the matrices,
s
11
Z
0 0
0 1
!
;s
22
Z
1 0
0 0
!
;
as well as the standard Pauli matrices s
x
and s
z
,
s
x
Z
0 1
1 0
!
;s
z
Z
1 0
0 K1
!
:
Using the identity (e.g.Louisell 1973)
e
iðus
C
Cvs
K
Þ
Zcos
ﬃﬃﬃﬃﬃ
uv
p
Ci
sin
ﬃﬃﬃﬃﬃ
uv
p
ﬃﬃﬃﬃﬃ
uv
p
ðus
C
Cvs
K
Þ;
we obtain as the solution of (2.3)
Gðx;t CDtÞ Z coshða
0
v
x
Þ Csinhða
0
v
x
Þ ns
K
C
1
n
s
C
Gðx;tÞ;ð2:4Þ
where a
0
Za=nZcDt=n.For any realanalytic function f(x),we have
coshðav
x
Þf ðxÞ Z
1
2
ðe
av
x
Ce
Kav
x
Þf ðxÞ Z
1
2
ðf ðx CaÞ Cf ðxKaÞÞ;
as well as
sinhðav
x
Þf ðxÞ Z
1
2
ðe
av
x
Ke
Kav
x
Þf ðxÞ Z
1
2
ðf ðx CaÞKf ðxKaÞÞ;
so that if we write xZja
0
,tZlDt with integer j,l,we obtain the following
difference scheme to integrate Maxwell’s equations:
Fðj;l C1Þ
Bðj;l C1Þ
!
Z
1
2
Fðj C1;lÞCFðjK1;lÞC
1
n
ðBðj C1;lÞKBðjK1;lÞÞ
nðFðj C1;lÞKFðjK1;lÞÞCBðj C1;lÞCBðjK1;lÞ
0
B
@
1
C
A
:ð2:5Þ
It is clear that equation (2.5) is an exact consequence of the Maxwell equations.
We now perform a check of validity of the above simple difference scheme,
namely,let us consider the initialvalue problem for the ﬁelds propagating in a
vacuum (nZ1):
Fðx;0Þ ZF
0
cosðkxÞ;Bðx;0Þ ZKF
0
cosðkxÞ;ð2:6Þ
for 0!x!l and Fðx;0ÞZBðx;0ÞZ0 for all other x.Contrary to appearances,this
problem is not entirely trivial,since it involves discontinuities at two spatial
points while we have derived (2.5) under the assumption of real analyticity of the
ﬁelds.We would like to check numerically whether the algorithm is able to treat
the discontinuities in the Cauchy data properly and that the rounding errors do
2929
Cellular automaton approach
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
not accumulate in a signiﬁcant way.The exact solutions are,naturally,
Fðx;tÞ ZF
0
cosðkðctKxÞÞ½QðtKðxKlÞ=cÞKQðtKx=cÞ;
Bðx;tÞ ZKF
0
cosðkðctKxÞÞ½QðtKðxKlÞ=cÞKQðtKx=cÞ:
The results of comparison of the above exact solution with the numerical one are
presented in ﬁgure 1 for F
0
Z100,kZ0.1.In this ﬁgure,the relative error of our
algorithm with respect to the above exact solution has been plotted with respect
to the coordinate x after 200 000 timesteps.
The two solutions are completely indistinguishable even for large numbers of
timesteps—the relative error has not exceeded 10
K7
.Thus,there are no
difﬁculties with discontinuities in the Cauchy data.We have also observed that
the rounding errors do not accumulate.
If we want to write down a similar difference scheme for dielectrics with an
xdependent refractive index,we have to sacriﬁce the exact nature of the
automaton and introduce approximations.
Let n(x) be a positive (the case of lefthanded materials is not covered here)
refractive index varying in space,let
^
A denote the operator v
x
and let
^
B denote
ð1=n
2
ðxÞÞv
x
.We have the following formal solution for the vector G:
Gðx;t CDtÞ Zexp½cDtð
^
As
K
C
^
Bs
C
ÞGðx;tÞ:ð2:7Þ
We do not know of any simple way to represent explicitly the action of an
operator of the type
expðconst:!
^
A
^
BÞ;
on a function of x.Operators of this form appear explicitly when we evaluate the
expression on the righthand side of equation (2.7).Therefore,we must apply an
approximation to this to obtain a difference scheme.The simplest one relies on
0
1×10
–8
2×10
–8
3×10
–8
4×10
–8
5×10
–8
6×10
–8
7×10
–8
199800
200000
200200
200400
200600
200800
201000
201200
relative error
x/a
Figure 1.Relative error in the numerical solution shown for the propagation of an initially
discontinuous signal (2.6) after 200 000 timesteps.
M.W.Janowicz and others
2930
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
the requirement that n(x) is a slowly varying function of x,so that we can ignore
the xdependence of n(x) when we evaluate the exponential function in (2.7).
Then we end up with a scheme (2.5),where the refractive index n is not constant
but is taken at the point x corresponding to the cell j.
In order to estimate the error and discover the limitations of the above
procedure,we calculate the expression RZGðx;tCDtÞKGðx;tÞ as well as a
similar expression R
neg
obtained from R by neglecting the derivatives of n(x)
over x,both in the lowest order with respect to Dt.We ﬁnd
R Z cDt v
x
s
K
C
1
n
2
ðxÞ
s
C
CðcDtÞ
2
v
x
1
n
2
ðxÞ
v
x
s
11
C
1
n
2
ðxÞ
v
2
x
s
22
COððDtÞ
3
Þ
Gðx;tÞ
and
R
neg
Z cDt v
x
s
K
C
1
n
2
ðxÞ
s
C
CðcDtÞ
2
1
n
2
ðxÞ
v
2
x
s
11
C
1
n
2
ðxÞ
v
2
x
s
22
COððDtÞ
3
Þ
Gðx;tÞ:
The difference RKR
neg
is equal to
ðcDtÞ
2
½K2ðn
0
ðxÞ=n
3
ðxÞÞs
11
vGðx;tÞ
vx
;
which gives conditions on the time increment Dt in terms of the variation of the
refractive index and the ﬁelds themselves.In particular,one might think that for
a stepwise change in the n characteristic for physically interesting systems like
photonic crystals,BBA is not applicable at all,because there are points where
n
0
(x) does not exist.However,we will provide numerical evidence to argue that
the situation is not that bad,at least for systems with piecewise constant
refractive indices.In that case,our CA propagation in each particular layer is
still exact.In addition,we have observed that the ﬁeld values at neighbouring
cells on both sides of any boundary are reasonably close to each other for small
and moderate values of the difference of refractive indices.This suggests that the
continuity conditions are approximately fulﬁlled at the boundary.More
importantly,the amplitudes of the reﬂected waves and of the transmitted
waves agree remarkably well with the Fresnel expressions for those amplitudes if
the difference of indices is smaller than about 5.In table 1,we have shown the
relative errors for the amplitude of reﬂected and transmitted waves when an
initial sinusoidal pulse falls on a single boundary between a vacuum and
transparent dielectric with refractive index n.We have performed simulations
with the wavelength of the pulse equal to 200 cells,i.e.lZ200a,and with the
length of the pulse equal to 10 wavelengths.
The data shown in table 1 are the ratios
A
rr
Z A
r
K
nK1
n C1
nK1
n C1
2931
Cellular automaton approach
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
for the reﬂected wave and
A
tr
Z A
t
K
2
n C1
2
n C1
for the transmitted wave,where A
r
and A
t
are the amplitudes of the reﬂected and
transmitted wave read from the numerical data provided by a F
ORTRAN
77
implementation of our algorithm (see electronic supplementary material).
We note that the error in the ratios of the reﬂected and transmitted amplitude
coefﬁcients behaves even better,e.g.the relative error of that ratio for nZ4 is as
small as 1.17!10
K10
.Nevertheless,we have observed that for n larger than
about 5 the algorithm breaks down,producing very big numbers for nz6.Thus,
in any dispersionless layered system with strong differences between the
refractive indices of the neighbouring layers,our cellular automaton is unreliable.
The simplest and most natural remedy would be to introduce thin ‘buffer’ layers
at the boundaries between two dielectrics with indices n
1
and n
2
in such a way
that the refractive index in the ‘buffer’ varies relatively slowly and interpolates
between the values n
1
and n
2
.It should be noted that the ‘buffer’ must be as thin
as possible in order to avoid the distortion of the wings of the reﬂected and
transmitted waves;moreover,the thinner the buffer,the less the relative error in
the amplitudes of these waves will be.We have found that with the buffer present
even for the dielectric constant of the medium as large as 900,the error does not
exceed 1%.Therefore,one might conclude that if n is smaller than 4,it is
advisable to retain the algorithmin its original form,while for larger n a buffer is
necessary.
The above considerations pertain to the least favourable case fromthe point of
view of our cellular automaton,namely,that which involves a discontinuous
initialvalue function and a discontinuous parameter n(x).However,if the initial
value is a smooth function,the situation is better.Let us again consider a
wavetrain falling from the vacuum side on the border between a vacuum and the
dielectric with the refractive index nZ5.The initial F and B functions are given
by (2.6),but modulated with a Gaussian function.Figure 2 illustrates the
comparison of the numerical and exact results for that case with ﬁgure 2a
showing almost indistinguishable plots of the ﬁeld F for the two cases,while
ﬁgure 2b displays the relative error that does not exceed 3.5%.
We now turn to the case of homogeneous conducting media for which we
propose the following BBA.Let s be the conductivity and let hZcm
0
s (where m
0
Table 1.Relative errors in the amplitudes of reﬂected and transmitted waves for various refractive
indices of the dielectric.
n A
rr
(10
K4
) A
tr
(10
K4
)
1.33 0.40 0.38
1.5 0.73 0.73
2.0 1.71 1.71
3.0 2.47 2.47
4.0 4.16 4.16
5.0 4.81 4.81
M.W.Janowicz and others
2932
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
is the magnetic permeability of the vacuum).Then Maxwell’s equations in (1C1)
dimensions can be written as
vG
vt
Zc v
x
s
x
K
h
2
ð1 Cs
z
Þ
h i
Gðx;tÞ;ð2:8Þ
with the meaning of G the same as before.We can thus write the formal
solution
Gðx;t CDtÞ Ze
Kð1=2ÞhcDt
exp cDt v
x
s
x
K
1
2
hs
z
h i
Gðx;tÞ:ð2:9Þ
The difﬁculty in obtaining the exact difference scheme follows from the
difﬁculty in evaluating the expression
exp
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
h
2
=4 Cv
2
x
q
gðxÞ;
–400
–300
–200
–100
0
100
200
300
400
5500
5550
5600
5650 5700
E/c
x/a
numerical results
exact results
0
0.005
0.010
0.015
0.020
0.025
0.030
0.035
relative error
l =1000a, k =0.1
l = 1000a, k = 0.1
(a)
(b)
Figure 2.(a) Comparison of numerical and exact solutions shown for the propagation of a
transmitted Gaussianmodulated wavetrain after crossing the boundary between the vacuum
(x!d) and the dielectric (xOd) with a refractive index of nZ5.The Gaussian envelope has a
width 50a and was centred at xZ500a.The boundary was located at dZ5000a.(b) Relative error
for solutions displayed in (a).
2933
Cellular automaton approach
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
where g(x) is any smooth function.It would be tempting to invoke the
Lie–Trotter formula here (see ch.VIII of Reed & Simon (1974)).In brief,the
Lie–Trotter theoremstates that if A,B and ACB are selfadjoint operators,then
sKlim
m/N
½expðiAðt=mÞÞexpðiBðt=mÞÞ
m
ZexpðitðACBÞÞ;ð2:10Þ
where sKlim denotes the limit in the strong (norm) sense.This means that if we
manage to write the evolution equation for G in the form
vG
vt
ZiðACBÞG;
with selfadjoint A,B and ACB,we can be sure that by taking the limit Dt/0
ðDtZt=mÞ,the approximate solution given by
GðtÞ Z½expðiADtÞexpðiBDtÞ
m
Gð0Þ;
i.e.that obtained by applying m times the operator expðiADtÞexpðiBDtÞ,
converges to the exact solution.This fact has been used as a basis for what is
called the ‘splitoperator’ or ‘splitting operator’ method (see Feit et al.1982).
This method is often used in the symmetric form,that is
e
iðACBÞDt
ze
ð1=2ÞiBDt
e
iADt
e
ð1=2ÞiBDt
;
and when combined with the fast Fourier transformation,it has recently been
extensively used in the context of the physics of atomic interactions with strong
light pulses (e.g.Pont et al.1991 or Dimou & Kull 2000),the physics of atomic
collisions (e.g.Schultz et al.2002),in molecular processes (Chu & Chu 2001) and
the theory of Bose condensates (Jackson & Zaremba 2002).
The difﬁculty which we immediately encounter,however,is that the operator
multiplying the vector G to produce the righthand side of equation (2.8) is not
selfadjoint (due to the presence of damping associated with the nonzero
conductivity).Therefore,we shall use the thesis of the Lie–Trotter theorem only
as a kind of heuristic device to write
Gðx;t CDtÞze
Kð1=2ÞhcDt
expðcDts
x
v
x
Þexp KcDt
h
2
s
z
Gðx;tÞ;ð2:11Þ
and postpone the proof of stability (which implies convergence in the considered
case) of the difference scheme obtained in this way to the ﬁrst part of appendix
A.If we compute the difference between the exact and approximate evolution
operators (i.e.the operators which multiply G on the righthand sides of
equations (2.9) and (2.11)),we obtain
KihðcDtÞ
2
v
x
s
y
COððDtÞ
3
Þ:
Accurate results may,therefore,be obtained only if Dt is much smaller than the
square root of 1=ðc
2
hjkjÞ,where k enumerates spatial harmonics of the Fourier
decomposition of G.That is,the propagation of harmonics with wavenumbers
greater that 1=ðc
2
ðDtÞ
2
hÞ will not be accurately covered by the algorithm.
M.W.Janowicz and others
2934
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
3.Propagation in Drudelike materials
In this section,we apply the cellular automaton method to the description of
propagation in (the layers of) homogeneous metals.The dielectric function of the
metal is modelled according to the Drude prescription,i.e.we have
eðuÞ Z1K
u
2
p
u
2
Cigu
:ð3:1Þ
The above dielectric function results from the following system of partial
differential equations for the scaled electric ﬁeld FðxÞZð0;0;FðxÞÞ,the magnetic
induction BZð0;BðxÞ;0Þ and an additional velocity ﬁeld VZð0;0;VðxÞÞ:
v
vt
Fðx;tÞ
Bðx;tÞ
Vðx;tÞ
0
B
@
1
C
AZ
cvBðx;tÞ=vxKcm
0
gðx;tÞVðx;tÞ
cvFðx;tÞ=vx
c
M
gðx;tÞFðx;tÞKgVðx;tÞ
0
B
B
B
B
@
1
C
C
C
C
A
;ð3:2Þ
where the constant M is equal to ðe
0
u
2
p
Þ
K1
.Indeed,by taking the Fourier
transformation of equation (3.2) with respect to time for a homogeneous time
independent medium (i.e.for gðx;tÞZ1),and by eliminating B and V,we obtain
a ‘onedimensional Helmholtz equation’:
v
2
F
vx
2
C
u
2
c
2
1K
m
0
c
2
MuðuCigÞ
F Z0;
so that the quantity
1K
m
0
c
2
MuðuCigÞ
ð3:3Þ
can be interpreted as a dielectric function of the medium (cf.Jackson (1982)).In
fact,the system of partial differential equations above was introduced precisely
in order to model a medium with that dielectric function.
If we now rescale the time and space variables and introduce the new velocity
ﬁeld W according to
t Zau
p
t;x Za
u
p
c
x;W Zðcm
0
=u
p
ÞV;
where the constant a speciﬁes the length of the timestep,we obtain the following
new system of equations,in which all the dependent variables have the same
dimension (of the magnetic induction) while the independent variables are
dimensionless:
v
vt
Fðx;tÞ
Bðx;tÞ
Wðx;tÞ
0
B
@
1
C
A
Z
vBðx;tÞ=vxKðgðx;tÞ=aÞWðx;tÞ
vFðx;tÞ=vx
ðgðx;tÞ=aÞFðx;tÞKðG=aÞWðx;tÞ
0
B
@
1
C
A
;ð3:4Þ
2935
Cellular automaton approach
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
where GZg=u
p
and where we have introduced the function gðx;tÞZgðx;tÞ to
specify the regions in space and time where the metal resides,i.e.where the
coupling between the velocity ﬁeld and the electromagnetic ﬁeld takes place.
To deal with the above system of differential equations,we propose the
following threestep scheme to obtain the triple (F,B,W) at time tCDt from
(F,B,W) at time t.
First,we propagate only the electromagnetic ﬁeld:
F
1
ðxÞ
B
1
ðxÞ
W
1
ðxÞ
0
B
@
1
C
AZ
1
2
½Fðx CDt;tÞ CFðxKDt;tÞ CBðx CDt;tÞKBðxKDt;tÞ
1
2
½Fðx CDt;tÞKFðxKDt;tÞ CBðx CDt;tÞ CBðxKDt;tÞ
Wðx;tÞ
0
B
B
B
B
@
1
C
C
C
C
A
:
ð3:5Þ
Second,we couple the electric ﬁeld with the velocity ﬁeld according to
F
2
ðxÞ
B
2
ðxÞ
W
2
ðxÞ
0
B
@
1
C
A
Z
F
1
ðxÞcosðgðx;tÞDt=aÞKW
1
ðxÞsinðgðx;tÞDt=aÞ
B
1
ðxÞ
F
1
ðxÞsinðgðx;tÞDt=aÞ CW
1
ðx;tÞcosðgðx;tÞDt=aÞ
0
B
@
1
C
A
:ð3:6Þ
The third step consists of multiplying W
2
(x) with the function
expðKgðx;tÞGDt=aÞ:
Fðx;t CDtÞ
Bðx;t CDtÞ
Wðx;t CDtÞ
0
B
@
1
C
A
Z
F
2
ðxÞ
B
2
ðxÞ
W
2
ðxÞexpðKgðx;tÞGDt=aÞ
0
B
@
1
C
A
:ð3:7Þ
The measure of the timestep in our scheme is given by a;this parameter must
be much larger than k,where again k is the wavenumber enumerating the spatial
Fourier components of G(harmonics with k larger than a will not be propagated
properly by our automaton).
It is interesting to compare the results obtained with the use of our CA with
those produced by other methods.As a benchmark we have chosen an algorithm
for the inversion of the Laplace transformation proposed by Hosono (see Hosono
1980) and improved in Wyns et al.(1989) (cf.Wyns et al.(1989) for discussion of
that algorithm).
An example of the comparison of the two methods above is shown in ﬁgure 3,
where the snapshots of solutions of the following initial boundaryvalue problem
have been displayed for tZ20 000.The space was assumed to be occupied by the
homogeneous (g(x)Z1) medium with the response function given by equation
(3.3) with gZ0.The ﬁelds at xZ0 are ﬁxed functions of t,tO0:
FðtÞZsin ðnt=aÞ,B(t)Z0.The other parameters of the system were chosen to
be:nZ1.2,aZ100.The solid line in ﬁgure 3 represents the results obtained with
the use of our automaton,while the dashed line represents those obtained from
the inversion of the Laplace transformation.
M.W.Janowicz and others
2936
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
So far we have only been able to prove the stability (which in our case
implied convergence) of the above scheme in the case of homogeneous media,
i.e.for gðx;tÞZ1.One may object that in view of the above difﬁculties in
describing layered systems with large n without a ‘buffer’ layer,a similar
buffer should be introduced in the case of plasmas.This may be the case,
although in all our numerical simulations we have observed excellent
behaviour of electromagnetic ﬁelds at the boundaries—the ﬁelds are as
‘continuous’ as possible for a system on a grid.In addition,there is a principal
difference between the plasma systems of this section and the perfectly
transparent dielectrics considered earlier.Whereas the response of those
transparent dielectrics to the incoming ﬁelds is ‘stiff’,instantaneous and static,
the Drudemedium can build its response in a dynamic way with the dielectric
function changing smoothly over space and time before it takes a well
established value in the bulk material,even though the region with nonzero
coupling between F and V has a welldeﬁned beginning and end.In the case of
necessity,one can easily create buffer layers by selecting groups of cells with
the values of g(x,t)Zg(x) interpolating between 0 and 1.The above comments
also apply to the case of the dispersive dielectrics discussed in §4.
4.Dispersive and lossy dielectrics
We now address the problemof dispersive dielectrics.We wish to obtain a simple
algorithm which would reduce to BBA in the absence of the dielectric.This can
naturally be achieved by adding two new ﬁelds representing the matter,namely,
the polarization ﬁeld XZ(0,0,X) and the associated velocity VZ(0,0,V).The
Maxwell equations together with the Newton equations for the polarization ﬁeld
–0.8
–0.6
–0.4
–0.2
0
0.2
0.4
0.6
0.8
0
5000
10000
15000
20000
E/c
x
t=20000
CA results
Laplacetransform
results
Figure 3.Comparison of results obtained with the use of CA with the results of the
Laplacetransform inversion,according to the algorithm of Wyns et al.(1989).
2937
Cellular automaton approach
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
in one spatial dimension read
v
vt
Fðx;tÞ
Bðx;tÞ
Xðx;tÞ
Vðx;tÞ
0
B
B
B
B
@
1
C
C
C
C
A
Z
cvBðx;tÞ=vxKcm
0
gðx;tÞVðx;tÞ
cvFðx;tÞ=vx
Vðx;tÞ
Ku
2
0
Xðx;tÞ C
c
M
gðx;tÞFðx;tÞKgVðx;tÞ
0
B
B
B
B
B
B
B
@
1
C
C
C
C
C
C
C
A
:ð4:1Þ
The function g(x,t) now speciﬁes the location of our dispersive dielectric.It is
equal to 1 in the spacetime regions where the dielectric is present and zero where
it is absent.The above model of the dielectric is usually associated with the name
of Lorentz (cf.Jackson (1982)).The constant Mis equal to ðbe
0
u
2
0
Þ
K1
,where b is
the static polarizability of the material while u
0
is the frequency of a single
resonance line (there is no problem in describing more complicated linear
dielectrics—the number of polarization ﬁelds has to be larger then).The
damping constant g characterizes the width of the resonance line.
The model considered here leads to the following dielectric function for the
medium:
eðuÞ Z1K
bu
2
0
u
2
Ku
2
0
Cigu
:ð4:2Þ
One can see this by performing the Fourier transformation of both sides of
equation (4.1) with respect to time (with g(x,t)Z1)),and eliminating all ﬁelds in
favour of F.The resulting equation is
v
2
F
vx
2
C
u
2
c
2
eðuÞF Z0;
where e(u) is given by (4.2),so that e(u) acquires its natural interpretation of the
dielectric constant or dielectric function.We immediately note that we can take
a limit:u
0
/0 with bu
2
0
Zu
2
p
kept constant to obtain a dielectric function
characteristic for the Drude metals considered in §3.In terms of equation (4.1)
this means that X and V remain constant in the absence of the electric ﬁeld.
Taking advantage of the existence of the natural timescale associated with u
0
,
we introduce new dimensionless time and coordinate variables:
t Zau
0
t;x Za
u
0
c
x;
where a is an arbitrary positive parameter to set the time and space scales.Thus,
our timestep is equal to T
0
=ð2paÞ,where T
0
is the period of free oscillations of
the polarization ﬁeld.Most of our trial numerical simulations have been obtained
for aZ100;however,we have checked that the results do not depend on a,at
least for 25!a!200.We now rescale the polarization variables in such a way
that all dependent variables have the dimensions of the magnetic induction:
X Z
1
cm
0
Y;V Z
u
0
cm
0
W:
M.W.Janowicz and others
2938
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
As a result,we obtain a system of equations simpler than (4.1):
v
vt
Fðx;tÞ
Bðx;tÞ
Yðx;tÞ
Wðx;tÞ
0
B
B
B
B
@
1
C
C
C
C
A
Z
vBðx;tÞ=vxKgðx;tÞWðx;tÞ=a
vFðx;tÞ=vx
Wðx;tÞ=a
KYðx;tÞ=aCbgðx;tÞFðx;tÞ=aKðG=aÞWðx;tÞ
0
B
B
B
B
@
1
C
C
C
C
A
;ð4:3Þ
where GZg=u
0
.
Owing to the favourable 4!4 matrix structure of equation (4.3),we can
propose the following scheme (‘automaton’) to integrate the above system.We
ﬁrst write (cf.(2.5) and (3.5))
F
1
ðxÞ
B
1
ðxÞ
Y
1
ðxÞ
W
1
ðxÞ
0
B
B
B
B
@
1
C
C
C
C
A
Z
1
2
½Fðx CDt;tÞ CFðxKDt;tÞ CBðx CDt;tÞKBðxKDt;tÞ
1
2
½Fðx CDt;tÞKFðxKDt;tÞ CBðx CDt;tÞ CBðxKDt;tÞ
Yðx;tÞcosðDt=aÞ CWðx;tÞsinðDt=aÞ
KYðx;tÞsinðDt=aÞ CWðx;tÞcosðDt=aÞ
0
B
B
B
B
B
B
B
B
B
B
@
1
C
C
C
C
C
C
C
C
C
C
A
;
ð4:4Þ
so that we obtain in the ﬁrst step the exact dynamics of the uncoupled
subsystems (except those for the damping of polarization ﬁelds).We believe that
this is quite an important feature;although the details of the dynamics of the
whole system are reproduced only approximately,the propagation of any
electromagnetic pulse from one spatial cell to another is fully causal.
In the second step,we take into account the interactions which result in
F
2
ðxÞ
B
2
ðxÞ
Y
2
ðxÞ
W
2
ðxÞ
0
B
B
B
B
@
1
C
C
C
C
A
Z
F
1
ðxÞcosð
ﬃﬃﬃ
b
p
gðx;tÞDt=aÞKW
1
ðxÞsinð
ﬃﬃﬃ
b
p
gðx;tÞDt=aÞ=
ﬃﬃﬃ
b
p
B
1
ðxÞ
Y
1
ðxÞ
F
1
ðxÞ
ﬃﬃﬃ
b
p
sinð
ﬃﬃﬃ
b
p
gðx;tÞDt=aÞ CW
1
ðx;tÞcosð
ﬃﬃﬃ
b
p
gðx;tÞDt=aÞ
0
B
B
B
B
B
@
1
C
C
C
C
C
A
:
ð4:5Þ
Equation (4.5) is valid,provided that bR0.However,if b!0,the trigonometric
functions should be replaced with hyperbolic functions and b with jbj with the
appropriate sign changes.This case corresponds to that of pumped media,which
will be discussed in more detail in a forthcoming paper.
To complete our algorithm,we take into account the damping of the
polarization ﬁeld,and hence we write
Fðx;t CDtÞ
Bðx;t CDtÞ
Yðx;t CDtÞ
Wðx;t CDtÞ
0
B
B
B
B
@
1
C
C
C
C
A
Z
F
2
ðxÞ
B
2
ðxÞ
Y
2
ðxÞ
W
2
ðxÞexpðKGDt=aÞ
0
B
B
B
B
@
1
C
C
C
C
A
:ð4:6Þ
2939
Cellular automaton approach
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
The above three steps are fully analogous to those of §3;now,however,there
exist nontrivial,oscillatory dynamics of the matter ﬁelds in the absence of
coupling with the electromagnetic ﬁeld.
The positive static polarizability b corresponds to the standard Lorentzian
dielectric.Dynamics of ﬁelds in such media have been the subject of several
fundamental papers by Oughstun and coworkers.To make our development
more speciﬁc and also more reliable,we provide examples of comparison of
results obtained by us with the use of our CA with those obtained using the
numerical procedure described in Wyns et al.(1989),as well as in comparison
with Oughstun’s asymptotic results.In this section,we follow Oughstun’s choice
of independent variable,which in virtually all his work is taken as qZct=xZt=x
with ﬁxed x (x).In Oughstun & Sherman (1989b),the following parameters
characterizing media have been chosen:
u
0
Z4!10
16
s
K1
;bu
2
0
Z20!10
32
s
K2
;gZ0:28!10
16
s
K1
;u
c
Z1!10
16
s
K1
;
where u
c
is the frequency of the applied ﬁeld.The above values correspond to a
highly dispersive and absorptive medium.We again consider the initial boundary
value problemas in §3 but for the Lorentzian mediumwith the above parameters
and for ﬁxed x
0
Z13 343,which corresponds to the ﬁxed value of x
0
Z10
K4
cm
found in Oughstun & Sherman (1989b).The boundary values at xZ0 are
Fð0;tÞZF
0
sinðntÞUðtÞ,Bð0;tÞZ0,where U(t) is the unitstep function,nZ
u
c
=ðau
0
Þ and F
0
Z2.This value of F
0
requires explanation,because the ‘second
canonical problem’ of Oughstun and coworkers involves the same formof F(0,t)
but with F
0
Z1.Oughstun and coworkers do not provide an explicit expression for
B(0,t),but it is given implicitly by their requirement that there is no propagation
for x!0 (x!0)—cf.a comment below eqn (9) in Wyns et al.(1989).The resulting
B(0,t) is fairly complicated,given by a Laplace transform which cannot be
inverted analytically.Therefore,in order to avoid the introduction of additional
errors due to the numerical inversion of the Laplace transformation,we have
chosen boundary values that give precisely the same Laplacetransformed electric
ﬁeld for xO0.That is,we have propagation for x!0,but for xO0 our results for the
electric ﬁeld are directly comparable with those of Oughstun.
In ﬁgure 4,we have shown the dynamics of the Sommerfeld forerunner
(Brillouin 1960) as obtained from our CA and from the numerical Laplace
transformation inverting algorithm developed in Wyns et al.(1989).
As can be clearly seen,excellent agreement has been found.Figure 4 can
be compared to both ﬁg.3 of Oughstun &Sherman (1989b) and ﬁg.5 of Wyns et al.
(1989) to see that very goodagreement with Oughstun’s and Sherman’s asymptotic
results has also been obtained for the Sommerfeld forerunner.
We now perform the same comparison for the Brillouin forerunner (Brillouin
1960) as well as for an initial part of the main signal just after the arrival of the
latter.The parameters characterizing the signal and the medium are the same as
before,but now x
0
Z10
K3
cm,which corresponds to x
0
Z133 426.The above
comparison is illustrated in ﬁgure 5.
It can once again be clearly seen that the agreement between the two numerical
methods is again very good,except for the amplitudes of the ﬁrst few peaks.It is
very difﬁcult to judge which method gives values closer to reality.This is
especially so because curiously enough the asymptotic results by Oughstun &
M.W.Janowicz and others
2940
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
Sherman—cf.ﬁg.10 of Oughstun & Sherman (1989b)—seem to give almost
exactly the arithmetic mean of the amplitudes of the main peak visible in our
ﬁgure 5.This interesting though not very important feature may be investigated
further in the future,but we still note the most obvious thing in ﬁgure 5:for the
arrival time of the Brillouin forerunner,the arrival time of the main signal and
the frequencies and phase relations,we obtain excellent agreement;the same is
true about the comparison of our ﬁgure 5 with ﬁg.10 of Oughstun & Sherman
(1989b).
5.Possible improvements and extensions
In this section,we would like to outline very brieﬂy how BBA can be extended or
improved in several directions.First,we observe that it can be employed to
–0.04
–0.02
0
0.02
0.04
0.06
0.08
0.10
0.12
1.45
1.50
1.55
1.60
1.65
1.70
1.75
1.80
E/c
x=133426
CA results
Laplacetransform results
q =t/x
Figure 5.Brillouin forerunner and a part of the main signal just after its arrival,as obtained from
our CA and by inverting the Laplace transformation using the algorithm by Wyns et al.(1989).
–0.004
–0.002
0
0.002
0.004
1.00
1.05
1.10
1.15
1.20
1.25
E/c
x=13343
CA results
Laplacetransform results
q=t/x
Figure 4.Timedependence of the Sommerfeld forerunner as obtained from our CA and by
inverting the Laplace transformation using the algorithm by Wyns et al.(1989).
2941
Cellular automaton approach
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
study (phenomenological models of) metamaterials with effective negative
refractive index (i.e.‘lefthanded materials’) in the spirit of Ziolkowski &
Heyman (2001).This is done quite simply by coupling the magnetic ﬁeld with the
magnetization ﬁeld (and associated velocity) in the same way as the electric ﬁeld
has been coupled here with polarization;then an appropriate splitting of the
evolution operator follows.In the case of homogeneous materials,the proof of
stability and convergence can again be provided without much difﬁculty.
Second,one can consider the extension to include the nonlinear Maxwell–
Bloch equations.Splitting of the timeevolution operator can be performed in
several ways,all of which are based on physical intuition.However,we have had
difﬁculty in proving the stability of our cellularautomaton difference schemes,
and this subject is still far from having any convincing treatment.
Third,following BiałynickiBirula (1994),we can extend CA to the case of
two and threedimensional propagation.For the sake of brevity,we consider two
spatial dimensions and are restricted to transverse electric (TE) polarization.
Namely,now let the ﬁeld F have only one nonzero component depending on x
and y,so that
F Zð0;0;F
z
ðx;yÞÞ;
while the magnetic induction has two components:
BZðB
x
ðx;yÞ;B
y
ðx;yÞ;0Þ:
Then the timedependent pair of Maxwell’s equations take the form
v
vt
F
z
B
x
B
y
0
B
@
1
C
A
Zc
vB
y
vx
K
vB
x
vy
!
K
vF
z
vy
vF
z
vx
0
B
B
B
B
B
B
B
B
B
@
1
C
C
C
C
C
C
C
C
C
A
;ð5:1Þ
which can be conveniently rewritten in matrix (or spinorial) form:
v
^
G
vt
Zcðs
x
v
x
Cs
y
v
y
Þ
^
G;ð5:2Þ
where
^
G Z
B
y
CiB
x
F
z
F
z
B
y
KiB
x
!
:ð5:3Þ
Equation (5.2) holds if and only if the condition of vanishing divergence of the B
ﬁeld is fulﬁlled.This equation forms the basis of a CA difference scheme on a
‘bodycentred’ grid in the xKy plane,obtained by using the Lie–Trotter theorem
to disentangle the exponential timeevolution operator that results from (5.2).
Interestingly,the difference scheme just outlined maintains the vanishing
divergence of B,provided that for every integration step G
12
remains real
(this is also a sufﬁcient condition for having G
12
ZG
21
at every step).The matrix
^
G can be augmented with the polarization ﬁeld and its associated velocity
M.W.Janowicz and others
2942
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
written as G
33
and G
44
;the linear coupling of G
44
with G
12
and G
21
provides us
with the integrator of twodimensional problems in dispersive dielectrics.
Finally,we note that our cellular automaton can be improved by using the
results of Suzuki (1990,1991).Motivated by applications to the QuantumMonte
Carlo simulations,Suzuki found approximations to the operator expðKxðACBÞÞ
in the form
e
t
1
A
e
t
2
B
e
t
3
A
e
t
4
B
.e
t
m
A
;
with explicitly constructed t
j
,such that the method is Oðx
mC1
Þ for small x.
Naturally,our problem ﬁts into Suzuki’s scheme well and his technique can be
applied to make the accuracy of our CA still higher.
6.Final remarks
In conclusion,we have proposed an extension of BiałynickiBirula’s algorithm
suitable for the numerical study of propagation in simple metals as well as
inhomogeneous and dispersive dielectrics.The algorithm is conceptually
extremely simple,since it only requires elementary properties of Pauli matrices
and it is very easy to implement.Indeed,a F
ORTRAN
77 code to implement BBA
contains just a few tens of lines,including declarations and input–output
instructions (see electronic supplementary material).Although from the point of
view of numerical analysis BBA forms just a speciﬁc explicit scheme to integrate
hyperbolic equations,it has distinct advantages in being an example of a
quantumcellular automaton and providing within another context a natural link
to the pathintegral representations of the Weyl and Dirac particles.BBA seems
to offer an alternative to the Laplacetransformation approach and other
methods of solving the wave propagation problems in dispersive materials,and
allows generalizations to multidimensional and nonlinear problems where the
Laplace transformation becomes inefﬁcient.
We have seen that the application of the algorithmto the case of one dielectric
interface leads to the amplitudes of reﬂected and refracted waves with small
values of the errors.In our future research,we plan to employ the algorithm
described above to the analysis of propagation in twodimensional media and in
systems with moving dielectric walls and layers.Work is already underway on
twodimensional as well as on nonlinear variants of BBA that are suitable for an
investigation of forerunners in Maxwell–Bloch systems.
M.W.J.gratefully acknowledges ﬁnancial support from the Alexander von Humboldt Foundation
and J.M.A.A.is supported by a Royal Commission for the Exhibition of 1851 Research Fellowship.
This work has also been supported in part (for A.O.) by the Polish State Committee for Scientiﬁc
Research,grant no.PBZKBN044/P03/2001.
Appendix A
The cellular automaton described in this paper can be thought of as just a certain
explicit difference scheme to numerically integrate a systemof partial differential
equations.This deﬁnition raises the very important question of whether that
scheme is actually stable.So far we are only able to provide the (rather simple)
2943
Cellular automaton approach
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
proofs for homogeneous media.We shall use the standard von Neumann &
Richtmyer analysis in terms of Fourier series as in Richtmyer (1957).Strictly
speaking,this method is applicable to the initialvalue problem with periodic
boundary conditions.However,we can embed the region of x or x for which we
obtain our numerical solution in a very large but ﬁnite domain L,so that it is
possible to impose periodic boundary conditions at the borders of L.This makes
it possible to use the expansion in terms of a Fourier series,
G
1
ðx;tÞ Z
X
k
e
ikx
~
G
1
ðk;tÞ;
where G
1
is a vector equal to (F,B)
T
in the case of propagation in a vacuum
or homogeneous dielectric or conducting medium;(F,B,W)
T
in the case of
propagation in Drudelike medium,or (F,B,Y,W)
T
in the case of propagation in
the singleresonance dielectric medium.
All cellular automata rules which have been deﬁned in the main body of the
paper can be rewritten for Fourier coefﬁcients in the form
~
G
1
ðk;t CDtÞ Z
^
H
~
Gðk;tÞ:ðA1Þ
The matrix
^
H which appears in the last equation is called the ampliﬁcation
matrix.
We start with a formulation of the necessary condition of stability (see
Richtmyer 1957 for details).Letting l
(i)
be the eigenvalues of
^
H,then the
inequality is
jl
ðiÞ
j%1 COðDtÞ;ðA2Þ
which must be fulﬁlled for every Dt,such that 0!Dt!T,for every k,and for all i
(where T is the total simulation time and the index i enumerates the eigenvalues)
constitutes the von Neumann necessary condition for stability.In the afore
mentioned book by Richtmyer,four sufﬁcient conditions for stability are
provided.We shall mainly use the second one,which can be formulated as
follows.If all the eigenvalues m
(i)
of the product PZ
^
H
^
H of the complex
conjugate of the ampliﬁcation matrix and ampliﬁcation matrix itself satisfy the
inequality
m
ðiÞ
%1 COðDtÞ;ðA3Þ
for all 0!Dt!T and all k,then the difference scheme is stable.In the case of
propagation in a vacuum or in the homogeneous dielectric,there is actually
nothing to prove:the scheme is exact,and as shown in ﬁgure 1,the rounding
errors seem to accumulate very slowly if at all.We therefore address the case of
propagation in homogeneous conducting media.
The matrix
^
H has the form
^
H Z
e
Kha
cosðkaÞ i sinðkaÞ
i e
Kha
sinðkaÞ cosðkaÞ
!
;ðA4Þ
where as always aZcDt.The above matrix has the eigenvalues
ð1=2Þe
Kha
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ð1 Ce
ha
Þ2 cos
2
ðkaÞK4 e
ha
q
Gð1 Ce
ha
ÞcosðkaÞ
:ðA5Þ
M.W.Janowicz and others
2944
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
The modulus of the eigenvalue with the ‘minus’ sign behind the square root is
bounded by either ð1=2Þð1CexpðKhaÞÞ or by expðKha=2Þ (both smaller than 1),
depending on the sign of the expression under the square root.The eigenvalue
with the ‘plus’ sign is bounded by either 1 or expðKha=2Þ,again depending on
whether the square root is real or imaginary.Thus,the moduli of both
eigenvalues are smaller than 1,so that von Neumann’s necessary condition for
stability is fulﬁlled.In order to prove the sufﬁcient condition,we multiply the
matrix
^
H by its complex conjugate and then expand eigenvalues of that product
in terms of Dt.Hence,we obtain the following pair:
ð1K2ha COða
2
Þ;1 COða
2
ÞÞ;
so that both eigenvalues are smaller or equal to 1CO(Dt),uniformly with respect
to k,so that the second sufﬁcient condition for stability is fulﬁlled.This concludes
the proof for the case of homogeneous conducting media.
We turn to the case of Drudelike media.Now,G
1
ZðF;B;WÞ
T
.In terms of
Fourier harmonics,equations (3.5)–(3.7) can again be written as
~
G
1
ðk;t CDtÞ Z
^
HG
1
ðk;tÞ;ðA6Þ
where the matrix
^
H is now the product of three matrices,
^
HZ
^
H
3
^
H
2
^
H
1
,where
^
H
3
Z
1 0 0
0 1 0
0 0 e
KGDt=a
0
B
B
@
1
C
C
A
;ðA7Þ
^
H
2
Z
cosðDt=aÞ 0 KsinðDt=aÞ
0 1 0
sinðDt=aÞ 0 cosðDt=aÞ
0
B
@
1
C
A
;ðA8Þ
^
H
1
Z
cosðkDtÞ i sinðkDtÞ 0
i sinðkDtÞ cosðkDtÞ 0
0 0 1
0
B
@
1
C
A
:ðA9Þ
The stability condition requires that there exists t
0
O0,such that for all k,all
Dt!t
0
and 0%mDt%T (T being the simulation time),the matrices
^
H
m
,are
uniformly bounded.We observe that the norm of
^
H is smaller or equal to the
product of the norms of matrices
^
H
i
,iZ1,2,3.However,the latter are normal
matrices;that is,they commute with their Hermitian conjugates.It follows that
their norms are equal to their spectral radii,i.e.their largest eigenvalues.The
moduli of the largest eigenvalues of
^
H
i
are exactly equal to 1,for iZ1,2,3.This
means that the norms of all the matrices
^
H
m
are also bounded by the constant 1,
that which was to be proved.
We now consider the scheme given by equations (4.4)–(4.6).The kth harmonic
of the vector GZðF;B;Y;WÞ
T
satisﬁes the equation
~
Gðk;tCDtÞ Z
^
H
3
^
H
2
^
H
1
~
Gðk;tÞ;ðA10Þ
2945
Cellular automaton approach
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
where
^
H
i
,iZ1,2,3 are 4!4 matrices:
^
H
3
Z
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 e
KGDt=a
0
B
B
B
B
@
1
C
C
C
C
A
;ðA11Þ
^
H
2
Z
cosð
ﬃﬃﬃ
b
p
Dt=aÞ 0 0 Ksinð
ﬃﬃﬃ
b
p
Dt=aÞ=
ﬃﬃﬃ
b
p
0 1 0 0
0 0 1 0
ﬃﬃﬃ
b
p
sinð
ﬃﬃﬃ
b
p
Dt=aÞ 0 0 cosð
ﬃﬃﬃ
b
p
Dt=aÞ
0
B
B
B
B
@
1
C
C
C
C
A
;ðA12Þ
^
H
1
Z
cosðkDtÞ i sinðkDtÞ 0 0
i sinðkDtÞ cosðkDtÞ 0 0
0 0 cosðDt=aÞ sinðDt=aÞ
0 0 KsinðDt=aÞ cosðDt=aÞ
0
B
B
B
B
@
1
C
C
C
C
A
:ðA13Þ
In order to estimate the norm of the matrix
^
H,we ﬁrst observe that
^
H
1
is
unitary.Thus,the norm of
^
H
m
is bounded by the norm of ð
^
H
3
^
H
2
Þ
m
.It is
therefore enough to check the necessary condition and one of the sufﬁcient
conditions for stability for the product
^
H
3
^
H
2
.We can observe that this product
has the following eigenvalues:1,1 and 2=ðAGBÞ,where
AZð1 Ce
GDt=a
Þcosð
ﬃﬃﬃ
b
p
Dt=aÞ;
B Z
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
K4 e
GDt=a
CA
2
p
:
The eigenvalues possess the following expansion with respect to Dt:
ð1;1;1Cð1=2aÞðKGK
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
G
2
K4b
p
ÞDt;1Cð1=2aÞðKGC
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
G
2
K4b
p
ÞDt COððDtÞ
2
ÞÞ;
uniform with respect to k,and for any mDt!T with integer m.Thus,the
necessary condition for stability is satisﬁed.Multiplication of
^
H
3
^
H
2
by its
complex conjugate (which is equal to the product itself),the computation of
eigenvalues and their expansion in terms of Dt lead to the result
ð1;1;1Kð1=aÞðGC
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
G
2
K4b
p
ÞDt;1Cð1=aÞðKGC
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
G
2
K4b
p
ÞDt COððDtÞ
2
ÞÞ;
so that our scheme satisﬁes again the second sufﬁcient condition for stability.
We should note here that in the apparently dubious (from the possibility of
exponential growth of the ﬁelds in time) case of media having a negative effective
static polarizability,we can still prove the stability of our CA—that proof,
however,will be presented elsewhere.
Finally,we observe that a theorem due to Lax asserts that for linear systems
with constant coefﬁcients stability implies convergence.This applies to all
M.W.Janowicz and others
2946
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
difference schemes introduced in this paper.Therefore,making a larger (so that
the timestep becomes smaller) in our system can only improve the results.
References
BiałynickiBirula,I.1994 Weyl,Dirac,and Maxwell equations on a lattice as unitary cellular
automata.Phys.Rev.D 49,6920.(doi:10.1103/PhysRevD.49.6920)
Brillouin,L.1960 Wave propagation and group velocity.New York:Academic Press.
Chopard,B.& Droz,M.1998 Cellular automata modelling of physical systems.Cambridge:
Cambridge University Press.
Chopard,B.,Luthi,P.O.& Wagen,J.F.1997 A lattice Boltzmann method for wave propagation
in urban microcells.IEE Proc.—Microwaves Antennas Propag.144,251.(doi:10.1049/ip
map:19971197)
Chu,X.& Chu,S.I.2001 Optimization of highharmonic generation by genetic algorithm and
wavelet timefrequency analysis of quantum dipole emission.Phys.Rev.A 64,063 404.(doi:10.
1103/PhysRevA.64.063404)
Dimou,L.& Kull,H.J.2000 Atomic and molecular processes in external ﬁelds—abovethreshold
ionization by strong anharmonic light pulses.Phys.Rev.A 61,043 404.(doi:10.1103/
PhysRevA.61.043404)
Dvorak,S.L.& Dudley,D.G.1995 Propagation of ultrawideband electromagnetic pulses through
dispersive media.IEEE Trans.Electromagn.Compat.37,192.(doi:10.1109/15.385883)
Dvorak,S.L.& Dudley,D.G.1996 A comment on propagation of ultrawideband electromagnetic
pulses through dispersive media—author’s reply.IEEE Trans.Electromagn.Compat.38,203.
Dvorak,S.L.,Ziolkowski,R.W.& Felsen,L.B.1998 Hybrid analytical–numerical approach for
modeling transient wave propagation in Lorentz media.J.Opt.Soc.Am.A 15,1241.
Feit,M.J.,Fleck,J.A.& Steiger,A.1982 Solution of the Schro¨dinger equation by a spectral
method.J.Comput.Phys.47,412.(doi:10.1016/00219991(82)900912)
Frisch,L.,Hasslacher,B.& Pomeau,Y.1986 Latticegas automata for the Navier–Stokes
equations.Phys.Rev.Lett.56,1505.(doi:10.1103/PhysRevLett.56.1505)
Hardy,J.,de Pazzis,O.& Pomeau,Y.1976 Molecular dynamics of a classical lattice gas:transport
properties and time correlation functions.Phys.Rev.A 13,1949.(doi:10.1103/PhysRevA.13.
1949)
Hosono,T.1980 Proc.1980 Int.Union of Radio Science Int.Symp.on Electromagnetic Waves.
Munich:International Union of Radio Science,FRG 1980 paper 112,pp.C1–C4.
Ilachinski,A.2001 Cellular automata.A discrete universe.Singapore:World Scientiﬁc.
Jackson,J.D.1982 Classical electrodynamics.Warsaw:Pan´stwowe Wydawnictwo Naukowe.
(Polish translation).
Jackson,B.& Zaremba,E.2002 Modelling Bose–Einstein condensed gases at ﬁnite temperatures
with Nbody simulations.Phys.Rev.A 66,033 606.(doi:10.1103/PhysRevA.66.033606)
Louisell,W.H.1973 Quantum statistical properties of radiation.New York:Wiley.
Oughstun,K.E.& Balictsis,C.M.1996 Gaussian pulse propagation in a dispersive,absorbing
dielectric.Phys.Rev.Lett.77,2210.(doi:10.1103/PhysRevLett.77.2210)
Oughstun,K.E.& Sherman,G.C.1989a Propagation of electromagnetic pulses in a linear
dispersive medium with absorption (the Lorentz medium).J.Opt.Soc.Am.B 5,817.
Oughstun,K.E.& Sherman,G.C.1989b Uniform asymptotic description of electromagnetic pulse
propagation in a linear dispersive medium with absorption (the Lorentz medium).J.Opt.Soc.
Am.A 6,1394.
Oughstun,K.E.& Sherman,G.C.1990 Uniform asymptotic description of ultrashort rectangular
optical pulse propagation in a linear,causally dispersive medium.Phys.Rev.A 41,6090.
(doi:10.1103/PhysRevA.41.6090)
2947
Cellular automaton approach
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
Pont,M.,Proulx,D.& Shakeshaft,R.1991 Numerical integration of the timedependent
Schro¨dinger equation for an atom in a radiation ﬁeld.Phys.Rev.A 44,4486.(doi:10.1103/
PhysRevA.44.4486)
Reed,M.& Simon,B.1974 Methods of modern mathematical physics,vol.I.New York:Academic
Press.
Richtmyer,R.D.1957 Difference methods for initialvalue problems.London:Interscience.
Rothman,D.H.& Zaleski,S.1994 Latticegas models of phase separation:interfaces,phase
transitions and multiphase ﬂow.Rev.Mod.Phys.66,1417.(doi:10.1103/RevModPhys.66.1417)
Schultz,D.R.,Reinhold,C.O.,Krstic,P.S.& Strayer,M.R.2002 Ejected electron spectrum in
lowenergy proton–hydrogen collisions.Phys.Rev.A 65,052 722.(doi:10.1103/PhysRevA.65.
052722)
Sebbah,P.& Vanneste,C.2002 Random laser in the localized regime.Phys.Rev.B 66,144 202.
(doi:10.1103/PhysRevB.66.144202)
Simons,N.R.S.,Bridges,G.E.& Cuhaci,M.1999 A lattice gas automaton capable of modeling
threedimensional electromagnetic ﬁelds.J.Comput.Phys.151,816.(doi:10.1006/jcph.1999.
6221)
Sornette,D.,Legrand,O.,Mortessagne,F.,Sebbah,P.& Vanneste,C.1993 The wave automation
for the timedependent Schro¨dinger,classical wave and Klein–Gordon equations.Phys.Lett.A
178,292.(doi:10.1016/03759601(93)91105E)
Suzuki,M.1990 Fractal decomposition of exponential operators with applications to manybody
theories and Monte Carlo simulations.Phys.Lett.A 146,319.(doi:10.1016/0375
9601(90)90962N)
Suzuki,M.1991 General theory of fractal path integrals with applications to manybody theories
and statistical physics.J.Math.Phys.32,400.(doi:10.1063/1.529425)
Vanneste,C.& Sebbah,P.2001 Selective excitation of localized modes in active random media.
Phys.Rev.Lett.87,183 903.(doi:10.1103/PhysRevLett.87.183903)
Wyns,P.,Foty,D.P.& Oughstun,K.E.1989 Numerical analysis of the precursor ﬁelds in linear
dispersive pulse propagation.J.Opt.Soc.Am.A 6,1421.
Ziolkowski,R.W.& Heyman,E.2001 Wave propagation in media having negative permittivity
and permeability.Phys.Rev.E 64,056625.(doi:10.1103/PhysRevE.64.056625)
M.W.Janowicz and others
2948
Proc.R.Soc.A (2006)
on November 30, 2013rspa.royalsocietypublishing.orgDownloaded from
Enter the password to open this PDF file:
File name:

File size:

Title:

Author:

Subject:

Keywords:

Creation Date:

Modification Date:

Creator:

PDF Producer:

PDF Version:

Page Count:

Preparing document for printing…
0%
Comments 0
Log in to post a comment