Higher-order sensitivity analysis of periodic 3-D eigenvalue problems for electromagnetic field calculations

An algorithm to perform a higher-order sensitivity analysis for electromagnetic eigenvalue problems is presented. By computing the eigenvalue and eigenvector derivatives, the Brillouin Diagram for periodic structures can be calculated. The discrete model is described using the Finite Integration Technique (FIT) with periodic boundaries, and the sensitivity analysis is performed with respect to the phase shift φ between the periodic boundaries. For validation, a reference solution is calculated by solving multiple eigenvalue problems (EVP). Furthermore, the eigenvalue derivatives are compared to reference values using finite difference (FD) formulas.


Introduction
Numerical algorithms for the solution of eigenvalue problems (EVP) have been known for many years.Different numerical methods can be applied to find a discrete formulation of this problem, including the Finite Integration Technique (FIT) or the Finite Element method (FE) to only name two of them (Weiland, 1996;Schuhmann and Weiland, 2006).The focus of this paper is on parametric formulations, where the algebraic formulation depends on some structural parameter, and the eigenvalue dependency on this parameter shall be calculated in an efficient way.
For the case of a periodic repeating structure, this parameter can be the phase shift between two opposite boundaries of the unit cell of the structure.Here, the parametric evaluation of the frequency as eigenvalue leads to dispersion diagrams.An early example of such an analysis can be found in the context of linear accelerator structures in Weiland (1986).In Bandlow (2011) and Bandlow et al. (2008), the calculation of the Brillouin diagram for periodic metamaterials is real-ized using a scattering matrix approach.An approach with a Taylor approximation is shown in Klindworth and Schmidt (2014) for band structures in photonic crystals.Also, for photonic crystals, the band structure is calculated using model order reduction (MOR) in Scheiber et al. (2011).A sensitivity analysis for a waveguide eigenvalue problem has finally also been demonstrated in Burschäpers et al. (2011), here with permittivity values as parameters in the context of a simple inverse problem.
In this paper, we present a sensitivity approach for a periodic eigenvalue problem.For the sensitivity, one is interested in the derivative of the eigenvalue and the eigenvector.The first efficient method for the calculation of the eigenvector derivative, which can be adapted for large matrices, has been published by Nelson (1976) and was extended for multiple eigenvalues by Dailey (1989).We extend this method for higher-order derivatives and apply it to periodic structures.
The paper is organized as follows: in Sect. 2 the discretization method (FIT) and the formulation of the electromagnetic eigenvalue problem with periodic boundaries are shortly reviewed.Section 3 describes the theory of the eigenvalue derivatives, followed by some details of our implementation.Some comments on the required matrix derivatives are given in Sect.4, and finally, an application example is presented in Sect. 5.

Periodic boundaries in FIT
We choose the FIT to describe the discrete eigenvalue problem.The state variables are the grid voltages e and h, defined as line integrals of the electrical field and the magnetic field along the edges of the primary and the dual Grid G, G, respectively.Additionally, the grid fluxes and b occur as Published by Copernicus Publications on behalf of the URSI Landesausschuss in der Bundesrepublik Deutschland e.V. P. Jorkowski and R. Schuhmann: Higher-order sensitivity analysis surface integrals of the electric and magnetic flux densities on the primary and dual facets, respectively.
From this, we can transform Maxwell's equations, without currents and space charges, into a set of matrix-vector equations, the so-called Maxwell's Grid Equations in frequency domain: The matrices C and S represent the discrete curl and divergence operator, respectively, S is the dual divergence matrix.In a Cartesian grid with n p grid nodes, the matrix C is of size 3n p × 3n p .Using the first and second equation in Eq. (1a) leads to the discrete curl-curl eigenvalue equation, while the first equation in Eq. (1b) represents the electric divergencefree condition.For the material, we use the linear relations = M ε e and b = M µ h, where the matrices M ε and M µ are diagonal and have the size 3n p , such as the variables.Finally, ω denotes the angular frequency and j the imaginary unit.

