Figge, M.T.

Kole, J.S.

Raedt, H. De

English

We present a family of unconditionally stable algorithms, based on the Suzuki product-formula approach, that solve the time-dependent Maxwell equations in systems with spatially varying permittivity and permeability. Salient features of these algorithms are discussed. As an illustration we compute, as a function of cluster size, the spectrum of electromagnetic modes in a cluster of photonic bandgap material.</p

Dissertations of the University of Groningen

   University of GroningenSolving the time-dependent Maxwell equations by unconditionally stable algorithmsKole, J.S.; Figge, M.T.; Raedt, H. DePublished in:EPRINTS-BOOK-TITLEIMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite fromit. Please check the document version below.Document VersionPublisher's PDF, also known as Version of recordPublication date:2003Link to publication in University of Groningen/UMCG research databaseCitation for published version (APA):Kole, J. S., Figge, M. T., & Raedt, H. D. (2003). Solving the time-dependent Maxwell equations byunconditionally stable algorithms. In DP. Landau, SP. Lewis, & HB. Schuttler (Eds.), EPRINTS-BOOK-TITLE (Vol. 90, pp. 205-210). (SPRINGER PROCEEDINGS IN PHYSICS; Vol. 90). University ofGroningen, The Zernike Institute for Advanced Materials.CopyrightOther than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of theauthor(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).The publication may also be distributed here under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license.More information can be found on the University of Groningen website: https://www.rug.nl/library/open-access/self-archiving-pure/taverne-amendment.Take-down policyIf you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediatelyand investigate your claim.Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons thenumber of authors shown on this cover page is limited to 10 maximum.Download date: 03-06-2022Solving the Time-Dependent Maxwell Equationsby Unconditionally Stable AlgorithmsJ.S. Kole, M.T. Figge, and H. De RaedtCentre for Theoretical Physics and Materials Science CentreUniversity of Groningen, Nijenborgh 4,NL-9747 AG Groningen, The NetherlandsE-mail: j.s.kole@phys.rug.nl, m.t.figge@phys.rug.nl, deraedt@phys.rug.nlhttp://www.compphys.rug.nlVersion: 6 March 2002Abstract. We present a family of unconditionally stable algorithms,based on the Suzuki product-formula approach, that solve the time-dependent Maxwell equations in systems with spatially varying permit-tivity and permeability. Salient features of these algorithms are discussed.As an illustration we compute, as a function of cluster size, the spectrumof electromagnetic modes in a cluster of photonic bandgap material.1 IntroductionThe Maxwell equations describe the evolution of electromagnetic (EM) fields inspace and time [1]. They apply to a wide range of different physical situations andplay an important role in a large number of engineering applications. In manycases, numerical methods are required to solve Maxwell’s equations, either in thefrequency or time domain. For the time domain, a well-known class of algorithmsis based on a method proposed by Yee [2] and is called the finite-difference time-domain (FDTD) method. These algorithms owe their popularity mainly due totheir flexibility and speed while at the same time they are easy to implement. Alimitation of Yee-based FDTD techniques is that their stability is conditional,depending on the mesh size of the spatial discretization and the time step of thetime integration [3].In this paper we describe a family of unconditionally stable algorithms thatsolve the time-dependent Maxwell equations (TDME) [4]. The key to the con-struction of these algorithms is the observation that orthogonal approximationsto the (orthogonal) time evolution operator automatically yield unconditionallystable algorithms. The Lie-Trotter-Suzuki product formulae [5–7] provide themathematical framework to construct orthogonal approximations to the time-evolution operator of the Maxwell equations.2 TheoryWe consider EM fields in a three-dimensional medium with spatially varyingpermittivity and/or permeability, surrounded by a perfectly conducting box. Inthe absence of free charges and currents, the EM fields in such a system satisfyMaxwell’s equations [3]∂∂tB = −∇× E and ∂∂tD = ∇× H , (1)divB = 0 and divD = 0 , (2)where H = (Hx(r, t),Hy(r, t),Hz(r, t))T = µB denotes the magnetic field, µ =µ(r) is the permeability, E = (Ex(r, t), Ey(r, t), Ez(r, t))T = D/ε the electricfield, and ε = ε(r) denotes the permittivity. For simplicity of notation, we willomit the spatial dependence on r = (x, y, z)T unless this leads to ambiguities. Onthe surface of the perfectly conducting box the EM fields satisfy the boundaryconditions n × E = 0 and n · H = 0, with n denoting the vector normal to aboundary of the surface [1].Some important symmetries of the Maxwell equations (1)-(2) can be madeexplicit by introducing the fields X(t) =√µH(t) and Y(t) =√εE(t). In termsof the fields X(t) and Y(t), the TDME read∂∂t(X(t)Y(t))=(0 − 1√µ∇× 1√ε1√ε∇× 1√µ 0)(X(t)Y(t))≡ H(X(t)Y(t)). (3)Writing Ψ(t) = (X(t),Y(t))T , Eq. (3) becomes ∂∂tΨ(t) = HΨ(t) It is easy toshow that H is skew-symmetric, i.e. HT = −H, with respect to the inner prod-uct 〈Ψ |Ψ ′〉 ≡ ∫VΨT · Ψ ′ dr, where V denotes the volume of the enclosing box.The formal solution of the TDME is given by Ψ(t) = U(t)Ψ(0) = etHΨ(0)where Ψ(0) represents the initial state of the EM fields. The operator U(t) =etH determines the time evolution. By construction ‖Ψ(t)‖2 = 〈Ψ(t)|Ψ(t)〉 =∫V[εE2(t) + µH2(t)]dr, relating the length of Ψ(t) to the energy density w(t) ≡εE2(t) + µH2(t) of the EM fields [1]. As U(t)T = U(−t) = U−1(t) = e−tH itfollows that 〈U(t)Ψ(0)|U(t)Ψ(0)〉 = 〈Ψ(t)|Ψ(t)〉 = 〈Ψ(0)|Ψ(0)〉. Hence the time-evolution operator U(t) is an orthogonal transformation, rotating the vector Ψ(t)without changing its length ‖Ψ‖. In physical terms this means that the energydensity of the EM fields does not change with time, as expected on physicalgrounds [1].3 Unconditionally Stable AlgorithmsA numerical procedure that solves the TDME necessarily starts by discretizingthe spatial derivatives. This maps the continuum problem described by H onto alattice problem defined by a matrix H. Ideally, this mapping should not changethe basic symmetries of the original problem. The underlying symmetry of theTDME suggests to use matrices H that are real and skew-symmetric. Adoptingthe staggered Yee grid and central differences for the spatial derivatives [2, 3],the resulting matrix H is indeed skew-symmetric and the EM fields satisfy therequired boundary conditions [4]. As the discretization procedure is not essentialfor what follows, we omit these technicalities here and refer the reader to Ref. [4].Formally the time evolution of the EM fields on the lattice is given by Ψ(t+τ) = U(τ)Ψ(t) = eτHΨ(t). In practice, a numerical procedure solves the TDMEby making use of an approximation Ũ(t) to the true time evolution U(t). Anecessary and sufficient condition for an algorithm to be unconditionally stable isthat ‖Ũ(t)Ψ(0)‖ ≤ ‖Ψ(0)‖. In other words, the length of Ψ(t) should be bounded,for arbitrary initial condition Ψ(t = 0) and for any time t [8]. By choosing forΨ(0) the eigenvector of Ũ(t) that corresponds to the largest eigenvalue of Ũ(t),it follows that the algorithm will be unconditionally stable by construction ifand only if the largest eigenvalue of Ũ(t) (denoted by ‖Ũ(t)‖) is less or equalthan one [8]. If the approximation Ũ(t) is itself an orthogonal transformation,then ‖Ũ(t)‖ = 1 and the numerical scheme will be unconditionally stable.A systematic approach to construct orthogonal approximations to matrixexponentials is to make use of the Lie-Trotter-Suzuki formula [5, 6] and general-izations thereof [7, 9]. For instance if H =∑pi=1 HiU1(τ) = eτH1 . . . eτHp , (4)might be a good approximation to U(τ) if τ is sufficiently small and, most im-portantly, if all the Hi are real and skew-symmetric, U1(τ) is orthogonal byconstruction. Therefore, by construction, a numerical scheme based on Eq. (4)will be unconditionally stable. The Taylor series of U(τ) and U1(τ) are identicalup to first order in τ . The product-formula approach provides simple, system-atic procedures to improve the accuracy of the approximation to U(τ) withoutchanging its fundamental symmetries. For example the orthogonal matrixU2(τ) = U1(−τ/2)T U1(τ/2) = eτHp/2 . . . eτH1/2eτH1/2 . . . eτHp/2 , (5)is correct up to second-order in τ [7, 9]. Suzuki’s fractal decomposition ap-proach [7] gives a general method to construct higher-order approximationsbased on U1(τ) or U2(τ). A particularly useful approximation, correct up tofourth-order in τ , is given by [7]U4(τ) = U2(aτ)U2(aτ)U2((1− 4a)τ)U2(aτ)U2(aτ) , (6)where a = 1/(4−41/3). The approximations Eqs. (4), (5), and (6) have proven tobe very useful in many applications [6, 9–18] and turn out to be equally usefulfor solving the TDME. In practice, an efficient implementation of the first-orderscheme is all that is needed to construct the higher-order algorithms Eqs. (5)and (6).4 Simulation ResultsNumerically the existence of photonic bandgaps is most easily studied by adopt-ing periodic boundary conditions. However, in nature periodic boundary condi-tions are very hard to realize. Therefore it is of interest to study the dependenceFig. 1. Left: The 3D system: group of connected spheres of air (grey) in a dielectricbackground (white). The centers of the spheres form an FCC structure Right: Localdensity of states as function of frequency for 3D systems consisting of n spheres. PBCdenotes the reference system with periodic boundary conditions.of photonic bandgaps as a function of system size with physically realizableboundary conditions. The method we outline above is well-suited for this pur-pose. It does not rely on periodic boundary conditions nor on the use of theFourier transform of the unit cell. We consider three-dimensional clusters oftouching spheres filled with air in a dielectric background. The centers of thespheres form an FCC lattice, see Fig.1. The procedure we use to compute the(Local) Density of States (L)DOS is described in detail in Ref. [19].First we compute the DOS for the system with periodic boundary conditionsfor a reference system containing 8 x 8 x 8 FCC unit cells of linear size λ.The FCC lattice constant is set at a = 3λ and, combined with a mesh sizeδ = 0.1λ, ensures that the spheres of air are discretized with sufficient detail. Thepermittivity of the dielectricum ε = 11.9ε0, where ε0 is the dielectric constantof vacuum. The DOS of this system clearly exhibits one photonic bandgap, seeFig.1 (dashed line) at ωa/2πc = 0.8 where c is the speed of light in vacuum.Then we compute the LDOS of the FCC cluster of n spheres filled with air(without periodic boundary conditions) by using as an initial state, a randomwave localized within the central sphere. In Fig.1 we also show the results forn = 19, 249, 1409. It is clear that the finite-sized clusters also exhibit a photonicbandgap and that the strength of the bandgap strongly depends on n (note thelogarithmic scale of the y-axis). Roughly speaking the minimum of the LDOSapproaches its limiting value as e−α√n where α is a material specific constant.5 ConclusionWe have described a new family of algorithms to solve the time-dependentMaxwell equations. Salient features of these algorithms are [4, 20]:– rigorously provable unconditional stability for 1D, 2D and 3D systems withspatially varying permittivity and permeability,– the order of accuracy in time and space can be systematically increasedwithout affecting the unconditional stability,– the exact conservation of the energy density of the electromagnetic field,– easy to implement in practice.In view of the generic character of the approach discussed in this paper, it canbe used to construct unconditionally stable algorithms that solve the equationsfor e.g. sound, seismic and elastic waves as well.AcknowledgementsThis work is partially supported by the Dutch ‘Stichting Nationale ComputerFaciliteiten’ (NCF).References1. M. Born and E. Wolf, Principles of Optics, (Pergamon, Oxford, 1964).2. K.S. Yee, IEEE Transactions on Antennas and Propagation 14, 302 (1966).3. A. Taflove and S.C. Hagness, Computational Electrodynamics - The Finite-Difference Time-Domain Method, (Artech House, Boston, 2000).4. J.S. Kole, M.T. Figge and H. De Raedt, Phys. Rev. E 64, 066705 (2001).5. H.F. Trotter, Proc. Am. Math. Soc. 10, 545 (1959).6. M. Suzuki, S. Miyashita, and A. Kuroda, Prog. Theor. Phys. 58, 1377 (1977).7. M. Suzuki, J. Math. Phys. 26, 601 (1985); ibid 32 400 (1991).8. G.D. Smith, Numerical solution of partial differential equations, (Clarendon, Ox-ford, 1985).9. H. De Raedt and B. De Raedt, Phys. Rev. A 28, 3575 (1983).10. H. De Raedt, Comp. Phys. Rep. 7, 1 (1987).11. H. Kobayashi, N. Hatano, and M. Suzuki, Physica A 211, 234 (1994).12. H. De Raedt, K. Michielsen, Comp. in Phys. 8, 600 (1994).13. A. Rouhi, J. Wright, Comp. in Phys. 9, 554 (1995).14. B.A. Shadwick and W.F. Buell, Phys. Rev. Lett. 79, 5189 (1997).15. M. Krech, A. Bunker, and D.P. Landau, Comp. Phys. Comm. 111, 1 (1998).16. P. Tran, Phys. Rev. E 58, 8049 (1998).17. K. Michielsen, H. De Raedt, J. Przeslawski, and N. Garcia, Phys. Rep. 304, 89(1998).18. H. De Raedt, A.H. Hams, K. Michielsen, and K. De Raedt, Comp. Phys. Comm.132, 1 (2000).19. A. Hams and H. De Raedt, Phys. Rev. E 62, 4365 (2000).20. J.S. Kole, M.T. Figge and H. De Raedt, submitted to Phys. Rev. E

Solving the time-dependent Maxwell equations by unconditionally stable algorithms

https://pure.rug.nl/ws/files/10065781/2003CompSimStudKole.pdf

Solving the time-dependent Maxwell equations by unconditionally stable algorithms

Abstract

Similar works

Full text

Available Versions

Dissertations of the University of Groningen