Efficient determination of the left-eigenvectors for the Method of Lines

The efficient determination of left eigenvectors in the method of lines (MoL) is described in this paper. The electromagnetic fields are expanded into eigenmodes and the eigenmodes are determined from an explicit matrix eigenvector problem. To study complicated structures with a moderate numerical effort, the analysis is done with a reduced set of these eigenmodes. The enforcements of the continuity of the transverse electric and magnetic fields at interfaces leads to expressions with rectangular matrices. Now left eigenvectors can be considered as inverse of these rectangular matrices. Until now, the left eigenvectors were determined from a second explicit eigenvalue problem. Here, it is shown how they can be determined with simple matrix products from previously determined right eigenvectors. This is done by utilizing the relation between the transverse electric and magnetic fields. The derived formulas hold for structures with Dirichlet, Neumann or periodic boundary conditions and the materials may be lossy. Open structures are modeled with perfectly matched layers (PML). To verify the expressions, various devices that contain such PMLs and lossy metals were studied. In all cases, error measures show that the algorithm derived in this paper works very well.


Introduction
Complicated waveguide structures can be analyzed analytically only in exceptional cases.Therefore, usually numerical methods are used for this purpose.There are various possibilities to classify these methods, e.g. according to their analytical/numerical part.In the well known finite difference time domain method (FDTD, see e.g.Taflove, 1995) all derivatives that occur in Maxwell's equations are discretized with finite differences.Since all points in space and time are dis-cretized, the simulation of complicated structures can be numerically demanding.
On the other side there are eigenmode algorithms (see e.g.Sudbø, 1993;Sztefka, 1992).Here, analytic expressions are used at least in the direction of wave propagation.Continuously varying structures are typically modeled with a stair-case approximation.Besides, usually the analytic expressions exist of infinite series, which must be truncated in practical applications.Therefore, in spite of its analytic formulation, approximations are also introduced in case of eigenmode methods.There exist various ways in computing the eigenmodes and it is not always easy to determine all important ones with analytic expressions in case of complex structures.
Here, we use the Method of Lines (MoL) where the eigenmodes are determined after discretization with finite differences from an explicit eigenvalue problem.
Since the MoL is well documented in the literature we just mention some book chapters (Pregla and Pascher, 1989) (Pregla, 1995) and the book (Pregla, 2008), here.For threedimensional structures a two-dimensional discretization is required and the occurring matrices become very large.To keep the numerical effort moderate, it was proposed in the past to use a reduced number of eigenmodes (Gerdes, 1994).At interfaces between different sections the transverse field components have to be continuous.When a reduced set of eigenmodes is used one has to determine inverses of rectangular (not square) matrices.For this reason, (Schneider, 1999) proposed the use of a pseudo inverse (Strang, 1986).In contrast, the utilization of left-eigenvectors was suggested in (Helfert et al., 2003).The given methods work quite well, but the numerical cost for determining the pseudo inverse or the left eigenvectors is quite high.Here, we will also apply lefteigenvectors but show how they can be determined ple matrix products from previously determined right eigenvectors.The derivations can also be understood as proof of the orthogonality of eigenmodes in closed structures.
The paper is organized as follows: we start with a repetition of the principles of the MoL.Following is a description of the left-eigenvectors, and we will particularly show how they can be determined in an efficient way.After that, numerical results are shown where we demonstrate that the developed expressions also work in case of perfectly matched layers and for lossy materials.The paper ends with a summary.

Method of Lines
In the Method of Lines (MoL) the wave propagation is described with eigenmodes.It allows the analysis of complicated devices.As example the polarization converter taken from (Mustieles et al., 1993) is shown in Fig. 1.It consist of a dielectric waveguide at the input followed by a region with a periodic modification of the core.With a suitable choice of the length , a polarization rotation is possible (see Mustieles et al., 1993).
The eigenmodes in the MoL are determined after a discretization with finite differences.In this section, we describe the basic algorithm.To model waveguide structures with the MoL, we first divide it into homogeneous sections with respect to the direction of propagation as shown in Fig. 2.Then, solutions for the homogeneous sections and the interfaces are determined.