Periodic boundaries
From these basic equations, we can define a periodic boundary condition using the Floquet condition as described e.g. in Weiland (1986).This can be done in all three coordinate directions, for the sake of simplicity we restrict here the periodic boundary to the z direction.Then, for e 1 and e 2 at the lower and the upper boundary in the z direction, respectively, we have e 1 = e 2 exp(j ϕ) with the phase shift ϕ.Thus, the tangential electric components on one side of the mesh are no independent degrees of freedom anymore and can be eliminated in the formulation, leading to a new vector e red with reduced dimension.
Formally, a sparse matrix L ϕ with e = L ϕ e red , as described in Schuhmann (2002) or Bandlow (2011), is defined to provide this transformation.It only contains 0, 1 and exp(j ϕ) entries and maps all grid voltages on themselves except for those at the upper z boundary, which are copied from the opposite side, multiplied by the phase shift exp(j ϕ).Additionally, the permittivity matrix (including averaged permittivity values ε as well as the sizes of primary edges and dual facets) has to be slightly changed in order to compensate the possible difference between the permittivity of the maped primary edges and dual faces and the original primary edges and dual faces.The new matrix is called M ε,per .Finally, the system matrices for the discrete eigenvalue formulation read: The matrix A ccp directly corresponds to the standard curlcurl formulation for electric eigenvalue problems, applied to the setup with periodic boundaries.It is well-known that this formulation contains a large number of zero eigenvalues, so-called static modes, corresponding to charges on all nonmetallic mesh nodes.To avoid these static modes, it is a standard procedure to use the grad-div-matrix A gdp as an additional operator.Due to the exact relation SC T = 0 it does not affect the searched dynamic modes, but all static solutions are shifted to non-zero eigenfrequencies within the spectrum.In order to adjust the magnitude of both expressions the matrix A gdp contains a scaling matrix D N which is described in detail in Schmitt et al. (1994) and Weiland (1985).
The shifted eigenmodes of the total matrix appear as unphysical, so-called ghost modes and have to be identified in the final results using an a-posteriori divergence check.
As a consequence, the matrix of the final EVP has full rank.In Sect.3.2 we show, why a full rank matrix is needed.Furthermore, this makes it easier to compute the fundamental modes or in general the modes with lowest (non-zero) frequencies.For the calculation of the lowest eigenvalue exist powerful algorithms.Combining the two matrices leads to the eigenvalue problem where I denotes the unity matrix.It should be noted that this formulation using FIT yields a simple EVP due to the invertibility of the permittivity matrix M ε, per (unlike a corresponding formulation using FE, where a generalized eigenvalue problem results).Furthermore, the matrix from Eq. ( 3) can be symmetrized by the similarity transformation ε, per .Therefore, only real eigenvalues are obtained.This can be justified by the fact that a loss-free structure is considered.

Eigenvalue sensitivity
We start with the right (or left) Eigenvalue problem (EVP) where A and I are the system matrix and unity matrix, respectively.Like in Nelson (1976) we choose the normalization such that the left eigenvector y H and the right eigenvector x fulfill the relation y H x = 1.

Derivative of the eigenvalue problem
The right EVP from Eq. ( 4) can be differentiated with respect to a parameter where λ and A denote the derivatives of the eigenvalue λ and the matrix A with respect to p, respectively.Multiplying Eq. ( 5) by the left Eigenvector leads to: To generalize the derivative of the eigenvalue, we differentiate n times w.r.t. the paramater p: In full analogy to the first derivative, this formula can be multiplied by the left eigenvector in order to eliminate the term with the highest order derivative of the eigenvector, which does not have to be calculated.Note, that all previous derivatives of the eigenvector are necessary.After reordering to the highest order derivative of the eigenvalue, it reads which is the algebraic equation to calculate the eigenvalue derivative.

Derivative of the eigenvectors
As mentioned above, the previous derivatives of the eigenvector are necessary.Reordering Eq. ( 5) w.r.t. the eigenvector derivative leads to: This step can be repeated recursively to obtain the nth order derivative of the eigenvector from the nth order eigenvalue and the (n − 1)th order eigenvector derivative.To obtain the linear system we need to solve Eq. ( 7) and reorder it w.r.t. the highest order derivative of the eigenvector Note that all these equations are linear systems with rank defect, e.g.A in Eq. ( 9) has rank n, but B has rank n − 1 (or less, depending on the multiplicity of λ).A strategy to cope with that rank defect is already described in Nelson (1976) for single eigenvalues and Dailey (1989) for multiple eigenvalues.

