Euler-Lagrange equations for full topology optimization of the Q-factor in leaky cavities

TL;DR

This paper develops a mathematical framework using Euler-Lagrange equations for optimizing the Q-factor in lossy optical cavities, enabling better design of cavities with minimal decay rates.

Contribution

It introduces a novel set of nonlinear Maxwell equations derived from topology optimization principles, incorporating Pareto optimality and multi-parameter perturbation theory.

Findings

Derivation of Euler-Lagrange equations for cavity decay rate optimization

Formulation of nonlinear Maxwell systems with switching functions

Discussion of parallels with optimal control theory

Abstract

We derive Euler-Lagrange equations for the topology optimization of decay rate in 3-d lossy optical cavities. This leads to a new class of time-harmonic differential or integro-differential equations, which can be written as nonlinear Maxwell systems with switching functions of special types. Our approach is based on the notion of Pareto optimal frontier and on the multi-parameter perturbation theory for eigenfrequencies. Parallels with optimal control theory are discussed.

Equations

$$\mathrm{i}\nabla \times \mathbf{H}(\mathbf{x})$$

$$-\mathrm{i}\nabla \times \mathbf{E}(\mathbf{x})$$

$$0<{ϵ}_{-}\le {\epsilon }_{\mathrm{r}}(\mathbf{x})\le {ϵ}_{+},\mathbf{x}\in D_{\mathrm{opt}},$$

$$\begin{gathered} \text{find feasible permittivity functions εr} \\ \text{ that generate ωmin(α) as their eigenfrequency.} \end{gathered}$$

$$\frac{|{\Phi }_{\mathbf{E}}(\mathbf{x})|}{{ϵ}_0{\mu }_0}\nabla \times \nabla \times \mathbf{E}(\mathbf{x})={\omega }^2\left(\frac{{ϵ}_{+}+{ϵ}_{-}}{2}|{\Phi }_{\mathbf{E}}(\mathbf{x})|+\frac{{ϵ}_{+}-{ϵ}_{-}}{2}{\Phi }_{\mathbf{E}}(\mathbf{x})\right)\mathbf{E}(\mathbf{x}),$$

$${\Phi }_{\mathbf{E}}(\mathbf{x})=\Im (\mathbf{E}(\mathbf{x})\cdot \mathbf{E}(\mathbf{x})),\text{where E⋅E=E⊤E=∑j=13(Ej)2}$$

$${\Phi }_{\mathbf{E}}(x)={\int }_{c_{-}}^{c_{+}}\Im (\mathbf{E}(\mathbf{x})\cdot \mathbf{E}(\mathbf{x}))\mathrm{d}{x}_3,$$

$$\mathbf{n}\times \mathbf{E}-\mathcal{Z}{\mathbf{H}}_{\tau }=0,\text{mbox}for\mathbf{x}\in \partial D,$$

$${ϵ}_{-}\le {\epsilon }_{\mathrm{r}}(\mathbf{x})\le {ϵ}_{+}$$

$$\omega (\zeta )$$

$$\text{where }{\mathcal{C}}_1(p)$$

$$A=\{C(1/\epsilon -1/{\epsilon }_{\min }):\epsilon \in F,C\ge 0\}$$

$$T=\{{\mathcal{C}}_1({\zeta }_1{p}_1+{\zeta }_2{p}_2):{\zeta }_1>0,{\zeta }_2>0,{\zeta }_1+{\zeta }_2<1\},$$

$$\text{ {∫Dpεmin2E⋅E dx : p∈A} is a subset of the upper closed half-plane C+}.$$

$$\mathrm{i}|\Im (\mathbf{E}\cdot \mathbf{E})|\nabla \times \mathbf{H}={\omega }_0{ϵ}_0\left(\frac{{ϵ}_{+}+{ϵ}_{-}}{2}|\Im (\mathbf{E}\cdot \mathbf{E})|+\frac{{ϵ}_{+}-{ϵ}_{-}}{2}\Im (\mathbf{E}\cdot \mathbf{E})\right)\mathbf{E}$$

$$\mathrm{i}\nabla \times \mathbf{H}(\mathbf{x})-\mathrm{i}\sigma (\mathbf{x})\mathbf{E}(\mathbf{x})$$

$$-\mathrm{i}\nabla \times \mathbf{E}(\mathbf{x})$$

$${\mathcal{M}}_{{\epsilon }^{-1},{\mu }^{-1},\sigma }\left(\begin{array}{c}\mathbf{E}\\ \mathbf{H}\end{array}\right)=\left(\begin{array}{cc}-\frac{\mathrm{i}\sigma (\mathbf{x})}{\epsilon (\mathbf{x}){ϵ}_0}& \frac{\mathrm{i}}{\epsilon (\mathbf{x}){ϵ}_0}\nabla \times \\ -\frac{\mathrm{i}}{\mu (\mathbf{x}){\mu }_0}\nabla \times & 0\end{array}\right)\left(\begin{array}{c}\mathbf{E}\\ \mathbf{H}\end{array}\right)$$

$$\begin{gathered} ⟨\Psi ,\Psi {⟩}_{{\epsilon }_{\mathrm{r}}}={\int }_D\left({ϵ}_0⟨{\epsilon }_{\mathrm{r}}(\mathbf{x})\mathbf{E}(\mathbf{x}),\mathbf{E}(\mathbf{x}){⟩}_{{\mathbb{C}}^3}+{\mu }_0⟨\mathbf{H}(\mathbf{x}),\mathbf{H}(\mathbf{x}){⟩}_{{\mathbb{C}}^3}\right)\mathrm{d}\mathbf{x}= \\ ={\int }_D\left({ϵ}_0{\epsilon }_{\mathrm{r}}(\mathbf{x})\mathbf{E}(\mathbf{x})\cdot \overline{\mathbf{E}(\mathbf{x})}+{\mu }_0\mathbf{H}(\mathbf{x})\cdot \overline{\mathbf{H}(\mathbf{x})}\right)\mathrm{d}\mathbf{x}, \end{gathered}$$

$${\mathbb{M}}_0^{\ast }{\Psi }_{\star }=\overline{{\omega }_0}{\Psi }_{\star }.$$

$${\epsilon }_{\mathrm{r}}(\mathbf{x},\zeta )={\left({\epsilon }_{\mathrm{r}}(\mathbf{x}{)}^{-1}+\zeta p(\mathbf{x})\right)}^{-1}$$

$$\mathbb{M}(\zeta )\Psi (\zeta )=\omega (\zeta )\Psi (\zeta ),\Psi (\zeta )={\Psi }_0+\zeta {\Psi }_1+o(\zeta ).$$

$$({\mathbb{M}}_0+\zeta {\mathbb{M}}_1)({\Psi }_0+\zeta {\Psi }_1+o(\zeta ))=({\omega }_0+\zeta {\mathcal{C}}_1+o(\zeta ))({\Psi }_0+\zeta {\Psi }_1+o(\zeta )).$$