Solution in homogeneous sections
For the numerical analysis, we start with generalized transmission line (GTL) equations, which relate the transverse electric and magnetic fields.In what follows, the coordinates are normalized with the free space wave number k 0 = ω √ ε 0 µ 0 = 2π/λ 0 : u = k 0 u (u = x, y, z).Further, the magnetic field is normalized with the free space wave impedance Z 0 = √ µ 0 /ε 0 = 120π : H u = Z 0 H u .Then, the GTL-equations read (see e.g.Pregla, 1999;Pregla, 2008): with Note: the continuous physical vectors were put in brackets here, to distinguish them from mathematical vectors that originate from discretization.These mathematical vectors are written in bold.
In the analysis, anisotropic materials of the following form are considered: Adv. Radio Sci., 13, 19-29, 2015 www.adv-radio-sci.net/13/19/2015/With these tensors for the material parameters the operators R H und R E are given as: x Note: in the following the derivatives in x, y direction will be treated differently than the derivative with respect to z.To indicate this, the abbreviation D v = ∂/∂v (with v = x, y) was introduced for these derivatives.
By combining the GTL-equations (1), the following coupled wave equations are obtained: with The next step of the analysis is the discretization of Eq. ( 4) with finite differences in the direction of the cross-section, while the z-dependency remains.Figures 2-3 show this discretization process from the top and in the cross-section.
As indicated, the components of the electric or magnetic fields (i.e.E x , E y resp.H x , H y ) are determined at different positions.With this shift, the coupling between the components that is described by the product of first derivatives in x and y direction (see off-diagonal elements of R E,H ) can be modeled in an optimal way.[More details can be found e.g. in Pregla, 2008].Now, the fields are combined in vectors, and the operators Q, R becomes operator matrices that contain the approximations of the derivatives with finite differences and the material parameters.
In what follows, we describe the solution of the wave equation for the electric and magnetic field in parallel.However, these fields are coupled by Eq. (1).Therefore, in practice, we have to solve the wave equation for only one of the fields and determine the remaining one with the help of Eq. (1).After discretization, Eq. (4) becomes: The expressions in Eq. ( 6) are coupled linear differential equation systems.As known from mathematics, such equation systems can be solved by transformation to principal axes with T E,H and 2 being the eigenvectors and eigenvalues of Q E,H .
Here, the matrices Q E and Q H were determined from Eq. ( 5) and then discretized with finite differences.However, identical matrices are obtained if the operators R E and R H are discretized first and the multiplication is done for these discrete matrices.Now, as known from mathematics (see e.g.Zurmühl and Falk, 1984, pp. 163) the eigenvalues of matrices that are determined from the product of two square matrices in reversed order are identical.For this reason Q E and Q H have the same eigenvalues and the subscript E or H was omitted for 2 .
After transforming the fields according to: a system of decoupled equations is developed: whose solution can be given immediately as: Equations ( 7)- (11) show that the solution of the wave equation is given in terms of eigenmodes.Each column of the eigenvector matrices T E,H represents the electric resp.magnetic fields of these eigenmodes, is a diagonal matrix with the propagation constants and the overlined quantities E f,b H f,b represent the amplitudes of the forward and backward propagating modes.By introducing the solution for the electric or magnetic field Eqs. ( 10), (11) into the GTLexpressions in discretized domain Eq.(1) one can obtain following relation between the amplitudes: In this case, the transformation matrices are related according to:

Interfaces
After solving the wave equation in homogeneous sections, interfaces must be considered.Here, the transverse electric and magnetic field components have to be continuous.transverse fields, a relation for the amplitudes of the eigenmodes is obtained, where Eqs. ( 8)-( 12) were considered From Eq. ( 14) we can e.g.compute the amplitudes of the eigenmodes in region II from the ones in region I in the following way: with In principle, we are now in position to analyze whole devices like the polarization converter shown in Fig. 1.For this purpose, we must introduce conditions at the input and output of the device and apply the expressions derived for the homogeneous sections and the interfaces.However, the exponential increasing terms in Eqs. ( 10)-( 11) can cause numerical problems if e.g.transfer matrices are used.Therefore, numerical stable expressions have been developed in the past, to avoid the numerical computation of these exponential increasing terms.These are algorithms with impedances/admittances (see e.g.Rogge and Pregla, 1993;Pregla, 2008) or with scattering parameters (e.g.Helfert and Pregla, 2002).At this point, we do not go into details, but just refer to the relevant literature.