Implementation
Due to the need of the previous derivatives of the eigenvalues and eigenvectors for the following eigenvalue derivative, we obtain a recursive algorithm.Figure 1 shows the schematic procedure to calculate the dispersion diagram of the unit cell of a periodic structure.The first step is to set the phase shift ϕ 0 , the expansion point in the dispersion diagram.From this, the matrices for the eigenvalue problem can be created.
After solving this eigenvalue problem, a divergence check needs to be performed in order to find the ghost modes which are caused by the A gdp matrix.For the non-ghost modes, we can compute the derivatives of the eigenvalues recursively.We start with Eq. ( 8) to calculate the first derivative of the eigenvalue and continue with the derivative for the eigenvector if higher-order derivative are desired.The derivative for the eigenvector is calculated by solving Eq. ( 10).After the derivative of the eigenvector is known, the higher-order derivative of the eigenvalue can be calculated from Eq. ( 8).One can repeat this procedure until the desired derivative of the eigenvalue has been found.Next, the dispersion diagram can be evaluated using a Taylor expansion around the expansion point ϕ 0 .
For the nth order derivative of the eigenvalue, this requires to solve one EVP (formulation 4 with the matrix from Eq. 3), (n − 1) linear systems (Eq.10, and n algebraic Eq. 8).The third row in Table 1 shows the number of operations needed to calculate the dispersion diagram.
The main effort of this numerical algorithm arises from the number of linear equations which have to be solved in order to obtain the eigenvector derivative.From the formulas in Eqs. ( 9) and (10) one can see that all these equations share the same system matrix, which thus has to be factorized (LU decomposition) only once.Also, the new method is the more efficient compared to the other methods, the less different modes are considered.

Model order reduction
As mentioned in the previous section, the main numerical effort results from solving a high number of linear systems or performing the matrix decomposition.In order to avoid this numerical cost a MOR for the solution of the linear system in Eq. ( 10) is performed.Therefore we assume two n × p matrices V and W such that a matrix A of dimension n × n is approximated by W AV H with a p × p matrix A with p n. Accordingly, the vector x is approximated by the vector V x where x has the length p.
For numerical stability the matrices are chosen such that V = W and V H V = I.The matrices V and W are built by an Lanczos algorithm, which is for example described in Lanczos (1950) or in Saad (2000).After the projection matrices for the lefthandside of Eq. ( 9) are known, the reduced vector x n is calculated by The accuracy which can be achieved with this MOR approach depends on the dimension of its subspace.If the solution is not yet satisfactory, it can be used as starting vector for an iterative solver to increase the accuracy of the solution of Eq. ( 10).An Conjugate Gradient (CG) solver is used to increase the accuracy of the solution.A describtion for the CG solver can be found in Saad (2000) or Shewchuk (1994).

Derivative of the matrix
To calculate Eqs. ( 8) and ( 10) it is necessary to know the derivative of the matrix A p .Obviously, this matrix derivative can be calculated analytically and easily implemented.From Eq. ( 2a) we obtain for the A ccp matrix and analogously for the A gdp matrix The matrix L ϕd denotes the L ϕ -operator without the 1entries and present the derivative with respect to ϕ (without the factor j or −j for the hermitian case).
Note, that this derivative of the matrices is with respect to the parameter ϕ which only takes place in the L ϕ matrix.The algorithm presented in Sect.3.3 allows more parameter dependencies than the one which is reflected in these formulas.When changing the parameter, only the derivative of the matrix needs to be exchanged.Also, note here, that since a divergence check is performed, such that λ only depends on A ccp , the derivative of the A gdp matrix is not necessary in theory.But ignoring the A gdp matrix results in numerical errors.Alternatively, the eigenvalue needs to be recalculated after the divergence check using only A ccp , with an extra effort, to prevent this numerical error.