$${\mathcal{C}}_1{\Psi }_0$$

$${\mathcal{C}}_1⟨{\Psi }_0,{\Psi }_{\star }{⟩}_{{\epsilon }_{\mathrm{r}}}$$

$$0\ne ⟨{\Psi }_0,{\Psi }_{\star }{⟩}_{{\epsilon }_{\mathrm{r}}}={\int }_D\left({ϵ}_0{\epsilon }_{\mathrm{r}}{\mathbf{E}}^0\cdot {\mathbf{E}}^0-{\mu }_0{\mathbf{H}}^0\cdot {\mathbf{H}}^0\right)\mathrm{d}\mathbf{x}$$

$${\mathcal{C}}_1=\frac{⟨{\mathcal{M}}_{p,0,0}{\Psi }_0,{\Psi }_{\star }{⟩}_{{\epsilon }_{\mathrm{r}}}}{⟨{\Psi }_0,{\Psi }_{\star }{⟩}_{{\epsilon }_{\mathrm{r}}}}=\frac{{\int }_D⟨{\epsilon }_{\mathrm{r}}\mathrm{i}p{ϵ}_0^{-1}\nabla \times {\mathbf{H}}^0,\overline{{\mathbf{E}}^0}{⟩}_{{\mathbb{C}}^3}\mathrm{d}\mathbf{x}}{⟨{\Psi }_0,{\Psi }_{\star }{⟩}_{{\epsilon }_{\mathrm{r}}}}$$

$$=\frac{{\int }_{D}p⟨{\epsilon }_{\mathrm{r}}({\omega }_0{\epsilon }_{\mathrm{r}}+\mathrm{i}\sigma {ϵ}_0^{-1}){\mathbf{E}}^0,\overline{{\mathbf{E}}^0}{⟩}_{{\mathbb{C}}^3}\mathrm{d}\mathbf{x}}{⟨{\Psi }_0,{\Psi }_{\star }{⟩}_{{\epsilon }_{\mathrm{r}}}}=\frac{{\int }_{D}p({\omega }_0{\epsilon }_{\mathrm{r}}^2+\mathrm{i}\sigma {\epsilon }_{\mathrm{r}}{ϵ}_0^{-1})({\mathbf{E}}^0\cdot {\mathbf{E}}^0)\mathrm{d}\mathbf{x}}{⟨{\Psi }_0,{\Psi }_{\star }{⟩}_{{\epsilon }_{\mathrm{r}}}}.$$

$${\mathcal{C}}_1(p)={\omega }_0\frac{{\int }_{D}p(\mathbf{x}){\epsilon }_{\mathrm{r}}^2(\mathbf{x}){\sum }_{j=1}^3(E_j^0(\mathbf{x}){)}^2\mathrm{d}\mathbf{x}}{{\int }_D\left({ϵ}_0{\epsilon }_{\mathrm{r}}(\mathbf{x}){\sum }_{j=1}^3(E_j^0(\mathbf{x}){)}^2-{\mu }_0{\sum }_{j=1}^3(H_j^0(\mathbf{x}){)}^2\right)\mathrm{d}\mathbf{x}}.$$

$$\sigma (\mathbf{x})=0$$

$${\epsilon }_{\pm }(\mathbf{x})={ϵ}_{\pm }$$

$$\frac{|{\Phi }_{\mathbf{E}}|}{{ϵ}_0{\mu }_0}\nabla \times \nabla \times \mathbf{E}-\frac{\mathrm{i}\omega \sigma (\mathbf{x})}{{ϵ}_0}\mathbf{E}={\omega }^2\left(\frac{{\epsilon }_{+}(\mathbf{x})+{\epsilon }_{-}(\mathbf{x})}{2}|{\Phi }_{\mathbf{E}}|+\frac{{\epsilon }_{+}(\mathbf{x})-{\epsilon }_{-}(\mathbf{x})}{2}{\Phi }_{\mathbf{E}}\right)\mathbf{E}.$$

Full text

Matthias Eller${}^{\text{ a}}$ and Illya M. Karabash${}^{\text{ b,c}}$

${}^{\text{a}}$ Department of Mathematics and Statistics, Georgetown University, Washington, DC 20057

${}^{\text{b}}$ Mathematical Institute, Rheinische Friedrich-Wilhelms Universität Bonn, Endenicher Allee 60, D-53115 Bonn, Germany

${}^{\text{c}}$ Institute of Applied Mathematics and Mechanics of NAS of Ukraine, Dobrovolskogo st. 1, Slovyans’k 84100, Ukraine

E-mails: mme4@georgetown.edu, i.m.karabash@gmail.com

OCIS codes: (120.4570) Optical design of instruments; (140.3945) Microcavities; (230.5750) Resonators

MSC-classes: 35Q61, 35B34, 78M50, 49R05, 58E17, 90C29, 47B44, 47J10, 47A55

Keywords: spectral optimization, eigenvalue optimization, Maxwell equations, Pareto optimization, dissipation frequencies, photonic crystal, high-Q optical resonator, resonance optimization, qausi-normal-eigenvalue, quasi-normal level, multi-parameter perturbations

1. 1 Introduction

2. 1.1 Eigenproblem and Pareto optimization settings

3. 1.2 Main results and methods of the paper

4. 2 Euler-Lagrange eigenproblem in the case of impedance boundary conditions

5. 3 Model with absorbing conductivity layers

6. 3.1 One-parameter perturbations of eigenfrequencies

7. 3.2 EL-equations with absorbing layers and local optimizers

8. 3.3 Optimization of 2d and 1d photonic crystals embedded into a 3d outer medium

9. 4 Conclusions and discussion

10. A Appendix: 2-parameter perturbations of eigenvalues

1 Introduction

Optical cavities (or resonators) with high quality factor (Q-factor) are important components in contemporary optical engineering with a number of applications [40, 37, 38, 34] ranging from quantum electrodynamics [13, 41] to bio-imaging [36]. Depending on a particular purpose, various designs are employed including Fabry-Pérot cavities composed of mirrors or micromirrors of special shape [13, 34, 24], toroidal resonators supporting whispering gallery modes [40, 37], and photonic crystal cavities based on photonic bandgap effect [2, 32, 38].

Many applications and, in particular, quantum computing, require strong localization of light in a small volume. Since radiation loss is, roughly speaking, in inverse proportion to the resonator size, the task of fabrication of microresonators with high Q-factors is very nontrivial. A substantial attention was attracted by 2-d and 1-d photonic crystal designs because the parallel development of nanotechnology allowed researchers to fabricate sophisticatedly tailored cavities [2, 32, 9] and apply numerical methods of structural (or topology) optimization [29, 7, 28, 30, 4].

The numerical topology optimization may be based either on parametrization of certain class of designs by a handful of variables and running over them optimization algorithms [2, 32, 9, 30], or, in the case of the “full topology optimization”, on the consideration of every “pixel” of a discretized version of the resonator as a degree of freedom [28].

