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 ﬁ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,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 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 ﬁnite-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 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 S-matrix formulation of the wave propagation

was used,which turned out to be equivalent to a ﬁnite-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 ﬁrst 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 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 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 ﬁ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 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

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

x-dependent 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 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 ﬁ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

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 ﬁ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 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 difﬁculty 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 ﬁ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 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 ﬁ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 ‘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

ﬁ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 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 ﬁeld and the electromagnetic ﬁeld 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 ﬁ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 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 ﬁgure 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 ﬁ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 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-deﬁ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

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 speciﬁes 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 ﬁ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 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 ﬁ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 non-trivial,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 co-workers.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 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

ﬁ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

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 ﬁ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 time-evolution 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 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 ﬁeld 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

ﬁeld is fulﬁlled.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 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 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 ﬁ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ł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 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 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 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 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 ﬁ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.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 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 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 ﬁ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 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 deﬁned in the main body of the

paper can be re-written 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 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

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 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 ﬁelds—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 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 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 ﬁ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 initial-value problems.London:Interscience.

Rothman,D.H.& Zaleski,S.1994 Lattice-gas 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

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 ﬁ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 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 ﬁ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

## Σχόλια 0

Συνδεθείτε για να κοινοποιήσετε σχόλιο