The Method of Lines in the time domain

The Method of Lines (MoL) is a semi-analytical numerical algorithm that has been used in the past to solve Maxwell’s equations for waveguide problems. It is mainly used in the frequency domain. In this paper it is shown how the MoL can be used to solve initial value problems in the time domain. The required expressions are derived for onedimensional structures, where the materials may be dispersive. The algorithm is verified with numerical results for homogeneous structures, and for the concatenation of standard dielectric and left handed materials.


Introduction
For developing waveguide structures numerical methods are widely used, because only in exceptional cases Maxwell's equations can be solved analytically.For studying e.g. the propagation of short pulses, time domain methods are best suited.
Recently there has been a huge interest in so called left handed materials (LHM see e.g.Eleftheriades andBalmain, 2005 Caloz andItoh, 2006; other names found in the literature are e.g.double negative index material or negative index material).
In such materials the phase velocity is negative and the wave vector and the Poynting vector are anti-parallel.Due to the negative phase velocity it is not obvious that information is transmitted in forward direction.Therefore, time domain computations are performed to check if causality is violated.
As an example for such time domain methods we would like to mention the well known FDTD (finite difference time domain method see e.g.Taflove, 1995Ziolkowski and Heyman, 2001, Ziolkowski, 2003).In this method all derivatives in time and space are approximated with finite differences.
An established technique in frequency domain is the Method of Lines (MoL) where besides the discretization also analytic expressions are used.Details can be found e.g. in the book articles Pregla and Pascher, 1989;Pregla, 1995 or the book Pregla, 2008.Though the MoL is used in electromagnetics mainly as frequency domain method also time-domain computations have been made with this algorithm.In (Schiesser, 1991) the time dependency is treated with finite differences.
Another approach was followed by the group of Itoh (Nam et al., 1989a(Nam et al., ,b, 1988a(Nam et al., ,b, 1989c)).A time dependency according to sin ωt and cos ωt was used and mainly the cut-off frequencies of waveguides were determined.
Recently, the MoL in frequency domain was combined with the discrete Fourier transformation to obtain timedependent solutions (see Gerdes, 2007).
In the time domain MoL (TD-MoL) that is described in this paper, a different approach is followed.We will solve Maxwell's equation as initial value problem (as e.g. in the FDTD).
In the next section, we will derive suitable expressions.The basic principles are described in one dimension for dispersion free materials.Following, we show how material dispersion is taken into account.After that we present numerical results.To test the algorithm, we begin with a homogeneous structure.Once the correct implementation of the algorithm has been established, the concatenation of standard dielectric material with left materials is studied as further example.It should be mentioned that LHMs must be strongly dispersive (see e.g.Smith and Kroll, 2000).The paper ends with a short summary.

Dispersion free materials
In this section we will show the principles of the Method of Lines in the time domain.To develop suitable expressions, we start with Maxwell's equations: Published by Copernicus Publications on behalf of the URSI Landesausschuss in der Bundesrepublik Deutschland e.V.

S. F. Helfert: TD-MoL
For materials without dispersion (at least in the considered frequency domain) the relative permittivity resp.permeability are introduced, leading to the well known expressions: Next, the magnetic field is normalized according to with the free space impedance Z 0 and the speed of light in vacuum c 0 determined as Then, Maxwell's Equations ( 1) can be writtten as For 3-D-considerations, we would use Eqs.
In the following, however, we restrict ourselves to onedimensional problems.This has two reasons: (a) The basic principles of the algorithm are identical for 1-D-3-D problems and (b) obviously the numerical effort is lower for 1-D so that we start with this most simple case.
For one-dimensional problems with ∂/∂x = ∂/∂y = 0 we take the field components E x and H y and introduce the following abbreviations: Then, Maxwell's equations (3) become: As usual in the MoL, we discretize some of the derivatives with finite differences and use analytic expressions for the remaining ones.For the time domain algorithm, we discretize in z-direction and treat the time dependency analytically.The situation is shown in Fig. 1, where we see lines in t-direction, which may be seen as justification for the name "MoL" here.
As can be seen from Eq. ( 4) the electric and magnetic fields are coupled by first derivatives with respect to z.Therefore, as usual in finite difference methods, they are determined on different positions, i.e. shifted by half a discretization distance, as shown in Fig. 1.Since the computational window in z-direction must be finite, perfectly matched layers, as described e.g. in (Werner and Mittra, 1997) were introduced to avoid (resp.minimize) reflections at the boundaries.
For the further analysis the discretized fields are combined into a vector [F ] and the discretized derivatives into the operator matrix Q.Then the expressions in Eq. (4)become: Fig. 1.Discretization in z-direction with finite differences Note: the mathematical vectors containing the fields at different positions z were written in brackets to distinguish them from the physical vectors, which contain the spatial components of the fields and which were written in boldface.
With the exception that an analytic derivative with respect to t occurs, Equation ( 5) looks very similar to the discretized generalized transmission line equations (see e.g.Pregla, 2008 pp. 19-20) which are often used as starting point in the MoL in frequency domain.
The following steps, however, are different from the ones in frequency domain described e.g. in (Pregla, 2008).Particularly, we want to use Eq. ( 5) to solve an initial value problem with analytic expressions for E 0 (t) as input parameters.To obtain suitable expressions, we split the matrix Q into two parts as shown below: i.e. the first column (with which we multiply E 0 ) is separated from the remaining ones and the first row of Q is omitted.Then, Eq. ( 5) is rewritten in the following way: Equation ( 6) describes a coupled, linear, inhomogeneous differential equation system (IDE), where E 0 (see Fig. 1) can be considered as source.
As well known from mathematics, the general solution of linear inhomogeneous differential equations is composed of the general homogeneous solution and a particular inhomogeneous one: To obtain the homogeneous part, we must solve As usual in the MoL, this system is diagonalized by transformation to principal axes: where T and are the eigenvector/eigenvalue matrix of Q 1 .
By transforming the fields according to we obtain with the solution: Next, we must find a particular solution of the inhomogeneous differential equation system.For this purpose, we transform Q 0 and [F i ] according to: Then, we obtain from Eq. ( 6): As is known from mathematics, a particular solution can be constructed by variation of the constants, i.e. with the ansatz: The first factor comes from the solution of the homogeneous differential equation and the second one has to be determined.Introducing Eq. ( 14) into Eq.( 13) leads to hence: It should be mentioned that Q 0 just causes constant factors.Hence, analytic expressions can be given if the product between e t and E 0 (t) can be integrated in closed form.Since the first factor (e t ) only depends on the structure, the choice for E 0 (t) decides, if we are able to find such a closed form for the solution.
Then, the time dependent fields are computed with the following procedure: we begin with determining the eigenvectors/eigenvalues of the structure.It is important to note that this has to be done only once.After that, a particular solution of the inhomogeneous differential is constructed for the first time interval.Since we formulated an initial value problem, the fields at time t = 0 are known, e.g., [F in (t = 0)] = [0].With these initial values and the particular solution of the IDEs we are able to determine the homogeneous part for t = 0: Then, the field distribution at the end of this time interval (i.e. for t = τ 1 ) is determined analytically with Eqs. ( 12) and ( 14) and we can write for the whole fields: (The superscript "I" corresponds to the first time interval.)At the beginning of the second time interval, again a solution of the inhomogeneous problem is constructed and the homogeneous part is used to ensure the continuity of the fields for t = τ 1 , resulting in: This procedure (construction of a solution of the IDE, and determining the homogeneous part from the continuitycondition of the fields) is repeated for all intervals.The remaining computations to determine the field distribution at an arbitrary time are also done analytically, with very low numerical effort.For this purpose we have to store the homogenous and inhomogeneous solutions at the beginning of each time interval, i.e., [ and apply Eqs. (12 and 14).The storage itself is not very memory demanding since the homogeneous solutions are vectors, whereas the inhomogeneous ones are stored as matrices but with very few columns.For the example given in this paper (see Sect. 3) we have matrices with maximum seven columns (a sinusoidal wave was turned on/off with a polynomial of the sixth degree.) Note: Since the geometric structure does not change with time, we could have written the Eq.(17-19) also for the transformed (overlined) fields.

Dispersive material
Dispersive material is introduced into the algorithm similarly to what is done in the FDTD (see e.g.Ziolkowski and Heyman, 2001;Ziolkowski, 2003).The material dependencies are written with the polarization and the magnetization: In this paper, we use left handed material as example.The material parameters are described by the Drude model in the following way (see e.g.Ziolkowski, 2003): As before, we restrict ourselves in the following to the 1D case and introduce the following abbreviations: Then, the material parameters in Eqs.(21-22) relate E and P resp.H and M in the following way: With the help of an inverse Fourier-transform this can be rewritten to differential equations: Since Eqs.(25-26) contain derivatives with respect to t we could apply them to the TD-MoL scheme.However, only time derivatives of P and M occur and we can simplify (in terms of numerical effort) these expressions by introducing electric and magnetic currents Ziolkowski and Heyman, 2001;Ziolkowski, 2003: Then, Eqs.(25-26) together with Maxwell's equations result in the following coupled differential equation system Ziolkowski and Heyman, 2001;Ziolkowski, 2003: As can be seen, first derivatives with respect to t occur on the left hand side of these equations, whereas only derivatives with respect to z and the quantities themselves occur on the right hand side.Therefore, after discretization and combining the fields and currents in a vector we formally end up with the same system of coupled, linear, inhomogeneous differential Eq. ( 5) that we had in the dispersion free case.Obviously, [F ] and the operator Matrix Q differ here from this dispersion free case.However, for the further steps in the algorithm we may refer to Sect.2.1 and do not need to repeat them.
3 Numerical results

Input field
To validate the algorithm, we examined various structures, where a sinusoidal wave was excited, which was smoothly turned on/off.In principle, the input field E 0 is shown in Fig. 2. As can be seen, four intervals occur.Mathematically, the field in these intervals is given by the following expressions: For such input fields the integral in Eq. ( 16) can be given in closed form as is shown in the appendix.
In particular, the following parameters taken from Ziolkowski and Heyman, 2001;Ziolkowski, 2003 were used to turn on/off the sinusoidal wave: T p stands for the cycle period and is related to the frequency via ω 0 T p = 2π.The number of periods to turn on/off the signal are indicated by m on,off .In our studies we used m on = m off = 2, and a constant sinusoidal wave for m s = 4 periods so that the local time variables t 1,...,4 and the global one (t) are related in the following way:

Examined structures
As first test, the structure shown in Fig. 3 was examined.The field was injected into the homogeneous section which is terminated with perfectly matched layers (PML).Numerical results are shown in Fig. 4. The electric field is presented at different moments in times.As expected, the field travels from left to right (i.e. in positive z-direction).Since there is no dispersion, the shape of the field remains unchanged.
Furthermore, it can be seen that the fields are damped by the PMLs and that no visible reflections occur, neither at the interface between the homogeneous structure and the PML nor at the end of the device.These results suggest that the algorithm (incl.PMLs) works in principle and has been implemented correctly .
A more complicated device was studied in the following.It consists of the concatenation of standard dielectric-(i.e."right handed-(RHM)") and left handed-materials (LHM).In LHMs the material parameters ε r and µ r are negative in the interesting frequency regions.This leads to peculiar characteristics, e.g. the phase velocity is negative and waves travel "backwards".Due to physical reasons, LHMs must be strongly dispersive.In this paper, we will not discuss all these characteristic but would like to refer the interested reader e.g. to the books Eleftheriades and Balmain, 2005;Caloz andItoh, 2006 andto the Ziolkowski-papers Ziolkowski andHeyman, 2001;Ziolkowski, 2003.Since waves travel backwards in LHMs one might wonder if basic physical principles, like causality, are violated.Therefore, Ziolkowski used the FDTD and made time domain calculations Ziolkowski and Heyman, 2001; Zi- The structure under study is presented in Fig. 5.The left handed material is sandwiched between right handed material with ε r = µ r = 1.Again the structure is truncated by PMLs.For the LHM we take permittivity and permeability from the Drude-Model Eqs.(21)(22).Particularly, the following parameters from Ziolkowski and Heyman, 2001;Ziolkowski, 2003 were taken: Hence, for ω 0 (i.e. the frequency of the sine-function) the relative permittivity and permeability are: The numerical results are shown in Fig. 6.Let us begin with a few remarks.Generally, the time progresses from the top graph to the bottom one.In each of the graphs two curves are sketched, which show the situation in a short time distance.
In the top graph the electric field has just reached the first RHM-LHM interface.So, the wave has travelled to the right in the standard RHM.This behaviour is also recognized in the beginning of the LHM.In the middle graph the wave front has arrived at the second (LHM-RHM) interface.As before, the wave front travels to the right.However, it can also be observed that at the LHM input, the phase of the wave is already propagating backwards i.e. to the left.

S. F. Helfert: TD-MoL
Finally, the bottom graph shows the steady state.Here we observe a positive phase velocity in the two right handed media, but a negative phase velocity in the LHM.We can also see that the PML region damps the fields and does not cause reflections at that interface.These results agree with the ones from the FDTD and can be interpreted in two ways: 1.The FDTD findings (Ziolkowski and Heyman, 2001;Ziolkowski, 2003) are confirmed.Left handed materials do not violate the causality principle, because the wavefront propagates in positive direction.

(more important)
The TD-MoL is able to treat material dispersion.

Summary and conclusions
In this paper we presented a time domain Method of Lines.It could be shown that the algorithm works in principle, at least as long as we restrict ourselves to 1D-problems in space.Dispersive materials can be treated.A few potential advantages (compared to the FDTD) are: -The solutions in time are given in an analytic form.Therefore, we can look at field distributions at arbitrary (particularly late) times immediately.
-In the FDTD fine grid sizes (in space) also require very short time-steps to keep the algorithm stable.Since analytic expressions are used in the TD-MoL, we may look at the fields at arbitrary times without stability problems and without any loss of accuracy.
However, besides these advantages a few problems must be named as well.In principle we could apply the developed method also to higher (2-D-3-D) problems.Unfortunately, the numerical effort rises significantly when higher dimensions are considered.Therefore, the topic of the current work deals with reducing this numerical effort to be able to treat 2-D, 3-D problems in the future.

Appendix A
The turning on/off the sinusoidal waves was described by the functions given in Eq. ( 31).For computing a particular solution of the inhomogeneous differential equation system (13) we must determine [c(t)] with the integral given in Eq. ( 16).
The integral contains the product of a diagonal matrix (e t ) with a constant vector (Q 0 ).Therefore, we end up with a system of decoupled scalar equations, and we must compute integrals of the form t k e t sin ω 0 t dt e t sin ω 0 t dt For the first integral we obtain t k e t sin ω 0 t dt = t k e t e j ωt − e j ωt 2j dt

− j ω dt
As seen, the integral is solved recursively.The expression in the last row can be found in mathematical formularies.(Here we used Bronstein and Semendjajew (1979).) To compute e t sin ω 0 t dt we can introduce k = 0 into Eq.(A1).Alternatively, we find from Bronstein and Semendjajew (1979): e t sin ω 0 t dt = e t 2 + ω 2 0 ( cos ω 0 t + ω 0 sin ω 0 t) (A2) As can be shown very easily, the expressions in Eq. (A1) can be transformed into the one in Eq. (A2) in case k = 0.

Fig. 4 .Fig. 5 .
Fig. 4. Field determined at different times in the structure shown in Fig. 3

Fig. 6 .
Fig. 6.Electric field distribution in the structure shown in Fig. 5 , snapshots at three different times