While microcavities with very high Q-factor have been fabricated, the ratio between predicted by simulations and actually observed Q-factors reaches one or even two orders of magnitude (see [9] and references in [30]). This essential difference was attributed to the neglection of absorption in simulations and to fabrication errors [30].

While further progress may be achieved by the study of robustness of high-Q designs and by optimization in presence of uncertainties (see [9]), an essential obstacle is that the optimization of Q-factor in 3-d optical cavities is not understood well enough from mathematical point of view.

The aims of this note are theoretical study of the full topology optimization for eigenmode Q-factors in 3-d optical cavities and the derivation of an analogue of the Euler-Lagrange equation for this problem.

While the employed cost functions appear to be standard and are connected with exponential loss rate (and Q-factor), they are used in the settings of Pareto optimization [6, 19, 3, 22], which is more common in mathematical economics, than in Physics. These settings resolve the difficulties with ill-posedness of the Q-factor maximization (the issues with nonexistence of optimizers were discussed in detail in [28]). Our approach shifts the difficulty to the sensitivity analysis, which requires in the case of Pareto optimization a more advanced multi-parameter perturbation theory [20].

1.1 Eigenproblem and Pareto optimization settings

We employ the eigenproblem-based formulation, which involves the time-harmonic Maxwell system

[TABLE]

equipped with lossy boundary conditions responsible to the leakage of energy. The leakage leads to the appearance of eigenfrequencies $\omega$ in the lower complex half-plane $\mathbb{C}_{-}:=\{z\in\mathbb{C}:\Im z<0\}$. The Silver-Müller radiation condition, absorbing layers, or impedance boundary conditions can be used for the modeling [5, 28, 31] of the dissipation effects.

In the optimization process, the spatially dependent (relative electric) permittivity $\varepsilon_{\mathrm{r}}(\mathbf{x})$ is modified inside a fixed bounded domain $D_{\mathrm{opt}}$ assuming that the constraints

[TABLE]

are satisfied. Here $\epsilon_{-}$ and $\epsilon_{+}$ are minimal and maximal permittivities of the materials available for the fabrication, for example, $\epsilon_{-}=1$ (air) and $\epsilon_{+}=11.9716$ (silicon). The permittivity functions $\varepsilon_{\mathrm{r}}(\cdot)$ satisfying (1.3) will be called feasible.

The problems of the maximization of the Q-factor and the minimization of the loss rate have been intensively investigated experimentally [2, 25], theoretically [41, 19, 28, 33, 21, 4, 22], and numerically [29, 16, 18, 32, 33, 28, 30, 10]. The quality-factor associated with a particular eigenmode $[\mathbf{E}\ \mathbf{H}]$ of (1.1), (1.2) can be defined as $Q=\frac{|\Re\omega|}{-2\Im\omega}$, where the real part $\alpha=\Re\omega$ of $\omega$ is the real-frequency of the corresponding eigenoscillations $\exp(-\mathrm{i}\omega t)[\mathbf{E}(x)\ \mathbf{H}(x)]$, and the negative of imaginary part $\Gamma(\omega)=-\Im\omega$ is the exponential loss rate. In the case of radiation conditions, $\Gamma(\omega)$ is the half-bandwidth of a resonance $\omega$.

The theoretical study of the maximization of the Q-factor requires a careful choice of an objective and of optimization constraints because many “straightforward” formulations of the problem are ill-posed [28]. The sequence of reformulations suggested in [28] to resolve this difficulty employs a number of ideas including the minimization of the volume of cavity for a given quality factor, the single complex-frequency formulation, and the frequency-averaged local density of states as objectives. In the somewhat different context of a quantum cavity, the issue with nonexistence of physically interesting minimizers was noticed already in the pioneering paper [14], where additional constraints $\alpha_{1}\leq\Re\omega\leq\alpha_{2}$ were suggested as one of possible recipes to avoid this difficulty.

The goal of present paper is to return to the simple loss rate objective $\Gamma(\omega)$, but to impose extreme versions of the real-frequency and cavity volume constraints with the aim to resolve the ill-posedness issue. That is, we assume that both a certain real-frequency $\Re\omega=\alpha$ and the region $D_{\mathrm{opt}}$, where the optimization is performed, are fixed. We take the loss rate $\Gamma(\omega)=-\Im\omega$ as the cost function. So, for a given $\alpha>0$, $\Gamma_{\min}(\alpha)$ is defined as the infimum of all possible $\Gamma$ such that $\omega=\alpha-\mathrm{i}\Gamma$ is an achievable complex-frequency, i.e., an eigenfrequency generated by a certain feasible function $\varepsilon_{\mathrm{r}}$.

Denoting the best possible value of eigenfrequency by $\omega_{\min}(\alpha)=\alpha-\mathrm{i}\Gamma_{\min}(\alpha)$, we can formulate the optimization problem in the following way:

[TABLE]

We call such $\varepsilon_{\mathrm{r}}(\cdot)$ optimal. In other words, we are interested in the maximal value of Q-factor for a prescribed value $\alpha$ of real-frequency $\Re\omega$.

For the 1-d equation $y^{\prime\prime}(s)=-\varepsilon_{\mathrm{r}}(s)\epsilon_{0}\mu_{0}\omega^{2}y(s)$, where $s=x_{3}$ and where $y=E_{2}$ describes a TEM-mode in a layered cavity, the similar formulation was theoretically and computationally studied in [19, 21, 22]. Generally, it is difficult to guarantee that the real parts $\Re\omega$ of achievable $\omega$ cover the whole positive semi-line $\mathbb{R}_{+}=\{\alpha\in\mathbb{R}:\alpha>0\}$. Such real parts $\Re\omega$ will be called achievable and their set is denoted by $\operatorname{dom}\Gamma_{\min}$. The 1d results of [21, 22] suggest that all high enough real parts $\Re\omega$ are achievable.

The above theoretical formulation is natural from the point of view of the Pareto optimization [6]. Indeed, the curve $\{\alpha-\mathrm{i}\Gamma_{\min}(\alpha):\alpha\in\operatorname{dom}\Gamma_{\min}\}$ is the Pareto optimal frontier for achievable complex-eigenfrequencies $\omega$ [19, 20].

1.2 Main results and methods of the paper

Our main goal is to derive the nonlinear Maxwell-type equation

[TABLE]

where $\Phi_{\mathbf{E}}(\mathbf{x})$ plays the same role as switching functions do in optimal control theory. In the case of 3d full topology Pareto optimization the switching function can be written in the form

[TABLE]

and $\mathbf{E}=[E_{1}\ E_{2}\ E_{3}]=(E_{1}\ E_{2}\ E_{3})^{\top}$ is perceived as a $\mathbb{C}^{3}$-column vector. This equation is fulfilled inside the optimization region $D_{\mathrm{opt}}$ for the $\mathbf{E}$-component of one of eigenmodes $[\mathbf{E}\ \mathbf{H}]$ of (1.1), (1.2) corresponding to the optimal $\omega=\omega_{\min}(\alpha)$.

