## Abstract

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.

## 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 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 1989*a*,*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 (1+1) 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):(2.1)where *ϵ* and *μ* 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 and the magnetic field vector as (the other polarization may be considered in an analogous way). It is convenient to introduce a new variable ** F**, such that

**=**

*E**c*

**.**

*F*Then Maxwell's equations take the form(2.2)where . If we introduce the vector , (2.2) takes the form(2.3)whereTogether with *σ*_{−} and *σ*_{+}, we shall also use the matrices,as well as the standard Pauli matrices *σ*_{x} and *σ*_{z},Using the identity (e.g. Louisell 1973)we obtain as the solution of (2.3)(2.4)where . For any real-analytic function *f*(*x*), we haveas well asso that if we write *x*=*ja*′, with integer *j*, *l*, we obtain the following difference scheme to integrate Maxwell's equations:(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 (*n*=1):(2.6)for 0<*x*<*l* and 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 not accumulate in a significant way. The exact solutions are, naturally,The results of comparison of the above exact solution with the numerical one are presented in figure 1 for *F*_{0}=100, *κ*=0.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^{−7}. 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 denote the operator ∂_{x} and let denote . We have the following formal solution for the vector ** G**:(2.7)

We do not know of any simple way to represent explicitly the action of an operator of the typeon 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 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 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 Δ

*t*. We findandThe difference

**−**

*R*

*R*_{neg}is equal towhich gives conditions on the time increment Δ

*t*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*′(

*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.

*λ*=200

*a*, and with the length of the pulse equal to 10 wavelengths.

The data shown in table 1 are the ratiosfor the reflected wave andfor 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 Fortran 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 *n*=4 is as small as 1.17×10^{−10}. Nevertheless, we have observed that for *n* larger than about 5 the algorithm breaks down, producing very big numbers for *n*≈6. 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 algorithm in its original form, while for larger *n* a buffer is necessary.

The above considerations pertain to the least favourable case from the 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 *n*=5. 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 2*a* showing almost indistinguishable plots of the field *F* for the two cases, while figure 2*b* 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 *σ* be the conductivity and let (where *μ*_{0} is the magnetic permeability of the vacuum). Then Maxwell's equations in (1+1) dimensions can be written as(2.8)with the meaning of ** G** the same as before. We can thus write the formal solution(2.9)

The difficulty in obtaining the exact difference scheme follows from the difficulty in evaluating the expressionwhere *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 theorem states that if *A*, *B* and *A*+*B* are self-adjoint operators, then(2.10)where *s*−lim denotes the limit in the strong (norm) sense. This means that if we manage to write the evolution equation for ** G** in the formwith self-adjoint

*A*,

*B*and

*A*+

*B*, we can be sure that by taking the limit , the approximate solution given byi.e. that obtained by applying

*m*times the operator , 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 isand 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(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

**on the right-hand sides of equations (2.9) and (2.11)), we obtainAccurate results may, therefore, be obtained only if Δ**

*G**t*is much smaller than the square root of , where

*k*enumerates spatial harmonics of the Fourier decomposition of

**. That is, the propagation of harmonics with wavenumbers greater that will not be accurately covered by the algorithm.**

*G*## 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(3.1)The above dielectric function results from the following system of partial differential equations for the scaled electric field , the magnetic induction and an additional velocity field :(3.2)where the constant is equal to . Indeed, by taking the Fourier transformation of equation (3.2) with respect to time for a homogeneous time-independent medium (i.e. for ), and by eliminating *B* and *V*, we obtain a ‘one-dimensional Helmholtz equation’:so that the quantity(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 towhere the constant *α* 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:(3.4)where and where we have introduced the function 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 from (*F*, *B*, *W*) at time *τ*.

First, we propagate only the electromagnetic field:(3.5)Second, we couple the electric field with the velocity field according to(3.6)The third step consists of multiplying *W*_{2}(*ξ*) with the function :(3.7)

The measure of the time-step in our scheme is given by *α*; 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

*α*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 *τ*=20 000. The space was assumed to be occupied by the homogeneous (*g*(*x*)=1) medium with the response function given by equation (3.3) with *γ*=0. The fields at *ξ*=0 are fixed functions of *τ*, *τ*>0: , *B*(*τ*)=0. The other parameters of the system were chosen to be: *ν*=1.2, *α*=100. 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.

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 . 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*)=*g*(*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 problem of 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 ** X**=(0, 0,

*X*) and the associated velocity

**=(0, 0,**

*V**V*). The Maxwell equations together with the Newton equations for the polarization field in one spatial dimension read(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 is equal to , where

*β*is the static polarizability of the material while

*ω*

_{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

*γ*characterizes the width of the resonance line.

The model considered here leads to the following dielectric function for the medium:(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*)=1)), and eliminating all fields in favour of *F*. The resulting equation iswhere *ϵ*(*ω*) is given by (4.2), so that *ϵ*(*ω*) acquires its natural interpretation of the dielectric constant or dielectric function. We immediately note that we can take a limit: with 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 *ω*_{0}, we introduce new dimensionless time and coordinate variables:where *α* is an arbitrary positive parameter to set the time and space scales. Thus, our time-step is equal to , where *T*_{0} is the period of free oscillations of the polarization field. Most of our trial numerical simulations have been obtained for *α*=100; however, we have checked that the results do not depend on *α*, at least for 25<*α*<200. We now rescale the polarization variables in such a way that all dependent variables have the dimensions of the magnetic induction:As a result, we obtain a system of equations simpler than (4.1):(4.3)where .

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))(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(4.5)Equation (4.5) is valid, provided that *β*≥0. However, if *β*<0, the trigonometric functions should be replaced with hyperbolic functions and *β* with |*β*| 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(4.6)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 *β* 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 with fixed *x* (*ξ*). In Oughstun & Sherman (1989*b*), the following parameters characterizing media have been chosen:where *ω*_{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 problem as in §3 but for the Lorentzian medium with the above parameters and for fixed *ξ*_{0}=13 343, which corresponds to the fixed value of found in Oughstun & Sherman (1989*b*). The boundary values at *ξ*=0 are , , where *U*(*τ*) is the unit-step function, and *F*_{0}=2. This value of *F*_{0} requires explanation, because the ‘second canonical problem’ of Oughstun and co-workers involves the same form of *F*(0, *τ*) but with *F*_{0}=1. Oughstun and co-workers do not provide an explicit expression for *B*(0, *τ*), but it is given implicitly by their requirement that there is no propagation for *x*<0 (*ξ*<0)—cf. a comment below eqn (9) in Wyns *et al*. (1989). The resulting *B*(0, *τ*) 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 *ξ*>0. That is, we *have* propagation for *ξ*<0, but for *ξ*>0 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 (1989*b*) and fig. 5 of Wyns *et al*. (1989) to see that very good agreement 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 , which corresponds to *ξ*_{0}=133 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 & Sherman—cf. fig. 10 of Oughstun & Sherman (1989*b*)—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 (1989*b*).

## 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 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 thatwhile the magnetic induction has two components:Then the time-dependent pair of Maxwell's equations take the form(5.1)which can be conveniently rewritten in matrix (or spinorial) form:(5.2)where(5.3)Equation (5.2) holds if and only if the condition of vanishing divergence of the

**field is fulfilled. This equation forms the basis of a CA difference scheme on a ‘body-centred’ grid in the**

*B**x*−

*y*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

**, provided that for every integration step**

*B**G*

_{12}remains real (this is also a sufficient condition for having

*G*

_{12}=

*G*

_{21}at every step). The matrix can be augmented with the polarization field and its associated velocity 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 Quantum Monte Carlo simulations, Suzuki found approximations to the operator in the formwith explicitly constructed *t*_{j}, such that the method is 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 Fortran 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 quantum cellular 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 algorithm to 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.

## Acknowledgments

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.

## Footnotes

The electronic supplementary material is available at http://dx.doi.org/10.1098/rspa.2006.1701 or via http://www.journals.royalsoc.ac.uk.

- Received December 2, 2005.
- Accepted February 21, 2006.

- © 2006 The Royal Society