Results
To test the algorithm, a periodic model from the CST-library (CST, 2016) is considered.This model has 46 787 grid points, resulting in n ≈ 140 000 degrees of freedom in the eigenvalue equation.Figure 2 shows the simulation setup.The matrices are exported to Matlab (MathWorks, 2015), where the algorithm has been implemented.For the full dispersion diagram, we solve the eigenvalue problem multiple times, to obtain a reference curve.Additionally, in order to validate the derivatives, a finite differences (FD) approximation with respect to the parameter dependency is used, applying a number of separate solutions of the eigenvalue problem with slightly modified parameter ϕ.
Figure 3 shows the dispersion curves for the first three modes which have been calculated.It can clearly be seen that the results of the new algorithm fit perfectly to the ones of the FD approach.For this diagram, a fourth order Taylor approximation has been chosen.The expansion point is ϕ = π/2 and a very good agreement to the original curve is obtained near the expansion point.As expected, the deviation gets larger at points far away from the expansion point.Also, note that the numerical error gets higher, the higher the order of the approximated derivative is.This is due to the recursiveness of the algorithm.
Figure 4 shows the same dispersion curve as in Fig. 3, this time with the results of the new method with and without solving Eq. ( 10) with MOR.Although the CG solver is stopped after 1000 steps the curves are nearly identical near to the expansion point ϕ 0 .Due to the natural deviation of the Taylor approximation a more accurate solution is not necessary.
The calculation time for a different amount of modes is shown in Fig. 5 for each method.The calculation has been done on a standard PC with an non-optimized Matlab implementation.As can be seen, the new method is most efficient for a small amount of calculated modes.Using MOR with the new method solves the problem of the high calculation time for a higher amount of modes.Finally Fig. 6 shows the calculation time for different amount of the grid points.The mesh resolution have been varied in each direction.For the calculation of 4 modes the total computation time varies considerably between the different methods.At least in this example, the new method clearly outperforms the FD approach, especially when combined with the MOR technique.The drawback is an only small deviation far away from the expansion point.
P. Jorkowski and R. Schuhmann: Higher-order sensitivity analysis

Conclusions
The numerical simulation demonstrates the accuracy of the algorithm to calculate the eigenvalue derivatives.The major advantage is the reduced simulation effort compared to the reference solutions using the FD method or a large number of single eigenvalue computations.This effect gets larger if it is possible to calculate the derivatives of the matrices analytically like in this example.Augmented by the expansion with model order reduction, the new method is also very efficient for a large amount of calculated eigenvalues.
Also, the presented method can be used for different parameter dependencies, where of course further implementation effort is needed to apply it e.g. to geometric parameter variations.For problems where the parameter dependency does not allow an analytical calculation of the matrix derivative, a small additional cost in computation time has to be expected to calculate it numerically, e.g. using again FD.
Possible further research efforts will focus on sensitivities with respect to multiple parameters.Also, there are some ideas to reduce the number of required single solutions to obtain accurate dispersion diagrams, e.g., one method with a forward and backward check has been presented in Klindworth and Schmidt (2014).Another way to compute the band structure is the solution of multiple EVPs at different points with MOR, which has been presented in Scheiber et al. (2011).Some of these points can be used as expansion points ϕ i for additional Taylor approximation to minimize the deviation away from the expansion point.Finally, the method will benefit directly from all kind of improvements in the implementation, especially concerning the solution of the basic eigenvalue problem at the expansion point and the following matrix inversions.Especially if higher-order derivatives are required, it is also an important issue to find some measures to control possible numerical errors.

Figure 2 .
Figure 2. Simulation setup with periodic boundaries at the front and back plain.

Figure 3 .
Figure 3.The first three modes of the model with the fully calculated diagram (•), reference from the FD approach (∇), and results from the new algorithm ( .)

Figure 4 .Figure 5 .
Figure 4.The first three modes of the model with the fully calculated diagram (•), the results from the new algorithm ( ) and the new algorithm with MOR (∇).

Figure 6 .
Figure 6.The calculation time for the fully calculated diagram ( * ), reference from the FD approach (∇), results from the new algorithm ( ) and the new algorithm with MOR (•) in dependence of the grid points.

Table 1 .
Operations performed for each method: n is the order of the maximum order of eigenvalue derivatives, and k is the number of modes to be considered.