We will say that (1.4) equipped with a leaky boundary condition is an Euler-Lagrange eigenproblem (EL-eigenproblem). In Section 2, the EL-eigenproblem is derived in the simplest possible settings, where $\omega$ is assumed to be a simple complex-eigenfrequency of the Maxwell system with a boundary condition of impedance type. Then, in Section 3, we show that the form of the EL-eigenproblem does not substantially depend on a particular choice of the type of lossy boundary conditions (while optimal eigenmodes and permittivity functions $\varepsilon_{\mathrm{r}}$ obviously do depend). To this end a more realistic model of a leaky cavity with an absorbing layer in the outer domain is used. In Section 4 we consider the connection between optimizers $\varepsilon_{\mathrm{r}}(\mathbf{x})$ and the eigenfields $\mathbf{E}(\mathbf{x})$ of the EL-eigenproblems, and discuss some of related mathematical questions, in particular, the parallels with optimal control theory.

Since 3d photonic crystals with sophisticatedly tailored design are difficult in fabrication, the most essential efforts in the experimental and numerical studies were concentrated on structures with 2d and 1d geometries [2, 16, 32, 19, 21, 10, 22].

Photonic crystals with 2d (or 1d) structure described in $D_{\mathrm{opt}}$ by $\varepsilon_{\mathrm{r}}(\mathbf{x})$ depending only on 2 components (or a 1 component) of the $\mathbf{x}$-vector, but surrounded by a 3d outer medium, are considered in Section 3.3. Using the same tools as in the 3d-geometry case, we show that corresponding EL-equations have the form (1.4) with different switching functions $\Phi_{\mathbf{E}}$ involving integrals over cross-sections of the optimization domain $D_{\mathrm{opt}}$, which is assumed now to be cylindrical. For example, in the case of optimization of a 2d photonic crystal we obtain the switching function

[TABLE]

where $\{\mathbf{x}=(x_{1},x_{2},x_{3})\in\mathbb{R}^{3}\ :\ c_{-}<x_{3}<c_{+}\}$ is the slab containing the optimization domain $D_{\mathrm{opt}}$. So, in the case of 1d and 2d geometries, equation (1.4) becomes of an integro-differential type.

The main tool for the derivation of the EL-eigenproblem is an analytical lemma about the sets covered by 2-parameter perturbations of zeros of analytic functions, which is a particular case of a results of [20] and is given in Appendix for a convenience of the reader.

To apply the 2-parameter lemma, we need as a preliminary step a formula for the first correction term of an eigenfrequency of (1.1)-(1.2) under 1-parameter perturbations of $\varepsilon_{\mathrm{r}}(\mathbf{x})$. This perturbation formula is derived in Section 3.1. While there is an abundance of studies on perturbation expansions for various types of Maxwell equations (see [4, 17, 19, 21, 26] and references therein), we were not able to find in the literature a corresponding perturbation result for a lossy cavity. For resonances corresponding to outgoing conditions, a somewhat less explicit formula can be found in [26], where the degenerate case is also included. On the other side, the formula of Section 3.1 can be seen as an extension to the case of 3d Maxwell equations of the perturbation results on the 3d Schrödinger equation [14] and on the 1d string case [19, 21]. In 1d settings, the leading correction term of Puiseux series for resonances of arbitrary multiplicity was studied analytically in [19, 21].

2 Euler-Lagrange eigenproblem in the case of

impedance boundary conditions

Let us consider a bounded 3d domain $D$ (i.e., a bounded open connected subset of $\mathbb{R}^{3}$) with a smooth boundary $\partial D$ and impose on $\partial D$ the impedance boundary condition of the simplest form

[TABLE]

where $\mathbf{n}$ is the unit outward normal vector to $\partial D$, $\mathbf{H}_{\tau}=(\mathbf{n}\times\mathbf{H})\times\mathbf{n}$ is the tangential component of the field $\mathbf{H}$, and the impedance coefficient $\mathcal{Z}(\mathbf{x})$ is a scalar uniformly positive $L^{\infty}(\partial D)$-function (uniformly positive means that there exists a constant $C>0$ such that $\mathcal{Z}(\mathbf{x})\geq C$ on the domain of its definition).

An optimization domain $D_{\mathrm{opt}}$ is a 3d subdomain of $D$ that represents the region where the scalar relative permittivity function $\varepsilon_{\mathrm{r}}(\cdot)$ can be changed in the process of optimization. The structure of the cavity in the outer region $D_{\mathrm{out}}=D\setminus D_{\mathrm{opt}}$ is assumed to be known, unchangeable during the optimization, however, not necessarily homogeneous. We assume that this outer structure is described by a uniformly positive permittivity function $\varepsilon_{\mathrm{out}}(\mathbf{x})$, $\mathbf{x}\in D_{\mathrm{out}}$.

The family of all $L^{\infty}(D)$-functions $\varepsilon_{r}$ such that

[TABLE]

will be called the family of feasible permittivities and denoted by $F$.

Assume that $\alpha>0$, that the optimal complex-eigenfrequency $\omega_{0}=\omega_{\min}(\alpha)$ (see Section 1.1) is achieved for a certain feasible permittivity function $\varepsilon_{r}(\mathbf{x})=\varepsilon_{\min}(\mathbf{x})\in F$, and that for this $\varepsilon_{\min}$ the eigenfrequency $\omega_{0}$ is simple (i.e. of algebraic multiplicity 1). Let $\Psi_{0}(\mathbf{x})=[\mathbf{E}^{0}(\mathbf{x})\ \mathbf{H}^{0}(\mathbf{x})]$, $\mathbf{x}\in D$, be a corresponding eigenfield, which is unique up to a multiplication on a complex constant.

We consider linear perturbations $\zeta p+\frac{1}{\varepsilon_{\min}}$ in the direction $p(\mathbf{x})$, $\mathbf{x}\in D$, of the inverse $1/\varepsilon_{\min}(\mathbf{x})$ to the permittivity function. A direction $p$, which is also an $L^{\infty}(D)$-function, is called admissible if $\varepsilon_{r}(\mathbf{x},\zeta)=\left((\varepsilon_{\min}(\mathbf{x}))^{-1}+\zeta p(\mathbf{x})\right)^{-1}$ is feasible for small enough positive numbers $\zeta$. The perturbed eigenfrequency $\omega(\zeta)$ corresponding to $\varepsilon_{r}(\mathbf{x},\zeta)$ is an analytic complex valued function in $\zeta$ for complex $\zeta$ close enough to [math] and satisfies

[TABLE]