Reduction of the eigenmode system
The presented algorithm works well, if the number of discretization lines is not too high.However, particularly if vectorial, 3-D problems are considered, the occurring matrices for the electric/magnetic field distribution become very large.
The determination of all eigenmodes in Eq. ( 7) is very time consuming and the memory requirement becomes very high.
For this reason, the analysis with a reduced number of eigenmodes was proposed in Gerdes (1994).However, in that paper all eigenmodes were determined first.Only after that, some of them were omitted for the further simulations.
In Schneider (1999) was described how this computation of all eigenmodes can be avoided.For this purpose the Arnoldi algorithm that is implemented in the MATLAB program package (MATLAB, 2014) was used.This was done here as well.Information about this method can be found in the literature (e.g.Arnoldi, 1951;Golub and van Loan, 1989) so that we do not go into details here.
With a reduced set of eigenmodes the continuity of the transverse fields at interfaces leads to expressions where the matrices have the following shape: cretization lines is not too high.However, particularly if vectorial, 3D problems are considered, the occurring matri-155 ces for the electric/magnetic field distribution become very large.The determination of all eigenmodes in Eq. ( 7) is very time consuming and the memory requirement becomes very high.For this reason, the analysis with a reduced number of eigenmodes was proposed in (Gerdes, 1994).However, in 160 that paper all eigenmodes were determined first.Only after that, some of them were omitted for the further simulations.
In (Schneider, 1999) was described how this computation of all eigenmodes can be avoided.For this purpose the Arnoldi algorithm that is implemented in the MATLAB pro-165 gram package (MATLAB, 2014) was used.This was done here as well.Information about this method can be found in the literature e.g.(Arnoldi, 1951), (Golub and van Loan, 1989) so that we do not go into details here.
With a reduced set of eigenmodes the continuity of the 170 transverse fields at interfaces leads to expressions where the matrices have the following shape:  15)-( 17) require the inversion of non-square matrices.The use of a pseudo-inverse (see Strang, 1986) was proposed in Schneider (1999) for this purpose.Here, we follow a different path and apply left eigenvectors as described in Helfert et al. (2003).
In principle the left eigenvectors Y of a matrix A are solutions of in contrast to the more known problem for the right eigenvectors The eigenvalues λ for both problems are identical.After suitable sorting and normalization of the left and right eigenvectors we may write: where I is the identity matrix.For our purpose, it is important that this relation is also true when a reduced set of eigenmodes is used.This is indicated in Fig. 5. Note: Eq. ( 20) is true as long as the eigenvalues are different.For degenerated modes with identical eigenvalues this condition does not have to be fulfilled.However as shown e.g. in Zurmühl and Falk (1984), pp.155 it is possible to transform the eigenvectors in such a way that condition Eq. ( 20) is enforced.A summary of this transformation is given in the Appendix.
Once the left eigenvectors (indicated by the superscript "L") have been computed, we determine the expressions given in Eq. ( 17) as: The superscript "r" indicates that a reduced set of eigenmodes was used.It is now worth noting that the matrices in Eq. ( 21) are sub-matrices of the ones given in Eq. ( 17).The illustration shown in Fig. 5 can be interpreted as a special case of this feature.
Unfortunately, the author is not aware of a standard algorithm that determines right-and left-eigenvectors at the same time.As consequence, in the past, the right and left eigenvectors were determined independently of each other (Helfert et al., 2003).This was done with Q E,H and its transposed.Obviously, the numerical effort is quite high.As second problem the numerical orthogonality of the left and right eigenvectors Eq. ( 20) can be quite bad when they are computed independently.

