An auxiliary field approach for computing optical resonances in dispersive media

We report on an auxiliary field approach for solving nonlinear eigenvalue problems occurring in nano-optical systems with material dispersion. The material dispersion can be described by a rational function for the frequency-dependent permittivity, e.g., by a Drude-Lorentz model or a rational function fit to measured material data. The approach is applied to compute plasmonic resonances of a metallic grating.


Introduction
Detailed knowledge on the resonant states of nano-optical systems is essential for understanding the physical properties of the systems and for designing related photonic devices [1][2][3].With numerical approaches it is possible to compute the resonant states, which are typically solutions to nonlinear eigenvalue problems (NLEVPs) arising from Maxwell's equations.The material dispersion described by the permittivity causes the nonlinearity of the eigenproblems.A multitude of numerical solution techniques are used for solving the NLEVPs, such as linearization, iterative projection methods and contour integral methods [4][5][6].
In nano-optics, linearization with physically derived auxiliary fields is a common approach [7][8][9][10][11][12].In this work, we report on an auxiliary field approach based on modeling the permittivity with rational functions.We implement the approach using an iterative projection method.Motivated by scatterometry applications, the numerical realization is applied to compute resonant states of a metallic line grating.

Auxiliary field approach for dispersive nano-optical systems
In the steady-state regime, the resonant states of nanooptical systems satisfy the time-harmonic Maxwell's equations in a source-free medium, given in the secondorder form by *Correspondence: burger@zib.de 1 Zuse Institute Berlin, Takustraße 7, 14195 Berlin, Germany 2 JCMwave GmbH, Bolivarallee 22, 14050 Berlin, Germany where E(r, ω) is the electric field.The permittivity tensor (r, ω), depending on the complex angular frequency ω and the position r, describes the material dispersion and the spatial distribution of materials.For optical frequencies, the permeability tensor μ(r, ω) typically equals the vacuum permeability μ 0 .Equation 1 becomes a non-Hermitian problem in the presence of open boundary conditions or lossy materials.
To obtain a numerical solution to Eq. (1), we apply the finite element method (FEM) [13,14].This discretization technique leads to an algebraic NLEVP of the form where A, B(ω) ∈ C n×n are the system matrices, ω ∈ C is an eigenvalue and u ∈ C n is the corresponding eigenvector.The problem is nonlinear through the eigenvaluedependence of the mass matrix B(ω), which is based on (r, ω).If the permittivity model (r, ω) is a rational function of the frequency with poles of order one, e.g., a Drude model [15] or a rational fit of measured material data, the matrix B(ω) has the form where B 0 , . . ., B N ∈ C n×n are matrices resulting from the partial fraction decomposition and ω 1 , . . ., ω N ∈ C are the poles of the rational function.Note that physical dispersion models have to satisfy Kramers-Kronig relations to ensure causality.
To compute eigenvalues ω and corresponding eigenvectors u, an implementation of the shift-and-invert Arnoldi method is applied [16].For this, the shifted eigenvalue ω = ω − σ , the shifted poles ωi = ω i − σ , i = 1, . . ., N, and the auxiliary fields are defined, where σ is the chosen shift.As the matrices B 1 , . . ., B N have only non-zero entries for degrees of freedom of the discretization corresponding to the dispersive object, the auxiliary fields u 1 , . . ., u N can be restricted to this subset.However, for the sake of a simpler notation, we define them on the entire domain.Using the auxiliary fields with Eq. ( 3) to reformulate Eq. ( 2) yields where I ∈ R n×n is the identity matrix.This is a linear eigenvalue problem of the form where Ã, B are augmented system matrices and ũ is an augmented field containing the original eigenvector u and the auxiliary fields u 0 , . . ., u N .The linear eigenvalue problem in Eq. ( 5) is solved by applying the Arnoldi method to assuming that ũ is suitably scaled.The Arnoldi iteration typically converges to the largest eigenvalue, i.e., to the smallest shifted eigenvalue ω = ω − σ .Thus, the eigenvalue ω of the NLEVP in Eq. ( 2) which is closest to the shift σ is obtained.Note that the auxiliary field approach increases the dimension of the eigenvalue problem with the number of poles of the rational function.
Remark For the computation of the Krylov subspace within the Arnoldi iteration for Eq. ( 6), the linear system Ã ũ = Bṽ is considered for the given input vector where ṽ ∈ C (N+2)n is an initial vector and f , f 0 , . . ., f N ∈ C n .The first n rows in Eq. ( 4) with the initial vector ṽ for the right-hand side lead to and substitution of the auxiliary fields yields Instead of solving the linear system Ã ũ = Bṽ to generate the Krylov subspace K m , the system Âu = ωf is solved yielding u and Eq. ( 7) is used to achieve u 0 , . . ., u N .This approach has the advantage that the matrix Â is equal to the matrix which is considered for solving Maxwell's equations in presence of a source.Such a scattering problem has the form where s(ω) is a source term.Setting ω = σ yields Â = A − σ 2 B(σ ) .Thus, the implementation of a scattering solver can also be used in the framework of solving eigenproblems.