and $\mathcal{K}_{1}=\left(\int_{D}\left(\varepsilon_{\min}\epsilon_{0}\mathbf{E}^{0}\cdot\mathbf{E}^{0}-\mu_{0}\mathbf{H}^{0}\cdot\mathbf{H}^{0}\right)\mathrm{d}\mathbf{x}\right)^{-1}\neq 0$ (see Section 3.1, where this formula is analytically derived in more general settings).

Note that any admissible direction $p$ satisfies $p(\mathbf{x})=0$ in $D_{\mathrm{out}}$, and so, $p$ can be considered as $L^{\infty}(D_{\mathrm{opt}})$-function. The first correction term $\mathcal{C}_{1}$ depends also on $\omega_{0}$, $\varepsilon_{\min}$, and $\Psi_{0}$, but since they are fixed this dependence is not important for our needs. Then the mapping $p\mapsto\mathcal{C}_{1}(p)$ can be considered as a complex-valued linear bounded functional on $L^{\infty}(D_{\mathrm{opt}})$.

For the convenience in the use of Convex Analysis techniques, we identify the complex plane $\mathbb{C}$ with the 2d real plane $\mathbb{R}^{2}$. So $p\mapsto\mathcal{C}_{1}(p)$ will be simultaneously considered as an $\mathbb{R}^{2}$-valued and as a $\mathbb{C}$-valued real-linear mapping defined on $L^{\infty}(D_{\mathrm{opt}})$.

Since the family $F$ of feasible permittivities is closed and convex in $L^{\infty}(D)$, the set of admissible directions

[TABLE]

and its image $\mathcal{C}_{1}[A]=\{\mathcal{C}_{1}(p)\ :\ p\in A\}$ under the mapping $p\mapsto\mathcal{C}_{1}(p)$ are closed convex cones in $L^{\infty}(D_{\mathrm{opt}})$ and $\mathbb{R}^{2}$, respectively.

Obviously, the origin [math] belongs to the image $\mathcal{C}_{1}[A]$ because one can take $p\equiv 0$ in $D_{\mathrm{opt}}$. Moreover,

[TABLE]

Otherwise, one can take two admissible directions $p_{1}(\mathbf{x})$ and $p_{2}(\mathbf{x})$ so that the open triangular domain

[TABLE]

which is the interior of the convex hull of [math], $\mathcal{C}_{1}(p_{1})$, and $\mathcal{C}_{1}(p_{2})$, contains a certain complex segment $\mathrm{i}(0,\delta)=\{\mathrm{i}\tau\ :\ 0<\tau<\delta\}$ of the positive imaginary semi-axis $\mathrm{i}\mathbb{R}_{+}$. Then the Maclaurin series (2.3) and the homotopy arguments (see Appendix A) easily lead to the conclusion that $\omega_{0}-\mathrm{i}\delta_{1}$ is an achievable eigenfrequency for a certain $\delta_{1}>0$, and so, to the contradiction with the assumption that $\omega_{0}$ has a minimal possible loss rate $\Gamma(\omega)$ among all achievable $\omega$ on the line $\{\Re\omega=\alpha\}$.

Since $\mathcal{C}_{1}[A]$ is a closed convex cone having [math] on its boundary, there exists $\theta\in\mathbb{R}$ such that the closed half-plane $e^{\mathrm{i}\theta}\overline{\mathbb{C}}_{+}=\{e^{\mathrm{i}\theta}z:\Im z\geq 0\}$ contains $\mathcal{C}_{1}[A]$.

Summarizing, we see from the form (2.4) of $\mathcal{C}_{1}(p)$ that replacing, if necessary, $\Psi_{0}$ by its multiple $[\mathbf{E}\ \mathbf{H}]=e^{\mathrm{i}\theta_{1}}\Psi_{0}$ with a certain phase $\theta_{1}\in\mathbb{R}$ one can ensure that

[TABLE]

This and the constraints $\epsilon_{-}\leq\varepsilon_{\mathrm{r}}(\mathbf{x})\leq\epsilon_{+}$ of the feasible family $F$ give for $\mathbf{x}\in D_{\mathrm{opt}}$ that $\varepsilon_{\min}(\mathbf{x})=\epsilon_{\pm}$ whenever $\pm\Im(\mathbf{E}(\mathbf{x})\cdot\mathbf{E}(\mathbf{x}))>0$.

Combining the last observation with the Maxwell system, one gets for $\mathbf{x}\in D_{\mathrm{opt}}$

[TABLE]

and, in turn, the EL-eigenproblem (1.4) with $\Phi_{\mathbf{E}}=\Im(\mathbf{E}\cdot\mathbf{E})$ and $\omega=\omega_{0}$.

Remark 2.1*.*

One can choose other forms of the switching function. For example, $\Phi_{E}=\Re(\mathbf{E}\cdot\mathbf{E})$ leads to essentially equivalent EL-eigenproblem with a solution $\mathbf{E}$ replaced by $e^{-\mathrm{i}\pi/4}\mathbf{E}$. Note that the EL-equation is linear with respect to (w.r.t) multiplication of $\mathbf{E}$ on real constants $C$, but is nonlinear w.r.t. multiplication on complex constants $C\in\mathbb{C}\setminus\mathbb{R}$.

3 Model with absorbing conductivity layers

A cavity in Optical Engineering is usually a part of a bigger complex structure. The leakage of energy occurs rather due to various dissipative effects in the medium near the cavity, than due to escape of waves to the far surrounding modeled by an idealized infinite vacuum. One can make the model of a lossy cavity more realistic and adjustable to various specific setting by the introduction of artificial absorbing region.

On the other hand, absorption layers are intensively employed in simulations because of their convenience from the point of view of numerical approximations [31].

A simplest way to model the absorption is to introduce a spatially varying positive conductivity $\sigma(\mathbf{x})$. This approach also allows one to take into account possible metallic and semiconductor structures in the surrounding region $D_{\mathrm{out}}$.

For a time being, let us impose minimal assumptions on $\sigma(\mathbf{x})$ supposing that it is a nonnegative $L^{\infty}(D)$-function. This model leads to the time-harmonic Maxwell system

[TABLE]

As before we assume that the impedance boundary condition (2.1) is imposed on $\partial D$.

We will denote the corresponding pseudo-Hamiltonian [12] by $\mathcal{M}_{\varepsilon_{\mathrm{r}}^{-1},1,\sigma}$, where

[TABLE]

is defined analogously to classical conservative Hamiltonians $\mathcal{M}_{\varepsilon^{-1},\mu^{-1},0}$ of [8, 17] with the only difference that, instead of conservative boundary conditions on $\partial D$, the dissipative impedance boundary condition participates in the definition of the domain of the operator [11, 27].

3.1 One-parameter perturbations of eigenfrequencies

In this subsection we derive the first correction term for perturbations of a simple isolated eigenvalue $\omega_{0}$ of the linear operator $\mathbb{M}_{0}=\mathcal{M}_{\varepsilon_{\mathrm{r}}^{-1},1,\sigma}$. Let $\Psi_{0}=[\mathbf{E}^{0}\ \mathbf{H}^{0}]$ be an associated eigenfield.