Efficient determination of the left eigenvectors
In this section we describe how we can avoid solving the eigenvalue problem twice.To develop suitable expressions, we start with the operators R E,H .For Dirichlet, Neumann, or periodic boundary conditions, we may write the discretized operators with matrix products of the following form: The matrices D x , D y and their transposed (indicated by the superscript "T") contain the finite difference approximation of the derivatives with respect to x and y. ε x,y,z and µ x,y,z are diagonal matrices containing the permittivities and permeabilities on the discretization points.As can be seen, the matrices are symmetric (not Hermitian) Next we look at the matrices that occur in the wave equation Eq. ( 6).As mentioned before, we may write If we now transpose e.g.Q H , we find, due to Eq. ( 24), the following relation: Next, this transposed matrix is diagonalized.We use the left eigenvector matrix (indicated by the superscript "L", as before) as inverse of T H and obtain Similar, the diagonalization of Q E reads: A comparison of Eq. ( 26) and Eq. ( 27) with consideration of Eq. ( 25) shows that T E , T H and the corresponding left eigenvectors (i.e.their inverse matrices) are related according to: Note: since eigenvectors can be determined up to an arbitrary scaling factor, Eq. ( 28) is only true, if the eigenvectors have been normalized accordingly.This is not a principle problem, but has to be kept in mind.
As final step we introduce Eq. ( 28) into Eq.( 13) and obtain: The products in Eq. ( 29) can be computed with verly low numerical effort, because R H and R E are sparse-matrices.Further, the eigenvalues are combined in the diagonal matrix so that its inverse can be computed very easily.Hence, Eq. ( 29) shows that the left eigenvectors can be determined with simple matrix products from the previously determined right eigenvectors.The numerical effort for this procedure is clearly lower than that for solving the eigenvalue/eigenvector problem twice or as determining a pseudo inverse.

Open structures
As mentioned before, the derived expression Eq. ( 29) only holds for special boundary conditions, for which the operator matrices R E,H are symmetric.Now outgoing waves are completely reflected at Dirichlet or Neumann boundaries ("hard boundaries") leading to modeling errors.There exist various ways to reduce such reflections, e.g.absorbing boundary conditions (ABC) as described in Pregla (1995) or Vassallo and Collino (1997).Unfortunately, with the ABCs mentioned above, the symmetry condition in Eq. ( 24) does not hold any longer.Hence, we need to apply methods that permit the inclusion of hard BCs.
Particularly, we decided to use perfectly matched layers (PMLs).These are lossy layers positioned at the boundary of the computational window.Waves that enter this region are damped and after being reflected at the hard wall, they are www.adv-radio-sci.net/13/19/2015/Adv.Radio Sci., 13, 19-29, 2015 damped further to negligible amplitudes before re-entering the original computational window.Generally, PMLs can be interpreted in various ways (see e.g.Berenger, 1994;Mittra and Pekel, 1995;Rappaport, 1995;Al-Bader and Jamid, 1998).Here, they are considered as lossy anisotropic layers, which are described by tensors that only contain non-zero elements on the main diagonal (see e.g.Sacks et al., 1995;Werner and Mittra, 1997;Cucinotta et al., 1999), i.e. tensors of the form given by Eq. ( 3), for which the formulas were derived.Hence, Eq. ( 29) is also true for PMLs.

Orthogonality of the eigenmodes
By considering the conservation of energy and from the reciprocity theorem one can derive orthogonality relations for the eigenmodes in waveguide sections (see e.g.Syms and Cozens, 1992;Snyder and Love, 1983).For the fields of two different eigenmodes (labels i, k) that propagate in z direction the orthogonality condition reads Note: as before, the physical vectors ([E], [H ],[e]) were written in brackets.
In the MoL, the discretized field distributions of the eigenmodes are combined in the columns of the matrices T E , T H . Therefore, the discretized integral in Eq. ( 30) contains the product which results in a diagonal matrix, as derived.Hence, the orthogonality of the discretized eigenmodes in a closed structure can be seen directly from Eq. ( 28).We should point out that this orthogonality relation does also hold for lossy structures.