Application to metallic grating
The presented approach is applied to a line grating consisting of gold struts surrounded by air.We revisit an experimentally realized setup supporting plasmonic resonances [17].This system has been recently numerically investigated [18].The geometry is sketched in Fig. 1.Grating structures are of interest in, e.g., scatterometry.It has been proposed to employ the resonant states of gratings for increasing the sensitivity in measurements of their spatial dimensions [19].We apply the auxiliary field approach using the FEM solver JCMsuite to compute the resonant state which corresponds to an absorption peak near the wavelength λ = 650 nm [18].For the relative permittivity of the gold grating, a one-pole Drude model Fig. 1 Sketch of a line grating consisting of gold surrounded by air.The structure is periodic in x direction and infinitely extended in y direction.The period is a = 482.5 nm, the rod width is w = 347.5 nm and the rod height is h = 130 nm is considered, where ω = 1.26e+16 s −1 is the plasma frequency and γ = 1.41e+14 s −1 is the damping coefficient.The permittivity is then given by (ω) = 0 r (ω), where 0 is the vacuum permittivity.The chosen shift is σ = 2πc/(650 nm), where c is the speed of light.Different finite element degrees p = 1, . . ., 6 and a fixed mesh containing about 1e+03 triangles are applied.Corners are a known issue considering systems containing metals.To deal with the occurring field singularities at the corners, refinements with a minimum edge length of about 0.016 nm are used.Bloch boundary conditions with a Bloch vector of [2π/(5a), 0, 0] enforce the periodicity in x direction.To realize the open boundary conditions in z direction, perfectly matched layers (PMLs) are used.Convergence of the PML method is ensured by applying an adaptive numerical realization of the PML method [13].The relative error of the eigenvalue ω is shown in Fig. 2, where the reference solution ω ref is the eigenvalue computed with p = 6.Convergence to the reference solution is observed.For the finite element degree p = 5, the eigenvalue ω = 2πc/(649.1397576+ 11.0601049i nm ± (6.2e−06 + 1.5e−06i nm)) is obtained.
In order to validate the results of the auxiliary field approach, eigenvalues are calculated using a fixed-point iteration.The same shift σ = 2πc/(650 nm) as before is used to initialize the mass matrix B(σ ).Equation 2 becomes linear and is solved with the shift-and-invert Arnoldi method.The resulting eigenvalue ω iter is then used to update B(ω iter ) and to repeat the procedure until ω iter does not change up to a chosen tolerance.In Table 1, the relative difference between the results from the auxiliary field approach, denoted by ω, and the results from the fixed-point iteration, denoted by ω iter , are shown.An abort condition for the fixed-point iteration with a tolerance of 1e−08 is chosen for the real and imaginary parts of ω iter .This leads to about 10 iterations.For all finite element degrees p = 1, . . ., 5, matching results for the two

Conclusions
We have reported an approach for computing eigensolutions to Maxwell's equations in dispersive media.Auxiliary fields are used to linearize the corresponding NLEVP.The resulting linear eigenvalue problem is then solved with the shift-and-invert Arnoldi method.The approach has been applied to a metallic line grating and the results for the eigenvalues have been validated by an implementation of a fixed-point iteration.

Fig. 2
Fig.2Convergence of the eigenvalue.Relative error of the eigenvalue computed with the auxiliary field approach with respect to the numerical resolution.The reference solution ω ref is computed with the finite element degree p = 6

Table 1
Comparison of eigenvalues computed with the auxiliary approach and with the fixed-point iteration, denoted by ω and ω iter , respectively