In the case of lossy cavity, $\mathbb{M}_{0}$ is not Hermitian and $\omega_{0}$ is not necessarily real. Therefore an additional step is needed in the derivation of the first correction term. The resulting formula does not exactly coincide with the well-known first-order correction for a lossless cavity [17]. This additional step requires the use of the operator $\mathbb{M}_{0}^{*}$ adjoint to $\mathbb{M}_{0}$ w.r.t. the energy form

[TABLE]

where $\langle\cdot,\cdot\rangle_{\mathbb{C}^{3}}$ is a standard sesquilinear inner product in the 3d complex vector space $\mathbb{C}^{3}$, and $\overline{\mathbf{E}}=[\overline{E_{1}}\ \overline{E_{2}}\ \overline{E_{3}}]$, $\overline{\mathbf{H}}=[\overline{H_{1}}\ \overline{H_{2}}\ \overline{H_{3}}]$ consist of complex conjugates to the components of $\mathbf{E}$ and $\mathbf{H}$.

We use the observation that the field $\Psi_{\star}$ defined by $\Psi_{\star}=\begin{pmatrix}\overline{\mathbf{E}_{0}}\\ -\overline{\mathbf{H}_{0}}\end{pmatrix}$ is an adjoint eigenstate of $\mathbb{M}_{0}^{*}$ in the sense that

[TABLE]

To derive the first correction term to $\omega_{0}$ under the perturbations

[TABLE]

of the permittivity function, we introduce the operator-valued function $\mathbb{M}(\zeta)=\mathbb{M}_{0}+\zeta\mathbb{M}_{1},$ where $\mathbb{M}_{1}=\mathcal{M}_{p,0,0}$ and $\zeta$ is a complex number in a vicinity of [math].

Since $\mathbb{M}(\zeta)$ is analytic and of Kato’s type (A) [19, Sect. VII.2], the perturbed eigenvalue $\Omega(\zeta)$ is analytic near $\zeta=0$ and $\omega(\zeta)=\omega_{0}+\mathcal{C}_{1}\zeta+o(\zeta).$ Moreover, there exists an eigenstate $\Psi(\mathbf{x},\zeta)=\begin{pmatrix}\mathbf{E}(\mathbf{x},\zeta)\\ \mathbf{H}(\mathbf{x},\zeta)\end{pmatrix}$ analytically dependent on $\zeta$ and satisfying

[TABLE]

From $\mathbb{M}(\zeta)\Psi(\zeta)=\omega(\zeta)\Psi(\zeta)$, we see that

[TABLE]

and so

[TABLE]

Due to the assumption that $\omega_{0}$ is a simple eigenvalue, one has

[TABLE]

and

[TABLE]

We see that the first correction term $\frac{\mathrm{d}\omega(\xi)}{\mathrm{d}\xi}=\mathcal{C}_{1}(p)$ is a linear functional of the perturbation direction $p$.

We apply (3.4) under additional assumption that $p(\mathbf{x})\sigma(\mathbf{x})=0$ everywhere in $D$ (cf. (3.6)), which leads to a simpler formula

[TABLE]

Remark 3.1*.*

For the lossless cavity, (2.4) can be reduced to the standard formula [17] since $\omega$ is real and so the the real and imaginary parts of the eigenfields $\mathbf{E}^{0}$ and $\mathbf{H}^{0}$ are also eigenfields.

3.2 EL-equations with absorbing layers and local optimizers

The procedure of derivation of EL-equation in the presence of conductivity is similar to Section 2. However, we extend the EL-eigenproblem also to the case of local extrema.

Assume that the function $\varepsilon_{\mathrm{r}}$ belongs to the family of feasible permittivities $F$. That is, it is fixed in $D_{\mathrm{out}}$ and is a subject of optimization inside the optimization domain $D_{\mathrm{opt}}$ satisfying the constraints (2.2). Suppose that $\sigma(\cdot)$ is fixed during the optimization in whole $D$ and additionally

[TABLE]

i.e., the medium is lossless inside the optimization domain $D_{\mathrm{opt}}$. So the loss of energy happens in a subset of $D_{\mathrm{out}}$ where $\sigma(\mathbf{x})>0$ and on the boundary $\partial D$ due to the impedance boundary condition.

Let $\alpha>0$ and $\Re\omega_{0}=\alpha$. Assume that $\omega_{0}$ is a simple eigenvalue of the Maxwell system (3.1), (3.2), (2.1) with a certain $\varepsilon_{\mathrm{r}}\in F$.

To define local minimizers, let us introduce the family $F(\alpha)$ of pairs $\{\omega,\varepsilon_{\mathrm{r}}\}$ that satisfy the following conditions: (i) $\Re\omega=\alpha$, (ii) $\omega$ is an eigenvalue of $M_{\varepsilon_{\mathrm{r}}^{-1},1,\sigma}$, (iii) $\varepsilon_{\mathrm{r}}(\mathbf{x})$ is a feasible permittivity function. That is, $F(\alpha)$ is a family of feasible pairs $\{\omega,\varepsilon_{\mathrm{r}}\}$ under the additional constraint $\Re\omega=\alpha$, and can be considered as a subset of the normed space $\mathbb{C}\times L^{\infty}(D)$ with any reasonable norm $\|\cdot\|$ (e.g., $\|\{\omega,\varepsilon_{\mathrm{r}}\}\|=\max\{|\omega|,\|\varepsilon_{\mathrm{r}}\|_{L^{\infty}}\}$).

Assume that $\{\omega_{0},\varepsilon_{\mathrm{r}}\}$ is a local minimizer in the sense that $\Gamma(\omega_{0})\leq\Gamma(\omega)$ for all pairs $\{\omega,\varepsilon\}\in F(\alpha)$ in a ball of the space $\mathbb{C}\times L^{\infty}(D)$ centered in $\{\omega_{0},\varepsilon_{\mathrm{r}}\}$ with a certain positive radius. Then applying the perturbation formula (3.5) and the arguments of Section 2, one sees that there exists an eigenfield $[\mathbf{E}\ \mathbf{H}]$ of (3.1), (3.2), (2.1) that satisfies the EL-equation (2.5) in the optimization domain $D_{\mathrm{opt}}$.

The form of EL-equation in $D_{\mathrm{opt}}$ does not depend on the way how the leakage of energy is modeled in the outer region. However, the corresponding eigenfield $[\mathbf{E}\ \mathbf{H}]$ in $D_{\mathrm{opt}}$ does generally depend and so does the local optimal permittivity $\varepsilon_{\mathrm{r}}$. Assuming that in the surrounding region $D_{\mathrm{out}}$ there exists an absorption subdomain $D_{\mathrm{abs}}$, where $\sigma(x)>0$ (and that $D_{\mathrm{abs}}$ has a positive volume), one sees that the EL-eigenfield $[\mathbf{E}\ \mathbf{H}]$ satisfies in $D_{\mathrm{out}}$ a system different from that of Section 2. The glue-type condition between the parts of eigenfield in $D_{\mathrm{opt}}$ and $D_{\mathrm{out}}$ is encoded in the domain of definition of the pseudo-Hamiltonian $\mathcal{M}_{\varepsilon_{\mathrm{r}}^{-1},1,\sigma}$ (namely, in the requirements that the distributions $\nabla\times\mathbf{E}$ and $\nabla\times\mathbf{H}$ are $L^{2}(D)$-vector-functions).