Numerical results
In this section we show some numerical examples to evaluate the derived expression.For comparison the left eigenvectors were (a) determined with Eq. ( 29) and (b) by solving a second eigenvalue/eigenvector problem as described in (Helfert et al., 2003).For a quantitative assessment, an error measure was introduced.Ideally, the product between left-and righteigenvectors results in an identity matrix.The deviation from this ideal result can be written as: In the following we will use ||M err || max i.e. the maximum element of this matrix as error measure.
There are a few factors that contribute to this error.The first one is caused by the finite machine precision so that the eigenvectors can be determined only up to a certain accuracy.A second contribution to the error comes from degenerated eigenmodes, as mentioned in Sect.2.2.As will be shown numerically in this section this error can be reduced by a transformation of these modes.
To see if the developed expressions work in principle and to study the behavior of the PMLs, we started with two simple 2-D structures (TE-polarization with the components E y , H x , H z ).For both cases Dirichlet conditions were introduced at the lateral boundaries and results obtained with and without PMLs were compared.As mentioned before, the PMLs were introduced as anisotropic materials (18 lines close to the boundaries) where the following values for the materials were chosen: µ z = ε y = 1 − 0.36j , µ x = 1/µ z .For details about the relation between these values see e.g.Sacks et al. (1995), Werner andMittra (1997), andCucinotta et al. (1999).
We should point out that the z direction is treated with analytic expressions.Hence, to model the infinite long structures no PMLs are needed, but we assume that only forward propagating modes occur for these regions in Eqs. ( 10)-( 11), thus: E b = H b = 0. Since we are dealing with 2-D problems, a reduction of the eigenmode system is not necessary, and all eigenmodes were used for the computations.Particularly, the number of discretization points (= number of eigenmodes) was 107.
First of all, a tilted Gaussian beam was injected into a homogeneous air region.The input field is given as: with λ 0 = 1.55 µm, w = 2 µm and ϕ = 45 • .To compute the wave propagation of the Gaussian beam, we discretize the input field in x direction and must invert Eq. ( 8) after that.As mentioned before, an infinite long region with only forward propagating modes is assumed.Hence, the inversion (written with left eigenvectors) reads: In the following, the propagation in the region is computed with Eq. ( 10) and the original fields are obtained from Eq. ( 8).As can be seen, even in this simple case left eigenvectors are used.
As second example the abrupt ending of a waveguide was studied.As can be seen in Fig. 6, we have the concatenation of a waveguide section and an infinite long air region.In the simulations the fundamental waveguide mode was injected on the left.Then, the wave propagation was simulated with the expressions given in Sect. 2. For this structure the left eigenvectors were needed to enforce the continuity of the fields at the interface between the waveguide and the air region.
Before looking at the fields, let us examine the error.The structure in Fig. 6 consists of two parts (waveguide, homogeneous air region).Therefore, for each of these regions a different error value is obtained.For the further considerations Table 1.Error determined with Eq. ( 31) for the structures shown in Fig. 6; in case "PML symmetric a)" the eigenmodes were taken as obtained by the eigenvalue procedure; in "PML symmetric b)" the degenerated eigenmodes were transformed.the worse one was taken.Since the homogeneous air region is already included in these studies, there are no additional error measures required for the Gaussian beam problem.
In the simulations, several parameters were varied.We should repeat that the expression in Eq. ( 29) was derived for a general case incl.losses.These losses can be caused by dielectric media (i.e. with complex permittivity) or by PMLs (anisotropic media with complex permittivity and permeability).Hence, the influence of the losses on the error is of particular interest.
The results for the error are summarized in Table 1.As can be seen, two cases were compared; for a lossless core (n core = 1.4) and a lossy one (n core = 1.4-0.01j).The lossless case (row "electric wall") leads to the small error
The next rows show the influence of the PMLs on the error.In row "PML symmetric a)" the determined error-value was computed directly after determining the right and left eigenvectors with Eqs. ( 7) and ( 29).Here, degenerated eigenmodes (i.e.pairs with identical propagation constant) occur, resulting in very high error values.This error can be reduced drastically by a transformation of the eigenmodes leading to the results "symmetric b)".
Since the error in case of PMLs is caused by degenerated modes, it was examined next, what happens if a small lossless case (row " electric wall") leads to the small error Å ÖÖ Ñ Ü ½ ¼ ½¾ , which deteriorates only moderately for a lossy core.
Error determined with Eq. 31 for the structure shown in Fig. 6; in case 'PML symmetric a)' the eigenmodes were taken as obtained by the eigenvalue procedure; in 'PML symmetric b)' the degenerated eigenmodes were transformed The next rows show the influence of the PMLs on the error.

