wave propagation in dispersive media Cellular automaton approach to electromagnetic

overwhelmedblueearthAI and Robotics

Dec 1, 2013 (3 years and 8 months ago)

198 views

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#ref-list-1
http://rspa.royalsocietypublishing.org/content/462/2074/2927.full.

This article cites 26 articles
Email alerting service
herethe box at the top right-hand 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,Carl-von-Ossietzky 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,02-668 Warsaw,Poland
Extensions of Białynicki-Birula’s cellular automaton are proposed for studies of the one-
dimensional propagation of electromagnetic fields in Drude metals,as well as in both
transparent,dispersive and lossy dielectrics.These extensions are obtained by
representing the dielectrics with appropriate matter fields,such as polarization together
with associated velocity fields.To obtain the different schemes for the integration of the
resulting systems of linear partial differential equations,split-operator ideas are
employed.Possible further extensions to two-dimensional propagation and for the
study of left-handed 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 so-called lattice-gas 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,02-668
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 efficiency of more traditional methods,such
as finite differences or finite 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 finite-difference
time-domain 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 CA-based integrators of the wave equations,the
authors decided to relax one of the most important defining 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 S-matrix formulation of the wave propagation
was used,which turned out to be equivalent to a finite-difference 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łynicki-Birula (1994),who has
been mostly interested in problems of a theoretical nature.One of our aims is to
demonstrate that in the case of one-dimensional propagation at least,the
Białynicki-Birula 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łynicki-Birula’s cellular automaton in (1D1) dimensions
Let us start with the first time-dependent 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 sufficient to represent
the electric field vector as EZð0;0;EðxÞÞ and the magnetic field 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
ffiffiffiffiffi
uv
p
Ci
sin
ffiffiffiffiffi
uv
p
ffiffiffiffiffi
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 real-analytic 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 initial-value problem for the fields 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
fields.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 significant 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 figure 1 for F
0
Z100,kZ0.1.In this figure,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 time-steps.
The two solutions are completely indistinguishable even for large numbers of
time-steps—the relative error has not exceeded 10
K7
.Thus,there are no
difficulties 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
x-dependent refractive index,we have to sacrifice the exact nature of the
automaton and introduce approximations.
Let n(x) be a positive (the case of left-handed 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 right-hand 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 time-steps.
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 x-dependence 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 find
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 fields 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 field 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 fulfilled at the boundary.More
importantly,the amplitudes of the reflected 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 reflected 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 reflected 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 reflected 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 reflected and transmitted amplitude
coefficients 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 reflected 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
initial-value 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 figure 2a
showing almost indistinguishable plots of the field F for the two cases,while
figure 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 reflected 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 difficulty in obtaining the exact difference scheme follows from the
difficulty in evaluating the expression
exp
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
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 Gaussian-modulated 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 self-adjoint 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 self-adjoint 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 ‘split-operator’ 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 difficulty which we immediately encounter,however,is that the operator
multiplying the vector G to produce the right-hand side of equation (2.8) is not
self-adjoint (due to the presence of damping associated with the non-zero
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 first 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 right-hand 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 Drude-like 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 field FðxÞZð0;0;FðxÞÞ,the magnetic
induction BZð0;BðxÞ;0Þ and an additional velocity field 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 ‘one-dimensional 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
field W according to
t Zau
p
t;x Za
u
p
c
x;W Zðcm
0
=u
p
ÞV;
where the constant a specifies the length of the time-step,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 field and the electromagnetic field takes place.
To deal with the above system of differential equations,we propose the
following three-step scheme to obtain the triple (F,B,W) at time tCDt from
(F,B,W) at time t.
First,we propagate only the electromagnetic field:
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 field with the velocity field 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 time-step 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 figure 3,
where the snapshots of solutions of the following initial boundary-value 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 fields at xZ0 are fixed 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 figure 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 difficulties 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 fields at the boundaries—the fields 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 fields is ‘stiff’,instantaneous and static,
the Drude-medium 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 non-zero
coupling between F and V has a well-defined 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 fields representing the matter,namely,
the polarization field XZ(0,0,X) and the associated velocity VZ(0,0,V).The
Maxwell equations together with the Newton equations for the polarization field
–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
Laplace-transform
results
Figure 3.Comparison of results obtained with the use of CA with the results of the
Laplace-transform 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 specifies the location of our dispersive dielectric.It is
equal to 1 in the space-time 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 fields 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 fields 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 field.
Taking advantage of the existence of the natural time-scale 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 time-step is equal to T
0
=ð2paÞ,where T
0
is the period of free oscillations of
the polarization field.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
first 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 first step the exact dynamics of the uncoupled
subsystems (except those for the damping of polarization fields).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ð
ffiffiffi
b
p
gðx;tÞDt=aÞKW
1
ðxÞsinð
ffiffiffi
b
p
gðx;tÞDt=aÞ=
ffiffiffi
b
p
B
1
ðxÞ
Y
1
ðxÞ
F
1
ðxÞ
ffiffiffi
b
p
sinð
ffiffiffi
b
p
gðx;tÞDt=aÞ CW
1
ðx;tÞcosð
ffiffiffi
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 field,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 non-trivial,oscillatory dynamics of the matter fields in the absence of
coupling with the electromagnetic field.
The positive static polarizability b corresponds to the standard Lorentzian
dielectric.Dynamics of fields in such media have been the subject of several
fundamental papers by Oughstun and co-workers.To make our development
more specific 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 fixed 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 field.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 fixed x
0
Z13 343,which corresponds to the fixed 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 unit-step 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 co-workers involves the same formof F(0,t)
but with F
0
Z1.Oughstun and co-workers 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 Laplace-transformed electric
field for xO0.That is,we have propagation for x!0,but for xO0 our results for the
electric field are directly comparable with those of Oughstun.
In figure 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 fig.3 of Oughstun &Sherman (1989b) and fig.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 figure 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 first few peaks.It is
very difficult 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.fig.10 of Oughstun & Sherman (1989b)—seem to give almost
exactly the arithmetic mean of the amplitudes of the main peak visible in our
figure 5.This interesting though not very important feature may be investigated
further in the future,but we still note the most obvious thing in figure 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 figure 5 with fig.10 of Oughstun & Sherman
(1989b).
5.Possible improvements and extensions
In this section,we would like to outline very briefly 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
Laplace-transform 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
Laplace-transform results
q=t/x
Figure 4.Time-dependence 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.‘left-handed materials’) in the spirit of Ziolkowski &
Heyman (2001).This is done quite simply by coupling the magnetic field with the
magnetization field (and associated velocity) in the same way as the electric field
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 difficulty.
Second,one can consider the extension to include the nonlinear Maxwell–
Bloch equations.Splitting of the time-evolution operator can be performed in
several ways,all of which are based on physical intuition.However,we have had
difficulty in proving the stability of our cellular-automaton difference schemes,
and this subject is still far from having any convincing treatment.
Third,following Białynicki-Birula (1994),we can extend CA to the case of
two- and three-dimensional propagation.For the sake of brevity,we consider two
spatial dimensions and are restricted to transverse electric (TE) polarization.
Namely,now let the field F have only one non-zero 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 time-dependent 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
field is fulfilled.This equation forms the basis of a CA difference scheme on a
‘body-centred’ grid in the xKy plane,obtained by using the Lie–Trotter theorem
to disentangle the exponential time-evolution 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 sufficient condition for having G
12
ZG
21
at every step).The matrix
^
G can be augmented with the polarization field 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 two-dimensional 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 fits 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łynicki-Birula’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 specific 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 path-integral representations of the Weyl and Dirac particles.BBA seems
to offer an alternative to the Laplace-transformation approach and other
methods of solving the wave propagation problems in dispersive materials,and
allows generalizations to multi-dimensional and nonlinear problems where the
Laplace transformation becomes inefficient.
We have seen that the application of the algorithmto the case of one dielectric
interface leads to the amplitudes of reflected 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 two-dimensional media and in
systems with moving dielectric walls and layers.Work is already underway on
two-dimensional 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 financial 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 Scientific
Research,grant no.PBZ-KBN-044/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 definition 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 initial-value 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 finite 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 Drude-like medium,or (F,B,Y,W)
T
in the case of propagation in
the single-resonance dielectric medium.
All cellular automata rules which have been defined in the main body of the
paper can be re-written for Fourier coefficients 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 amplification
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 fulfilled 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 sufficient 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 amplification matrix and amplification 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 figure 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
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ð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 fulfilled.In order to prove the sufficient 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 sufficient condition for stability is fulfilled.This concludes
the proof for the case of homogeneous conducting media.
We turn to the case of Drude-like 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
satisfies 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ð
ffiffiffi
b
p
Dt=aÞ 0 0 Ksinð
ffiffiffi
b
p
Dt=aÞ=
ffiffiffi
b
p
0 1 0 0
0 0 1 0
ffiffiffi
b
p
sinð
ffiffiffi
b
p
Dt=aÞ 0 0 cosð
ffiffiffi
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 first 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 sufficient
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ð
ffiffiffi
b
p
Dt=aÞ;
B Z
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
K4 e
GDt=a
CA
2
p
:
The eigenvalues possess the following expansion with respect to Dt:
ð1;1;1Cð1=2aÞðKGK
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
G
2
K4b
p
ÞDt;1Cð1=2aÞðKGC
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
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 satisfied.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
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
G
2
K4b
p
ÞDt;1Cð1=aÞðKGC
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
G
2
K4b
p
ÞDt COððDtÞ
2
ÞÞ;
so that our scheme satisfies again the second sufficient condition for stability.
We should note here that in the apparently dubious (from the possibility of
exponential growth of the fields 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 coefficients 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 time-step becomes smaller) in our system can only improve the results.
References
Białynicki-Birula,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 high-harmonic generation by genetic algorithm and
wavelet time-frequency 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 fields—above-threshold
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 ultra-wideband 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 ultra-wideband 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/0021-9991(82)90091-2)
Frisch,L.,Hasslacher,B.& Pomeau,Y.1986 Lattice-gas 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 Scientific.
Jackson,J.D.1982 Classical electrodynamics.Warsaw:Pan´stwowe Wydawnictwo Naukowe.
(Polish translation).
Jackson,B.& Zaremba,E.2002 Modelling Bose–Einstein condensed gases at finite temperatures
with N-body 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 time-dependent
Schro¨dinger equation for an atom in a radiation field.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 initial-value problems.London:Interscience.
Rothman,D.H.& Zaleski,S.1994 Lattice-gas models of phase separation:interfaces,phase
transitions and multiphase flow.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
low-energy 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
three-dimensional electromagnetic fields.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 time-dependent Schro¨dinger,classical wave and Klein–Gordon equations.Phys.Lett.A
178,292.(doi:10.1016/0375-9601(93)91105-E)
Suzuki,M.1990 Fractal decomposition of exponential operators with applications to many-body
theories and Monte Carlo simulations.Phys.Lett.A 146,319.(doi:10.1016/0375-
9601(90)90962-N)
Suzuki,M.1991 General theory of fractal path integrals with applications to many-body 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 fields 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