To include this dependence on the absorption into the EL-system, let us introduce in the whole $D$ the constraint functions $\varepsilon_{\pm}(\cdot)$ defined by

[TABLE]

and write the EL-equation in the whole $D$ as

[TABLE]

with $\Phi_{\mathbf{E}}=\Im(\mathbf{E}\cdot\mathbf{E})$.

Remark 3.2*.*

If we drop the assumption that $\sigma(\mathbf{x})=0$ in $D_{\mathrm{opt}}$, but suppose instead that $\sigma(\mathbf{x})$, $\mathbf{x}\in D$, is fixed in the process of optimization of $\varepsilon_{\mathrm{r}}$, the resulting equation (3.7) is the same, but another switching function of the form $\Phi_{\mathbf{E}}(\mathbf{x})=\Im\left(\left(\varepsilon_{\mathrm{r}}^{2}+\frac{\mathrm{i}\sigma\varepsilon_{\mathrm{r}}}{\epsilon_{0}\omega}\right)(\mathbf{E}\cdot\mathbf{E})\right)$ has to be used. However, if one would like to take into account small conductivity effects inside $D_{\mathrm{opt}}$, a reasonable optimization problem should include also optimization of $\sigma(\mathbf{x})$, $\mathbf{x}\in D_{\mathrm{opt}}$, over a certain range of values. Such a model leads to switching effects not only for $\varepsilon_{\mathrm{r}}$, but also for $\sigma$. The resulting EL-equation is more technically involved and will be considered elsewhere.

Let $[\mathbf{E}\ \mathbf{H}]$ be a solution to the linear Maxwell system (3.1), (3.2), (2.1) with a certain feasible $\varepsilon_{\mathrm{r}}$ such that the eigenfield $\mathbf{E}$ satisfies additionally the EL-eigenproblem (3.7). We will say that such $\mathbf{E}(\mathbf{x})$ is an EL-eigenfield for a nonlinear eigenvalue $\omega$.

Remark 3.3*.*

Generally, there is no reason to expect the uniqueness of optimal $\varepsilon_{\mathrm{r}}$ (see [14, 3, 22] for the discussion of uniqueness in the 1d case and in the case of a quantum cavity) and so there is no reason to expect the uniqueness of an EL-eigenfield.

3.3 Optimization of 2d and 1d photonic crystals embedded into a 3d outer medium

In this section, the outer medium described by $\varepsilon_{\mathrm{out}}(\mathbf{x})$ and $\sigma(\mathbf{x})$ has as before a 3d geometry, as well as the outer domain $D_{\mathrm{out}}$.

The optimization domain of the cavity is taken to be a cylinder

[TABLE]

where $D_{2}$ is a bounded 2d domain and $c_{\pm}\in\mathbb{R}$ are certain constatnts so that $c_{-}<c_{+}$.

The internal structure of the cavity in $D_{\mathrm{opt}}$ is assumed to have either 2d, or 1d geometry. To model this, suppose that the permittivity function $\varepsilon_{\mathrm{r}}(\mathbf{x})$ in the optimization domain $D_{\mathrm{opt}}$ either depends only on $(x_{1},x_{2})$, or only on $x_{3}$. That is, $\varepsilon_{\mathrm{r}}(\mathbf{x})$, $\mathbf{x}\in D_{\mathrm{opt}}$, belongs

• •

either to the feasible family $F_{2}$ defined by the previous constraints (2.2) and the additional assumption that for each $(x_{1},x_{2})\in D_{2}$ we have $\varepsilon_{\mathrm{r}}(x_{1},x_{2},x_{3})=\varepsilon_{\mathrm{r}}(x_{1},x_{2},x^{\prime}_{3})$ for all $x_{3},x^{\prime}_{3}\in(c_{-},c_{+})$,

• •

or, to the feasible family $F_{1}$ defined by (2.2) and the assumption that for $x_{3}\in(c_{-},c_{+})$ we have $\varepsilon_{\mathrm{r}}(x_{1},x_{2},x_{3})=\varepsilon_{\mathrm{r}}(x^{\prime}_{1},x^{\prime}_{2},x_{3})$ for all points $(x_{1},x_{2}),(x^{\prime}_{1},x^{\prime}_{2})\in D_{2}$.

Suppose that $\sigma(\cdot)$ is fixed and $\sigma(\mathbf{x})=0$ for $\mathbf{x}\in D_{\mathrm{opt}}$.

Similarly to the previous section, we assume that the pair $\{\omega,\varepsilon_{\mathrm{r}}\}$ is minimizer or local minimizer over $F_{2}$ or over $F_{1}$ for $\Gamma(\omega)=-\Im\omega$ under the real part constraint $\Re\omega_{0}=\alpha>0$.

Then, the arguments of Section 2 show that there exists an eigenfield $[\mathbf{E}\ \mathbf{H}]$ of the linear Maxwell system (3.1), (3.2), (2.1) that satisfies in $D_{\mathrm{opt}}$ the EL-equation

[TABLE]

with the switching function $\Phi_{\mathbf{E}}$ of the form

• •

$\Phi_{\mathbf{E}}(\mathbf{x})=\int_{c_{-}}^{c_{+}}\Im(\mathbf{E}(\mathbf{x})\cdot\mathbf{E}(\mathbf{x}))\mathrm{d}x_{3}$ in the case of optimization over the family $F_{2}$ of 2d photonic crystals;

• •

$\Phi_{\mathbf{E}}(\mathbf{x})=\int_{D_{2}}\Im(\mathbf{E}(\mathbf{x})\cdot\mathbf{E}(\mathbf{x}))\mathrm{d}x_{1}\mathrm{d}x_{2}$ in the case of optimization over the family $F_{1}$ of 1d photonic crystals.

So the EL-equation has a non-local character of an integro-differential equation in $D_{\mathrm{opt}}$. Similarly to Section 3.2, it can be written as a nonlinear eigenvalue problem (3.7), (2.1) in the whole domain $D$.

4 Conclusions and discussion

We have derived new Euler-Lagrange equations for the Pareto decay rate minimization in 3-d, 2-d, and 1-d photonic crystals embedded into a structured 3-d outer medium with dissipative effects.

The efficiency of this approach to computations of the Pareto optimal frontier of optimal eigen-frequencies was up to now demonstrated only for the 1-d case of TEM-modes in idealized layered cavities [21, 22]. Application of the same scheme in 3-d settings require a deeper study of a series of new mathematical questions, which are briefly discussed below.