315
In row "PML symmetric a)" the determined error-value was computed directly after determining the right and left eigenvectors with Eqs. ( 7) and ( 29).Here, degenerated eigen- asymmetry in the PMLs is introduced.In particular, on the lower boundary, the imaginary part of µ z was increased by the factor 1.0001 and the other values (ε y , µ x ) were modified accordingly.Then, the degeneration of the modes vanishes, leading to the error presented in row: "PML asymmetric".These errors are close to the lossless case.Keeping in mind that the values in row "PML symmetric a)" can be improved to "PML symmetric b)", it is found that all errors are quite small.Hence the formulas Eq. ( 31) work without problems for the 2-D-case incl.losses due to lossy materials or PMLs.
For comparison, the left eigenvectors were also determined from solving a second eigenvalue problem.Here, the error for the lossless case is reduced to 10 −15 .For all other cases, however, the error is similar or slightly higher than given in Table 1.Since the error in case of PMLs is caused by degenerated modes, it was examined next, what happens if a small asymmetry in the PMLs is introduced.In particular, on the lower boundary, the imaginary part of Þ was increased by the factor 1.0001 and the other values ( Ý , Ü ) were modified accordingly.Then, the degeneration of the modes vanishes, leading to the error presented in row: "PML asymmetric".These errors are close to the lossless case.
Keeping in mind that the values in row "PML symmetric a)" can be improved to "PML symmetric b)", it is found that all errors are quite small.Hence the formulas Eq. ( 31) work without problems for the 2D-case incl.losses due to lossy materials or PMLs.
For comparison, the left eigenvectors were also determined from solving a second eigenvalue problem.Here, the error for the lossless case is reduced to ½¼ ½ .For all other cases, however, the error is similar or slightly higher than given in Table-1  The computed electric field distributions are presented in Fig. 7 for the Gaussian beam and in Fig. 8 for the abrupt ending of the waveguide.The PML-graphs were obtained with symmetric PMLs.However, there are no visual differences recognizable when a small asymmetry is introduced.
For both structures it can be seen, how the reflections caused by the Dirichlet boundaries Fig. 7a (resp.Fig. 8a) can be suppressed with the PMLs Fig. 7b (resp.Fig. 8b).
The results presented so far show that the algorithm works in principle in 2-D including losses.However, the analysis of such a 2-D-structures is numerical not very demanding and (as was done) the full eigenmode system can be used.
Therefore, the real reduction of the numerical effort is expected for 3-D-structures.Two examples were considered here.The polarization converter shown in Fig. 2 had already been examined with the MoL (Helfert et al., 2003).In that study the left eigenvectors were determined from an independent eigenvalue/eigenvector computation.The structure was analyzed with N x ×N y ×2 = 100×75×2 = 15 000 discretization points (factor "2" because of the vectorial analysis) and 300 eigenmodes.When the left eigenvectors were determined with the algorithm described here, the value 10 −8 was obtained for the error Eq. ( 31).The separate left eigenvector determination (Helfert et al., 2003) results in the error value 10 −7 .Both values are satisfactory.We should point out here, that most of the CPU-time for analyzing the polarization converter is required for the determination of the eigenmodes.In particular, computing the left eigenvectors separately, required 3 min, which could be reduced to 4 s with the algorithm presented here (computations on a PC).This shows the significant reduction of the numerical effort.
The features of this polarization converter had been studied with the MoL in detail in (Helfert et al., 2003).Since these results do not change, when the eigenvectors are determined differently, they are not repeated here.Instead we refer to Helfert et al. (2003).
As second example for a 3-D-device, an array of hollow waveguides as shown in Fig. 9 was examined.The idea is to use such hollow waveguide arrays as polarization converting elements.For this purpose, the difference of the propagation constants of the vertically and horizontally polarized eigenmodes is utilized ("half-wave plate") and the length L z of the device was chosen in such a way that these wave experience a phase difference of π .Hence: A detailed description of the structure is given in Helfert et al. (2014) and Helfert et al. (2015).Therefore, we concentrate here on the results as far as the numerical algorithm is concerned and do not repeat all the features of this device.Since the permittivity of the metal is very different in optics and for THz-frequencies, studies for these two frequency regimes were performed.Particularly, silver was taken as metal.Its relative permittivity was determined with a Dudemodel, resulting in: -THz (λ 0 = 0.48 mm, f = 0.625 THz) As can be seen, losses in the metal are considered.
With the values for the permittivity given above, the eigenmodes were computed.Then, the following lengths were determined from Eq. (34): Now, plane waves were injected into the structures as shown in Fig. 10.The magnetic field at the output is sketched in Fig. 11.The rotation of the fields is clearly seen.In the THz case (Fig. 11a) we recognize that the whole field is inside the air-region, whereas a non-negligible part of the field is inside the metal in optics (Fig. 11b).