Let us consider the dependence of locally optimal $\varepsilon_{\mathrm{r}}$ on the EL-eigenfield $\mathbf{E}$ in more detail. Section 2 shows that, inside $D_{\mathrm{opt}}$, we have $\varepsilon_{\mathrm{r}}(\mathbf{x})=\epsilon_{\pm}$ whenever $\pm\Phi_{\mathbf{E}}(\mathbf{x})>0$. However, (2.5) does not impose any assumption on $\varepsilon_{\mathrm{r}}(\mathbf{x})$ in the set $S$ of points $\mathbf{x}\in D_{\mathrm{opt}}$ where $\Phi_{\mathbf{E}}(\mathbf{x})=0$. Such phenomenon is common for some optimal control problems, for example, for economic trading model and the moon landing problem [35]. Adapting optimal control terminology, it is natural to call $S$ the singular set for the switching function $\Phi_{\mathbf{E}}$ (for connections between 1d resonance optimization and optimal control see [22]).

Since $S$ is the zero-level set of the switching function, it is natural to expect that, generically, it is a 2d manifold and so has a zero volume (i.e., zero 3d Lebesgue measure). In such a case, the values of $\varepsilon_{\mathrm{r}}$ on $S$ can be discarded since any change of $L^{\infty}$-coefficient $\varepsilon_{\mathrm{r}}$ on a set of zero volume does not influence the solution of the Maxwell system. Then adjusting Optimal Control terminology, it is natural to say that optimal $\varepsilon_{\mathrm{r}}(\mathbf{x})$ has a bang-bang structure in the sense that it takes in $D_{\mathrm{opt}}$ only two extreme possible values $\epsilon_{\pm}$. From the Engineering point of view, this corresponds to the high contrast design, where air holes are sophisticatedly tailored in a silicon crystal.

However, it is difficult to rule out the existence of 3d settings where a switching function has a singular sets of positive volume. The main reason to expect that there exist pathological settings of such type is the fact that presently available unique continuation results require Lipschitz regularity of the coefficient $\varepsilon_{\mathrm{r}}$ [39, 15]. In the optimization problem under consideration $\varepsilon_{\mathrm{r}}$ has, generally, a lower regularity of the $L^{\infty}$-space.

Remark 4.1*.*

Since the dynamic electromagnetic field can be considered to be real-valued, the complexification of the fields is due solely to the time-harmonic setting. One observes that $\Phi_{\mathbf{E}}=\Im(\mathbf{E}\cdot\mathbf{E})=2\Re\mathbf{E}\cdot\Im\mathbf{E}$. So the equality $\Im(\mathbf{E}\cdot\mathbf{E})=0$ means that $\Re\mathbf{E}=[\Re E_{1}\ \Re E_{2}\ \Re E_{3}]$ and $\Im\mathbf{E}=[\Im E_{1}\ \Im E_{2}\ \Im E_{3}]$ are orthogonal on $S$.

Appendix A Appendix: 2-parameter perturbations of eigenvalues

Let us consider a function $Q(z;\zeta)=Q(z;\zeta_{1},\zeta_{2})$ of 3 complex variables $z$, $\zeta_{1}$, $\zeta_{2}$ assuming that this function is analytic in a polydisc $\mathbb{D}_{R}^{3}$, where $\mathbb{D}_{R}=\{\omega\in\mathbb{C}\ :\ |\omega|<R\}$. The variables $\zeta_{1}$ and $\zeta_{2}$ are united in the pair $\zeta=(\zeta_{1},\zeta_{2})\in\mathbb{C}^{2}$ and play the role of perturbation parameters. We are interested in the behavior of the set of zeros of the function $Q(\cdot;\zeta)$. That is, we are interested in the set $\Sigma_{Q}(\zeta)$ of complex numbers $z\in\mathbb{D}_{R}$ satisfying the equation $Q(z;\zeta)=0$, and in the local evolution of this set under small perturbations of $\zeta$.

We take real positive solid triangles $T_{\delta}:=\{\zeta\in\mathbb{D}_{R}^{2}\ :\ \zeta_{1},\zeta_{2}>0\ \text{ and }\ 0<\zeta_{1}+\zeta_{2}<\delta\}$ as families of feasible perturbation parameters. Then $\Sigma_{Q}[T_{\delta}]:=\bigcup_{\zeta\in T_{\delta}}\Sigma_{Q}(\zeta)$ is the set of $z$-zeros achievable over $T_{\delta}$.

Assume that

• •

$z_{0}=0$ is a simple zero of a function $Q(\cdot,0,0)$,

• •

$\eta_{j}:=-\frac{\partial Q(0,0,0)}{\partial\zeta_{j}}\left(\frac{\partial Q(0,0,0)}{\partial z}\right)^{-1}\neq 0,\quad j=1,2,$

• •

and $\theta_{j}=\arg\eta_{j}$ are ordered such that $\theta_{2}=\theta_{1}+\theta_{0}\ (\mod 2\pi)$ with a certain $\theta_{0}\in(0,\pi)$.

Then there exists a small enough $r>0$ and an analytic function $\omega(\zeta)$ defined for $\zeta\in\mathbb{D}_{r}^{2}$ such that $Q(\omega(\zeta);\zeta)=0$ everywhere in $\mathbb{D}_{r}^{2}$ and

[TABLE]

It can be seen from (A.1) that for triangles $T_{\delta}$ with small enough $\delta$ the images $\omega_{1}[T_{\delta}]=\{\omega_{1}(\zeta)\ :\zeta\in T_{\delta}\}$ and $\omega[T_{\delta}]=\{\omega(\zeta)\ :\zeta\in T_{\delta}\}$ lie inside a certain complex half-plane, for example, inside $\exp\left(\mathrm{i}(\theta_{1}-\delta_{0})\right)\mathbb{C}_{+}$ with $\delta_{0}$ small enough (we use here the assumptions that $0<\theta_{0}<\pi$ and $\theta_{1}<\theta_{2}<\theta_{1}+\pi$).

Now, considering the evolution as $\delta\to 0$ of images $\Sigma_{Q}[\partial T_{\delta}]$ of the triangle’s boundaries $\partial T_{\delta}$, it is easy to show the main statement of this section: for any small enough $\delta_{n}>0$, $n=1,2$, there exists $\delta_{3}>0$ such that

[TABLE]

The above arguments are the simplest particular case of a more general result [20] on multi-parameter perturbations of zeros of arbitrary finite multiplicity. The proof in [20] is more technically involved because Puiseux series have to be used for multiple zeros instead of the Taylor series (A.1).

Acknowledgments. ME and IK have been supported by the VolkswagenStiftung project “Modeling, Analysis, and Approximation Theory toward applications in tomography and inverse problems”. IK has been supported by the Alexander von Humboldt Foundation.