Summary and conclusions
In this paper was shown how left eigenvectors can be determined very efficiently.For this purpose the relation between the transverse electric and magnetic fields of the eigenmodes was utilized.It was found that the matrix of the left eigenvectors can be computed by simple matrix products from previously determined right eigenvalues in case of hard walls The computations were done with AE Ü ¢AE Ý ¢¾ ½ ¼¢ ½¿¼¢¾ ¾ ¼ ¼ discretization points.The number of eigen-395 modes was 20 i.e. much lower.With these parameters, the following error values occur: Å ÖÖ Ñ Ü ½ ¼ (THz), Å ÖÖ Ñ Ü ½ ¼ (optics).When the left eigenvectors were determined separately (see (Helfert et al., 2003)), errors of the same order were obtained.So, Eq. ( 29) does also work 400 very well for the considered full vectorial 3D-structures.

Summary and conclusions
In this paper was shown how left eigenvectors can be determined very efficiently.For this purpose the relation between the transverse electric and magnetic fields of the eigenmodes 405 was utilized.It was found that the matrix of the left eigenvectors can be computed by simple matrix products from previously determined right eigenvalues in case of hard walls (Dirichlet, Neumann) or periodic boundaries.Compared to the erlier determination from a second eigenvalue problem, 410 the numerical effort could be reduduced significantly.
Absorbing boundary condition that are often employed in the MoL, cannot be used here, because of symmetry require- (Dirichlet, Neumann) or periodic boundaries.Compared to the erlier determination from a second eigenvalue problem, the numerical effort could be reduduced significantly.
Absorbing boundary condition that are often employed in the MoL, cannot be used here, because of symmetry requirements.Therefore, open structures were modeled with perfectly matched layers where the PMLs were introduced as lossy anisotropic materials.The numerical results (incl.error measures) for various structures in 2-D and for 3-D show that the method works very well also when losses are considered.

Figure 1 .
Figure 1.Polarization converter as example of a complex threedimensional structure.

Figure 2 .Figure 3 .
Figure 2. Analysis with the MoL: dividing the structure into homogeneous sections and discretization.

Figure 5 .
Figure 5. Multiplication of left-and right-eigenvectors resulting in an identity matrix.

Figure 6 .
Figure 6.Abrupt ending of a waveguide.

Fig. 7 .
Fig. 7. Injection of a tilted Gaussian beam into a homogeneous section, computed electric field, a) without PMLs, b) with PMLs, the horizontal lines indicate the boundaries of the PMLs

Figure 7 .
Figure 7. Injection of a tilted Gaussian beam into a homogeneous section, computed electric field, (a) without PMLs, (b) with PMLs, the horizontal lines indicate the boundaries of the PMLs.

Fig. 8 .
Fig. 8. Electric field in the structure shown in Fig. 6; a) without PMLs, b) PMLs included, the horizontal lines indicate the boundaries of the PMLs .

Figure 8 .
Figure 8. Electric field in the structure shown in Fig. 6; (a) without PMLs, (b) PMLs included, the horizontal lines indicate the boundaries of the PMLs.

Figure 10 .
Figure 10.Injection of a plane wave into an array of hollow waveguides, magnetic field in front of the waveguides.

S
Magnetic field at the output of the hollow waveguides, a) THz-frequencies, b) optics

Figure 11 .
Figure 11.Magnetic field at the output of the hollow waveguides, (a) THz-frequencies, (b) optics.