Inverse scattering with a non self-adjoint variational formulation

Daniel L. Marks; David R. Smith

doi:10.1364/OE.26.007655

1. Introduction

The design of optical components typically uses simplifying approximations of the interaction of fields and matter, for example, ray optics or weakly scattering media, so that the change in the medium required to produce a change in the fields is a tractable problem suitable for optimization. Exactly modeling the fields using the Finite Element Method (FEM) [1,2] or the Finite Difference Time Domain (FDTD) methods is computationally intensive and therefore reserved for situations for which the simpler models do not satisfactorily predict behavior. The changes to the medium needed are often difficult to ascertain using these more comprehensive models, and often it is difficult to optimize more than a few properties of a medium with these methods. As the field is ultimately a solution to the wave equation, both the field and the medium must satisfy this equation. The variational method specifies the solution to the wave equation as the stationary point of a functional that any variations of the fields, propagation constants, or the medium must maintain as stationary for the wave equation to continue to be satisfied. Therefore the variational method is ideal for examining what changes to the field and medium may jointly occur when searching for solutions to the inverse scattering problem. In this work, we examine the inverse scattering problem from the standpoint of the variational formulation of Maxwell’s equations, and use these to derive an iterative method for simultaneously inferring which fields and media satisfy a set of incident and scattered fields. This method is demonstrated by designing dielectric media that couple incident and scattered beams.

The FEM is a versatile and efficient method of calculating the fields scattered from complex structures. The FEM achieves this with simple steps: partitioning the medium into subdomains, defining trial functions over each subdomain, inserting these trial functions into the functional, and finding the stationary point of the functional. The functional is dependent on the medium permittivity and permeability and therefore the solution is specific to the medium, especially because steps like the partitioning and the choice of trial functions also often depends on the medium. What can be learned from the FEM about the solutions of a family of media, each with slightly perturbed permittivity and permeability? Typically the medium is not allowed to vary in the FEM, however, the FEM may be incorporated into iterative methods that may vary the medium, with the FEM used to repeatedly solve for the new fields due to the changed medium.

For an iterative inverse scattering method, it is often necessary to calculate the updates to a medium that are consistent with the measured field, with consistency defined by a minimized penalty function. This is often in the form of a gradient or functional derivative describing the change in the penalty function given a change in a medium parameter. Directly calculating the gradient is computationally costly as that involves solving the wave equation for every proposed change to the medium. The adjoint state method [3,4] is a method to reduce the computational burden of calculating gradients so that the medium parameters may be updated more efficiently. This efficiency is realized by solving for an adjoint field, which comprises the adjoint state, that describes the error contribution at each point in the medium proportional to the difference between the measured field and the field scattered by the current trial medium, which ideally should be the same if the medium is consistent with the measurements. This adjoint state may be used to calculate the gradient with respect to any parameters that vary the medium while the adjoint state needs to be solved for only once each iteration. Details on how to calculate these gradients for inverse scattering problems using the adjoint state method may be found in [5–7]. This method has been used to find efficient solutions to electromagnetic inverse scattering problems [8–17] by calculating the gradients of a penalty function efficiently with respect to medium parameters or obstacle shape parameters. The adjoint state gradient usually arises from defining a squared-error penalty between the measured field and the trial field, however, in the FEM, the adjoint solution is a fundamental tool used to solve non self-adjoint problems. In this work, the variational formulation of Maxwell’s equation is used to find consistent solutions for a medium, with the field and its adjoint treated equally. Rather than pose the inverse scattering problem as minimizing a penalty functional, the solution is found as a stationary point of a functional derived from Maxwell’s equations.

The variational principle in electromagnetics [18–24] is a means to transform the wave equation into an integral form with which weak solutions may be found. The functional may be interpreted as the Lagrangian of the electromagnetic field [25] and derived from the stationary action principle. This has been generalized to non self-adjoint problems [26,27] by introducing a functional that is stationary with respect to a pair of fields, one satisfying the medium and boundary conditions and the other the adjoint medium and boundary conditions. We extend this approach by identifying these two fields with the incident and scattered waves, each satisfying their own respective boundary conditions that reproduce the incident and scattered fields at the boundaries. The permittivity and permeability of the medium are allowed to vary as well as the fields so that the solution to the inverse scattering problem is a medium and pair of incident and scattered waves that is consistent with the boundary conditions.

In Section 2, the functional that is stationary when the solution to the wave equation is found is derived. Section 3 poses the inverse problem based on this functional. The problem of varying the medium to find an inverse scattering solution is addressed in Section 4. This is applied in Section 5 to design media that couple beams together. Several different methods of design are demonstrated, designing media that couple either one or two pairs of incident and scattered fields. Finally, conclusions are drawn in Section 6.

2. Deriving the functional for the non-self-adjoint Maxwell’s equations

The objective of the inverse scattering problem is to estimate the unknown properties of a medium from measurements of the waves that are scattered from the medium. The medium, confined to a volume V with a surrounding surface S, is characterized by a relative permittivity ∊_r and a permeability μ_r. The medium is inferred as that which is most consistent with scattering one or more incident fields E⁽⁺⁾ to respective scattered fields E⁽⁻⁾ on the surface S. This variational solution poses two problems to be simultaneously solved: a forward problem which models the propagation of the incident field forward in time through the medium and an adjoint problem that models the propagation of the scattered field backwards in time or backpropagated through the medium. Because the fields are modeled as time-harmonic, the adjoint solution is complex conjugated relative to the forward solution to model its backwards propagation. Ideally, a solution for the medium should be selected such that the forward and adjoint solutions are identical and therefore consistent with a single medium.

For the variational solution, a single functional F is stationary when its variation δF = 0 and both problems are simultaneously solved. Following [26], the solution to the forward problem can be posed as the solution to an operator equation L̂f = s, with L̂ being a linear operator that models the field including homogeneous boundary conditions, f being a field to be solved for, and s being a source function that includes any inhomogeneous boundary conditions. Likewise, the adjoint problem is L̂^af^a = s^a where L̂^a is the adjoint operator to L̂ including the adjoint homogeneous boundary conditions, f^a being an adjoint field to be solved for, and s^a being an adjoint source function, however, s^a is not necessarily determined by s. If we denote an inner product of two functions a and b over V by $〈 a, b 〉 = \int_{V} a^{*} \cdot b d V$ with the usual properties 〈a, b〉 = 〈b, a〉^* and 〈a, k₁b+ k₂c〉 = k₁ 〈a, b〉 + k₂ 〈a, c〉 for any functions a,b,c and constants k₁ and k₂, then the variation δF = 〈δf^a, L̂f − s〉 + 〈L̂^af^a − s^a, δf〉 = 0 for any variations δf and δf^a with the corresponding functional F = 〈f^a, L̂f〉 − 〈f^a, s〉 − 〈s^a, f〉. In the subsequent analysis, the variation and functional of Maxwell’s equations are explicitly derived using this inner product definition, using integration-by-parts in the form of Green’s theorem to transform the wave equation into a product of the incident and scattered fields.

In the volume V, the electromagnetic field obeys the time-harmonic Maxwell’s equations

\nabla \times E = i ω μ_{0} μ_{r} H

\nabla \times H = J - i ω ∊_{0} ∊_{r} E

with either of the following inhomogeneous boundary conditions on the closed surface S of volume V

\hat{n} \times (\hat{n} \times E) - η_{0} (\hat{n} \times H) = - E^{(+)}

\hat{n} \times (\hat{n} \times E) + η_{0} (\hat{n} \times H) = - E^{(-)}

with

η_{0} = \sqrt{μ_{0} / ∊_{0}}

. These boundary conditions approach the Sommerfeld radiation conditions as S surrounds the scatterer at an infinite distance. Eliminating the magnetic field from these relations:

\begin{array}{r} \hat{L} E = - i k^{- 1} η_{0} J for \hat{L} E = ∊_{r} E - k^{- 2} \nabla \times (μ_{r}^{- 1} \nabla \times E) \\ \hat{n} \times (\hat{n} \times E) + i k^{- 1} [\hat{n} \times (\nabla \times E)] = - E^{(+)} \end{array}

with

k = ω \sqrt{μ_{0} ∊_{0}}

and n̂ is a unit normal to the boundary. It is assumed that E⁽⁺⁾ · n̂ = 0 so that this boundary condition can hold. Because the boundary conditions are generally not self-adjoint, and in general μ_r and ∊_r are not real-valued so that L̂ is not self-adjoint, the solution is the stationary point of a variation and does not necessarily correspond to a maximum or minimum principle. This is because, in general, both the incident and scattered fields are being specified, and the scattered field may have less power than the incident field and the medium represented by μ_r and ∊_r would be a lossy medium to account for the difference in power. To find this stationary solution, we pose the adjoint problem using an adjoint operator and adjoint boundary conditions:

\begin{array}{r} {\hat{L}}_{a} E^{a} = - i k^{- 1} η_{0} J^{a} for {\hat{L}}_{a} E^{a} = ∊_{r}^{a} E^{a} - k^{- 2} \nabla \times (μ_{r}^{- 1 *} \nabla \times E^{a}) \\ \hat{n} \times (\hat{n} \times E^{a}) - i k^{- 1} [\hat{n} \times (\nabla \times E^{a})] = - E^{{(-)}^{*}} \end{array}

Again, it is assumed that E⁽⁻⁾ · n̂ = 0. Alternatively, this can be specified in terms of the complex conjugate of the adjoint field E^a*, so that the relationship to the non-adjoint solution is more obvious:

\begin{array}{r} \hat{L} E^{a *} = - i k^{- 1} η_{0} J^{a} for \hat{L} E^{a} = ∊_{r} E^{a *} - k^{- 2} \nabla \times (μ_{r}^{- 1} \nabla \times E^{a *}) \\ \hat{n} \times (\hat{n} \times E^{a *}) + i k^{- 1} [\hat{n} \times (\nabla \times E^{a *})] = - E^{(-)} \end{array}

We see that E^a = E^* when E⁽⁺⁾ = E⁽⁻⁾ and J^a = −J^*. The consistent solution to the inverse scattering problem is a medium μ_r and ∊_r such that with a given pair of incident and scattered fields, the condition E^a = E^* is satisfied. However, a medium may not exist so this condition is satisfied perfectly.

The solution is the fields E and E^a that simultaneously satisfy L̂E = −ik⁻¹η₀J and L̂^aE^a = −ik⁻¹η₀J^a. To do this, we find a functional F for which the solution to the problem is stationary. The variation of F is δF which must be zero for any variations δE and δE^a using the definition of the inner product as given above:

δ F = \int_{V} δ E^{a *} \cdot (\hat{L} E + i k^{- 1} η_{0} J) d V + \int_{V} δ E \cdot {({\hat{L}}^{a} E^{a} + i k^{- 1} η_{0} J^{a})}^{*} d V

Note that finding the functional that is made stationary using the variation of Eq. (8) does not preclude the possibility of varying ∊_r and μ_r. Substituting Eq. (5) and Eq. (6) into Eq. (8):

\begin{array}{l} δ F = \int_{V} ∊_{r} δ E^{a} \cdot E + i k^{- 1} η_{0} δ E^{a *} - J - k^{- 1} δ E^{a} \cdot \nabla \times [μ_{r}^{- 1} \nabla \times E] d V + \\ \int_{V} ∊_{r} δ E \cdot E^{a *} - i k^{- 1} η_{0} δ E \cdot J^{a *} - k^{- 1} δ E \cdot \nabla \times [μ_{r}^{- 1} \nabla \times E^{a *}] d V \end{array}

Apply the first vector Green’s identity:

\int_{V} a \cdot (\nabla \times b) d V = \int_{V} (\nabla \times a) \cdot b d V - \oint_{S} (a \times b) \cdot \hat{n} d S

the variation of Eq. (9) is

\begin{array}{l} δ F = & \int_{V} ∊_{r} δ E^{a *} \cdot E + i k^{- 1} η_{0} δ E^{a *} \cdot J - μ_{r}^{- 1} k^{- 2} {(\nabla \times δ E^{a})}^{*} \cdot (\nabla \times E) d V + \\ \int_{V} ∊_{r} δ E \cdot E^{a *} - i k^{- 1} η_{0} δ E \cdot J^{a *} - μ_{r}^{- 1} k^{- 2} (\nabla \times δ E) \cdot {(\nabla \times E^{a})}^{*} d V + \\ \oint_{S} k^{- 2} [δ E^{a *} \times (\nabla \times E)] \cdot \hat{n} d S + \oint_{S} k^{- 2} [δ E \times {(\nabla \times E^{a})}^{*}] \cdot \hat{n} d S \end{array}

Using the oriented volume of a parallelpiped identity a · (b× c) = b· (c × a) = c · (a × b) the surface integrals of Eq. (11) become

\begin{array}{l} δ F = & \int_{V} ∊_{r} δ E^{a *} \cdot E + i k^{- 1} η_{0} δ E^{a *} \cdot J - μ_{r}^{- 1} k^{- 2} {(\nabla \times δ E^{a})}^{*} \cdot (\nabla \times E) d V + \\ \int_{V} ∊_{r} δ E \cdot E^{a *} - i k^{- 1} η_{0} δ E \cdot J^{a *} - μ_{r}^{- 1} k^{- 2} (\nabla \times δ E) \cdot {(\nabla \times E^{a})}^{*} d V - \\ \oint_{S} k^{- 2} [\hat{n} \times (\nabla \times E)] \cdot δ E^{a *} d S - \oint_{S} k^{- 2} [\hat{n} \times {(\nabla \times E^{a})}^{*}] \cdot δ E d S \end{array}

The boundary conditions of Eq. (5) and Eq. (6) are now inserted:

\begin{array}{l} δ F = & \int_{V} ∊_{r} δ E^{a *} \cdot E + i k^{- 1} η_{0} δ E^{a *} \cdot J - μ_{r}^{- 1} k^{- 2} {(\nabla \times δ E^{a})}^{*} \cdot (\nabla \times E) d V + \\ \int_{V} ∊_{r} δ E \cdot E^{a *} - i k^{- 1} η_{0} δ E \cdot J^{a *} - μ_{r}^{- 1} k^{- 2} (\nabla \times δ E) \cdot {(\nabla \times E^{a})}^{*} d V + \\ \oint_{S} i k^{- 1} [\hat{n} \times (\hat{n} \times E) + E^{(+)}] \cdot \hat{n} \times {(\hat{n} \times δ E^{a})}^{*} d S + \\ \oint_{S} i k^{- 1} {[\hat{n} \times (\hat{n} \times E^{(a)}) + E^{{(-)}^{*}}]}^{*} \cdot \hat{n} \times (\hat{n} \times δ E) d S \end{array}

The functional F with the variation δF given by Eq. (13) is

\begin{array}{l} F = & \int_{V} ∊_{r} E^{a *} \cdot E + i k^{- 1} η_{0} (E^{a *} \cdot J - E \cdot J^{a *}) - μ_{r}^{- 1} k^{- 2} {(\nabla \times E^{a})}^{*} \cdot (\nabla \times E) d V + \\ \oint_{S} i k^{- 1} [\hat{n} \times (\hat{n} \times E) + E^{(+)}] \cdot {[\hat{n} \times (\hat{n} \times E^{a}) + E^{{(-)}^{*}}]}^{*} d S \end{array}

which can be verified by substituting E + δE for E and E^a + δE^a for E^a and expanding, and retaining only terms linear in δE and δE^a*.

3. Posing the inverse scattering problem

Ordinarily for the finite element method, to find the field E, E is expressed as a sum of basis functions with unknown coefficients, inserted into Eq. (14), and then F is differentiated with respect to the coefficients to find a linear system that determines the coefficients. However, in the inverse scattering problem, ∊_r, μ_r, E, and E^a are all unknown. Only the incoming and outgoing waves are known which in turn determine the boundary conditions using the Kirchoff integral theorem. In general, many sets of incoming and outgoing waves are needed to simultaneously infer the medium and the fields generated within the medium, with the medium having a common ∊_r and μ_r for all of the wave pairs. We assume that there are N sets of outgoing waves $E_{n}^{(+)}$ and incoming waves $E_{n}^{(-)}$ with a medium responding with the corresponding fields E_n and $E_{n}^{a}$ . In addition we add a weight w_n for each field and each field may have a different frequency k_n. A total functional F_T then has a stationary point for the solutions for ∊_r, μ_r, E_n, and $E_{n}^{a}$ :

\begin{array}{l} F_{T} = & \sum_{n = 1}^{N} w_{n} \int_{V} ∊_{r} E_{n}^{a *} \cdot E_{n} + i k_{n}^{- 1} η_{0} (E_{n}^{a *} \cdot J_{n} - E_{n} \cdot J_{n}^{a *}) - μ_{r}^{- 1} k_{n}^{- 1} {(\nabla \times E_{n}^{a})}^{*} \cdot (\nabla \times E_{n}) d V + \\ w_{n} \oint_{S} i k_{n}^{- 1} [\hat{n} \times (\hat{n} \times E_{n}) + E_{n}^{(+)}] \cdot {[\hat{n} \times (\hat{n} \times E_{n}^{a}) + E_{n}^{{(-)}^{*}}]}^{*} d S \end{array}

Unfortunately, because ∊_r and μ_r is unknown, one can not directly solve for the unknown coefficients as in Galerkin’s method to simultaneously find the stationary solution. However, one can perform a successive alternating estimation, where a trial solution of ∊_r and μ_r is assumed, and the Galerkin method is, for example, used to find E_n and

E_{n}^{a}

, with these solutions E_n and

E_{n}^{a}

in turn used as trial solutions to infer the new ∊_r and μ_r. Applying the finite element method is a computationally intensive process because it involves generating and solving a linear system of equations for the unknown coefficients of the basis functions. For a successive alternating estimation process to find a solution, it is only required that the estimate of the field be slightly improved each iteration; the residuals do not need to be completely minimized. Instead, we suggest a fixed point iteration method which performs successive updates to the field and medium which is more naturally evocative of the process of the holographic exposure of a medium. We can examine if an error reduction fixed-point iteration may be used to find improved estimates of E and E^a. The equations to be solved for the fields are L̂_nÊ_n = −ik⁻¹η₀J_n and

{\hat{L}}_{n}^{a} E_{n}^{a} = - i k^{- 1} η_{0} J_{n}^{a}

. If all of the eigenvalues of the linear operators L̂_n and

{\hat{L}}_{n}^{a}

had a positive real part, the following iterations would find the solutions by subtracting the error in the equation given a small enough step size β:

{E_{n}}^{(p + 1)} = {E_{n}}^{(p)} - β ({\hat{L}}_{n} {E_{n}}^{(p)} + i k^{- 1} η_{0} J_{n})

E_{n}^{a (p + 1)} = E_{n}^{a (p)} - β ({\hat{L}}_{n}^{a} E_{n}^{a (p + 1)} + i k^{- 1} η_{0} J_{n}^{a})

However, in general the operators L̂_n and

{\hat{L}}_{n}^{a}

have eigenvalues with negative real parts that prevent these simple fixed point iterations from converging. To circumvent this problem, these iterations may be modified to be:

{E_{n}}^{(p + 1)} = {E_{n}}^{(p)} - β {\hat{B}}_{n}^{a} {\hat{L}}_{n}^{a} ({\hat{L}}_{n} {\hat{B}}_{n} {E_{n}}^{(p)} + i k^{- 1} η_{0} J_{n})

E_{n}^{a (p + 1)} = E_{n}^{a (p)} - β {\hat{B}}_{n} {\hat{L}}_{n} ({\hat{L}}_{n}^{a} {\hat{B}}_{n}^{a} E_{n}^{a (p + 1)} + i k^{- 1} η_{0} J_{n}^{a})

The compositions

{\hat{L}}_{n}^{a} {\hat{L}}_{n}

and

{\hat{L}}_{n} {\hat{L}}_{n}^{a}

are of an operator and its adjoint, so the eigenvalues of the combination are nonnegative and therefore the iteration is stable and converges with a sufficiently small β. The projection operators B̂_n and

{\hat{B}}_{n}^{a}

may be included to project the fields back to feasible boundary conditions. For example, if part of the boundary is a perfect electric conductor, the operator would set the tangential electric field to be zero at that part of the boundary. Because a projection operator has eigenvalues of zero or one, it does not increase the magnitude of any eigenvalues and therefore does not prevent convergence of the iterations. While the above iteration is an example that may be used to perform small updates to the fields each iteration that are consistent with the wave equation, in practice it is much more expedient to apply a small number of iterations of the conjugate gradient method to the equations

{\hat{B}}_{n}^{a} {\hat{L}}_{n}^{a} {\hat{L}}_{n} {\hat{B}}_{n} {E_{n}}^{(p + 1)} = - i k^{- 1} η_{0} {\hat{B}}_{n}^{a} {\hat{L}}_{n}^{a} J_{n}

and

{\hat{B}}_{n} {\hat{L}}_{n} {\hat{L}}_{n}^{a} {\hat{B}}_{n}^{a} E_{n}^{a (p + 1)} = - i k^{- 1} η_{0} {\hat{B}}_{n} {\hat{L}}_{n} J_{n}^{a}

, but it is not necesssary to solve the equations until a small residual is achieved.

Rather than directly enforcing the boundary condition, it is easier to use an equivalent current source for the inhomogeneous boundary terms with an absorbing medium or a perfectly matched layer at the boundary. The use of the absorbing layer is advantangeous if the entire incoming or outgoing field on the surface S is not known so that the unknown part of the boundary does not reflect any waves. If we examine Eq. (13), given that the equivalent current source is tangential to the boundary, it may be rewritten as:

\begin{array}{l} δ F = & \int_{V} ∊_{r} δ E^{a *} \cdot E - μ_{r}^{- 1} k^{- 1} {(\nabla \times δ E^{a})}^{*} \cdot (\nabla \times E) d V + \\ \int_{V} ∊_{r} δ E^{a *} - μ_{r}^{- 1} k^{- 2} (\nabla \times δ E) \cdot {(\nabla \times E^{a})}^{*} d V + \\ \oint_{S} i k^{- 1} [\hat{n} \times (\hat{n} \times E) + E^{(+)} + η_{0} J Δ] \cdot \hat{n} \times {(\hat{n} \times δ E^{a})}^{*} d S + \\ \oint_{S} i k^{- 1} {[\hat{n} \times (\hat{n} \times E^{a}) + E^{{(-)}^{*}} + η_{0} J^{a} Δ]}^{*} \cdot \hat{n} \times (\hat{n} \times δ E) d S \end{array}

The dimension Δ is thickness of a thin layer of current placed at the surface S, and is in practice the spatial sampling rate of the volume. Equivalent current sources

J Δ = η_{0}^{- 1} E^{(+)}

and

J^{a} Δ = η_{0}^{- 1} E^{{(-)}^{*}}

may be substituted for the incident and scattered electric fields at the boundary.

4. Varying the medium to find a solution to the inverse scattering problem

When varying ∊_r and μ_r to find the medium, we only allow variations of ∊_r and μ_r not on the boundary. Without the boundary terms, F_T is

F_{T} = \int_{V} ∊_{r} [\sum_{n = 1}^{N} w_{n} E_{n}^{a *} \cdot E_{n}] - μ_{r}^{- 1} [\sum_{n = 1}^{N} w_{n} k_{n}^{- 1} {(\nabla \times E_{n}^{a})}^{*} \cdot (\nabla \times E_{n})] - K d V

The function K is a constant function of position that has been included but does not influence the stationary point of the functional. To determine the functional derivatives of F_T with respect to μ_r and ∊_r, we examine the variation of F_T :

δ F_{T} = \int_{V} δ ∊_{r} [\sum_{n = 1}^{N} w_{n} E_{n}^{a *} \cdot E_{n}] + δ μ_{r} μ_{r}^{- 2} [\sum_{n = 1}^{N} w_{n} k_{n}^{- 2} {(\nabla \times E_{n}^{a})}^{*} \cdot (\nabla \times E_{n})] d V

One might consider performing gradient descent using these functional derivatives to minimize F_T as in the following iterations:

μ_{r}^{(p + 1)} = μ_{r}^{(p)} - γ_{μ} \frac{\partial F_{T}}{\partial μ_{r}}

∊_{r}^{(p + 1)} = ∊_{r}^{(p)} - γ_{∊} \frac{\partial F_{T}}{\partial ∊_{r}}

\frac{\partial F_{T}}{\partial μ_{r}} = μ_{r}^{- 2} \sum_{n = 1}^{N} w_{n} k_{n}^{- 2} {(\nabla \times E_{n}^{a})}^{*} \cdot (\nabla \times E_{n})

\frac{\partial F_{T}}{\partial ∊_{r}} = \sum_{n = 1}^{N} w_{n} E_{n}^{a *} \cdot E_{n}

where γ_μ and γ_∊ are step sizes. Unfortunately, in general this does not work, as the solution to the functional F_T is stationary, but it is not necessarily a minimum. If instead the problem is posed as a self-adjoint problem for which L̂^a = L̂ and therefore ∊_r and μ_r are real functions:

\begin{array}{r} \hat{L} E = - i k^{- 1} η_{0} J and \hat{L} E^{a} = - i k^{- 1} η_{0} J^{a} \\ \hat{n} \times (\hat{n} \times E) + i k^{- 1} [\hat{n} \times (\nabla \times E)] = - E^{(+)} \\ \hat{n} \times (\hat{n} \times E^{a}) - i k^{- 1} [\hat{n} \times (\nabla \times E^{a})] = - E^{(-)} \end{array}

The stationary solution for a medium consistent with both fields E and E^a is E^a = E rather than E^a = E^* as in the non self-adjoint case. Because in the self-adjoint case the stationary point of the functional F_T is a true minimum, using the above the gradient descent steps finds a minimum and yields a solution with the substitution

E_{n}^{a} \to E_{n}^{a *}

given the problem as originally posed and real-valued ∊_r and μ_r. In terms of the non self-adjoint problem defined by Eqs. (5) and (6), but with real-valued ∊_r and μ_r, the iterations of Eq. (24) are:

μ_{r}^{(p + 1)} = μ_{r}^{(p)} - γ_{μ} μ_{r}^{- 2} \sum_{n = 1}^{N} w_{n} k_{n}^{- 2} (\nabla \times E_{n}^{a}) \cdot (\nabla \times E_{n})

∊_{r}^{(p + 1)} = ∊_{r}^{(p)} - γ_{∊} \sum_{n = 1}^{N} w_{n} E_{n}^{a} \cdot E_{n}

Using the self-adjoint solution, an interpretation of Eq. (22) is that if $E_{n}^{a} = E_{n}$ , the term δ∊_rE_n^* · E_n is the change in the electric energy density stored in the medium with the variation δ∊_r and $δ μ_{r} μ_{r}^{- 2} {(\nabla \times E_{n})}^{*} \cdot (\nabla \times E_{n})$ is the change in the magnetic energy density stored in the medium with the variation δµ_r. The quantities $∊_{r} E_{n}^{a *} \cdot E_{n}$ and $μ_{r}^{- 1} {(\nabla \times E_{n}^{a})}^{*} \cdot (\nabla \times E_{n})$ may then be regarded as the electric and magnetic “mutual energy densities” which is the energy stored in the medium as it is being transferred between fields. The updates of Eq. (23) and Eq. (24) find the solution such that $E_{n}^{a} = E_{n}$ while minimizing the energy stored by the medium when transferring energy between the fields which is the stationary action principle applied to time-harmonic functions.

The solution obtained by applying the updates of Eqs. (28) and (29) tends to be a “diffractive” solution that obtains grating-like structures that diffract the incident to the scattered field. However, there is an alternative that produces more refractive-like structures which is called the “refractive” solution. Another way to make the variation δF_T zero is to make F_T a constant given E_n and $E_{n}^{a}$ by setting its integrand to zero. Solving Eq. (21) for the relationship between ∊_r and μ_r:

μ_{r} ∊_{r} = \frac{K μ_{r} + \sum_{n = 1}^{N} w_{n} k_{n}^{- 1} {(\nabla \times E_{n}^{a})}^{*} \cdot (\nabla \times E_{n})}{\sum_{n = 1}^{N} w_{n} E_{n}^{a *} \cdot E_{n}}

As the local wave number in a medium is

k \sqrt{μ_{r} ∊_{r}}

, this may be interpreted as finding the local wave number that couples the incident and scattered fields together for the case K = 0. For example, if

E_{n}^{a *} = E_{n} = exp (i k_{n} n u \cdot r)

, where |u| = 1, one can verify that μ_r∊_r = n², or the refractive index squared. There is a choice of whether to change μ_r, ∊_r, or both, so that this solution is incomplete. Furthermore, when the magnitude of the denominator of Eq. (30) is small, the solution is ill-conditioned. To mitigate this, we propose modifying Eq. (30) so that in the limit of small magnitude fields, the relation tends to a constant function of μ_r∊_r:

μ_{r} ∊_{r} = \frac{K μ_{r} + \sum_{n = 1}^{N} w_{n} E_{m}^{4} + \sum_{n = 1}^{N} w_{n} k_{n}^{- 2} {(\nabla \times E_{n}^{a})}^{*} \cdot (\nabla \times E_{n}) {(E_{n}^{a} \cdot E_{n})}^{*}}{\sum_{n = 1}^{N} w_{n} E_{m}^{4} + \sum_{n = 1}^{N} w_{n} {| E_{n}^{a *} \cdot E_{n} |}^{2}}

For |E| ≪ E_m and |E^a| ≪ E_m, Kμ_r may be selected to be a constant function of position which serves as a prior to determine μ_r∊_r in the absence of fields at a given point. Eq. (31) is consistent with Eq. (30) for a single field E_n, but changes how all N fields are combined to determine the medium properties. As the purpose of the constant E_m is to provide a stable solution when fields are not present, the magntiude of E_m should be selected to be a small fraction (e.g. 0.01 to 0.1) of that of the field strength expected to be present after the wave equation is solved. Too small E_m fails to provide a stable solution so that undetermined points are overly sensitive to numerical error, while too large E_m does not allow the medium to vary to satisfy the inverse scattering problem.

A crucial difference between the diffractive solution and the refractive solution is that while the diffractive solution may produce perturbations to the medium that change ∊_r and μ_r both positively and negatively, the refractive solution always decreases the wave number. While this may seem a disadvantage for the refractive solution, as it might afford fewer possibilities for optimization, it is an advantage when designing continuous media. A grating-like structure, when smoothed, tends to cancel positive and negative perturbations in the medium which reduces the scattering strength of the grating. Because the refractive solution only decreases the wave number, smoothing the refractive solution does not remove the medium perturbations. A medium with smooth features may be easier to fabricate than one with periodic fine features. However, if a limited modulation range is available in a medium, a diffractive solution is more likely to produce a solution with a higher efficiency, analogous to how a volume hologram, even with weak dielectric perturbations, can achieve a high diffraction efficiency that is unlikely to be achieved by a refractive medium with the same modulation range.

In practice, as we know that at the solution E^a = E^*, it is helpful to form affine combinations of the two so that iterations containing Eq. (31) are more stable. Instead, if $0 \leq α \leq \frac{1}{2}$ is a parameter, two fields $E_{n}^{A} = E_{n}^{a *} α + E_{n} (1 - α)$ and $E_{n}^{B} = E_{n}^{a} (1 - α) + E_{n}^{*} α$ may be inserted into Eq. (31) without changing the fixed point E^a = E^*:

μ_{r} ∊_{r} = \frac{K μ_{r} \sum_{n = 1}^{N} w_{n} E_{m}^{4} + \sum_{n = 1}^{N} w_{n} k_{n}^{- 2} {(\nabla \times E_{n}^{B})}^{*} \cdot (\nabla \times E_{n}^{A}) {(E_{n}^{B *} \cdot E_{n}^{A})}^{*}}{\sum_{n = 1}^{N} w_{n} E_{m}^{4} + \sum_{n = 1}^{N} w_{n} {| E_{n}^{B *} \cdot E_{n}^{A} |}^{2}}

Likewise, the updates in Eqs. (28) and (29) become:

μ_{r}^{(p + 1)} = μ_{r}^{(p)} - γ_{μ} μ_{r}^{- 2} \sum_{n = 1}^{N} w_{n} k_{n}^{- 1} (\nabla \times E_{n}^{B}) \cdot (\nabla \times E_{n}^{A})

∊_{r}^{(p + 1)} = ∊_{r}^{(p)} - γ_{∊} \sum_{n = 1}^{N} w_{n} E_{n}^{B} \cdot E_{n}^{A}

To determine how well a medium couples the incident and scattered fields together, we provide a figure of merit. Lorentz reciprocity states that the reaction of the incident field on the current distribution that generates the scattered field must equal the reaction of the scattered field on the current distribution that generates the incident field:

\int_{V} J \cdot E^{a *} d V = \int_{V} J^{a *} \cdot E d V

The volume V contains all sources so there is no surface term. The reaction is the total power coupled between the incident and scattered fields through the medium. Therefore we propose the following figure of merit:

P = | \sum_{n = 1}^{N} \int_{V} J_{n} \cdot E_{n}^{a *} + J_{n}^{a *} \cdot E_{n} d V |

Two media may be compared and the desired solution is that which couples the most power between the fields. Furthermore, the two reactions may be compared to determine if the solution is consistent with Lorentz reciprocity and therefore should be accepted.

5. A demonstration of medium design

A demonstration of finding a dielectric medium consistent with incident and scattered fields in a two-dimensional medium is described based on finding stationary solutions to the functional F. As a simplified problem, a non-magnetic three-dimensional medium and field is simulated that is uniform in the z-direction which can be modeled as a scalar problem in the transverse magnetic polarization with just the ẑ-component of the electric field, so that E = Eẑ. The medium volume is sampled at a spatial sampling rate of Δ. We describe a process that finds a medium and the fields within the medium consistent with the provided boundary conditions. To find the electric fields consistent with the current medium and boundary conditions, the updates of the fields E and E^a as given by Eqs. (18) and (19) are:

\begin{array}{l} {\hat{L}}_{n} E = ∊_{r} E + k_{n}^{- 1} [\frac{\partial^{2} E}{\partial x^{2}} + \frac{\partial^{2} E}{\partial y^{2}}] \\ {\hat{L}}_{n}^{a} E^{a} = ∊_{r}^{*} E^{a} + k_{n}^{- 2} [\frac{\partial^{2} E^{a}}{\partial x^{2}} + \frac{\partial^{2} E^{a}}{\partial y^{2}}] \\ {E_{n}}^{(p + 1)} = {E_{n}}^{(p)} - β {\hat{L}}_{n}^{a} ({\hat{L}}_{n} {E_{n}}^{(p)} + i {(k Δ)}^{- 1} E_{n}^{(+)}) \\ E_{n}^{a (p + 1)} = E_{n}^{a (p)} - β {\hat{L}}_{n} ({\hat{L}}_{n}^{a} E_{n}^{a (p)} + i {(k Δ)}^{- 1} E_{n}^{{(-)}^{*}}) \end{array}

These updates are accelerated by means of application of iterations of the conjugate gradient algorithm to the linear systems

{\hat{L}}_{n}^{a} {\hat{L}}_{n} {E_{n}}^{(p + 1)} = - i {(k Δ)}^{- 1} {\hat{L}}_{n}^{a} {E_{n}}^{(+)}

and

{\hat{L}}_{n} {\hat{L}}_{n}^{a} E_{n}^{a (p + 1)} = - i {(k Δ)}^{- 1} {\hat{L}}_{n} E_{n}^{{(-)}^{*}}

. The two-dimensional Laplacian is implemented by convolution with the 9-point stencil (kΔ)⁻² {

\frac{1}{6}

,

\frac{2}{3}

,

\frac{1}{6}

;

\frac{2}{3}

,

- \frac{10}{3}

,

\frac{2}{3}

;

\frac{1}{6}

,

\frac{2}{3}

,

\frac{1}{6}

}. The updates of Eq. (34) and Eq. (32) become

∊_{r}^{(p + 1)} = ∊_{r}^{(p)} - \frac{γ_{∊}}{E_{m}^{2}} \sum_{n = 1}^{N} w_{n} E_{n}^{B (p)} \cdot E_{n}^{A (p)}

\begin{array}{l} ∊_{r}^{(p + 1)} = (1 - γ_{∊}) ∊_{r}^{(p)} + \\ γ_{∊} = \frac{∊_{r}^{(0)} \sum_{n = 1}^{N} w_{n} E_{m}^{4} + \sum_{n = 1}^{N} w_{n} k_{n}^{- 2} [\frac{\partial E_{n}^{B {(p)}^{*}}}{\partial x} \frac{\partial E_{n}^{A (p)}}{\partial x} + \frac{\partial E_{n}^{B {(p)}^{*}}}{\partial y} \frac{\partial E_{n}^{A (p)}}{\partial y}] (E_{n}^{B {(p)}^{*}} E_{n}^{A {(p)}^{*}})}{\sum_{n = 1}^{N} w_{n} E_{m}^{4} + \sum_{n = 1}^{N} w_{n} {| E_{n}^{B^{*} (p)} E_{n}^{A (p)} |}^{2}} \end{array}

with

E_{n}^{A} = E_{n}^{a *} α + E_{n} (1 - α)

and

E_{b}^{B} = E_{n}^{a} (1 - α) + E_{n}^{*} α

. Note that the squared quantities in the numerator are not absolute values, so that the numerator unlike the denominator is not necessary real-valued. The constant γ_∊ is a step size chosen small enough so that the iterations approach a desired solution. E_m serves as a regularization constant determining the field strength needed to change ∊_r. Eq. (38) is the diffractive solution for ∊_r while Eq. (39) is the refractive solution. The function

∊_{r}^{(0)}

is the initial state of the permittivity before iterations start.

The foregoing iteration formulas do not describe an optimization process. As noted earlier, the solutions that simultaneously find the fields and the medium consistent with the wave equation are the stationary points of the functional F, but these stationary points do not necessarily correspond to minima of this functional. In the special case of a lossless medium, the operator L̂ is self-adjoint and therefore the functional has minima that may be found, for example, using the update of Eq. (38). However, in general this is not true, and a stationary solution corresponds to a saddle point, and furthermore this saddle point is usually on a constraint boundary which further complicates the stability analysis by introducing unknown Lagrange multipliers that influence the curvature of the functional at the saddle point. The iterative process finds these saddle points but this is an inherently unstable process without other constraints on the solution. Without additional constraints, the iterations indeed tend to approach a solution such that E^a = E^* but as this condition can never be perfectly attained, further iteration eventually overshoots and repels from the saddle point, and eventually a new saddle point is found and the process repeated. Smaller step sizes only mitigate this problem. Another criterion must be specified to distinguish between which of the multitude of saddle points is the desired solution. We use the sum of the reactions in Eq. (36) as a reward function to decide if the medium is efficient at transferring the energy from the incident to the scattered wave and vice versa.

At the perimeter of the volume, outside of the surface where the surface currents are applied to produce the incident and scattered wave boundary conditions, an absorbing medium is placed to prevent reflections at the boundary. Because the simulation is of an isotropic, nonmagnetic medium, a perfectly matched layer can not be used, however, a simple absorbing layer produced sufficient attenuation of the outgoing waves. In the simulation, a band of 20 samples is placed around the edge, with ∊_r varying linearly from 1 + 0.0i at the boundary where the surface currents are applied to 1 + 0.5i at the edge of the sampled area. At the actual edge of the sampled field, the Laplacian operator used Dirichlet boundary conditions. The absorbing layer was not modified during the iterative process.

As a first example, an incident beam at 45° from the lower left corner of a 14λ × 14λ volume is intended to be coupled as a scattered wave at 60° at the lower right corner, so the beam is diverted 105° in the space of 14λ. Between the left and right edges is a slab of dielectric of nominal ∊_r = 3 which is designed to coupled the incident to the scattered beam. Three solutions are demonstrated to this problem:

The diffractive solution, which tends to create gratings that couple the incident and scattered waves.
The refractive solution without smoothing, which also tends to create gratings, however, these tend to be somewhat more uniform in refractive index.
The refractive solution with smoothing, which creates a graded refractive-index medium.

For all of the simulations, E_m is set to 0.25 times the peak absolute value of either of the incident or scattered fields, α = 0.45, and γ_∊ = 0.05. An initial 3000 iterations of the conjugate gradient method are used to find the field in the dielectric slab, and then afterwards 200 iterations of the conjugate gradient method alternate with the updates of either Eq. (38) or (39).

The diffractive solution is shown in Fig. 1. The upper left panel is the real part of the incident/non-adjoint field and the upper right panel is the real part of the scattered/adjoint field. The lower left panel shows the permittivity distribution created by the algorithm. Finally, the lower right figure shows the magnitude of both the incident and scattered fields, with the incident in the red channel, and the scattered in the green channel. Note that the figure is mostly yellow, showing that the two have largely overlapped and that a solution has been found that couples the incident to the scattered wave. The holographic solution is a diffraction grating, as expected, however, because there are standing waves inside the dielectric slab, these waves affect the coupling. The grating incorporates the standing waves in the dielectric slab as the iterations progress. Note that the permittivty in the slab is both above and below the nominal value ∊_r = 3.

Fig. 1 The diffractive solution to the coupling an incident to a scattered beam. The beam is incident from the lower left at an angle of 45° and scattered to the lower right at an angle of 60°. The upper left panel is the real part of the incident/non-adjoint electric field, the upper right panel is the real part of the scattered/adjoint field, the lower left is the final permittivity of the medium. the lower right panel shows both field magnitudes, the red channel being the non-adjoint field and the green channel being the adjoint field. The animation of these iterations is shown in Visualization 1.

Download Full Size | PDF

The non-smoothed refractive solution is shown in Fig. 2. Like the holographic solution, the non-smoothed solution is a grating, but the permittivity in the slab does not exceed the nominal value ∊_r = 3. This solution also overlaps the incident and scattered fields well, as evidenced by the largely yellow color in the lower right panel. Because the permittivity modulation is in the negative direction, there is a refractive component to the scattering, as unlike in the holographic case where the incident wave sharply bends into the diffracted beam, the refractive solution more gently bends the wave into an arc through the slab.

Fig. 2 The refractive solution, without smoothing, to the coupling an incident to a scattered beam. The beam is incident from the lower left at an angle of 45° and scattered to the lower right at an angle of 60°. The upper left panel is the real part of the incident/non-adjoint electric field, the upper right panel is the real part of the scattered/adjoint field, the lower left is the final permittivity of the medium. the lower right panel shows both field magnitudes, the red channel being the non-adjoint field and the green channel being the adjoint field. The animation of these iterations is shown in Visualization 2.

Download Full Size | PDF

The smoothed refractive index solution is shown in Fig. 3. To smooth the refractive index, a Lagrange multiplier was applied to the permittivity to maintain the total squared gradient at a particular value. As can be seen, the medium is indeed much smoother, and the incident beam travels in a smooth arc across the slab. However, the overlap between the fields is not as good as in the other two cases, and so the efficiency is reduced. The blue region in the permittivity graph corresponds to a low permittivity, and the wave may be seen as conforming to the edge of this region, as the wave is pinned to the boundary between high and low refractive index, similar to propagation in a graded refractive index waveguide.

Fig. 3 The refractive solution, with smoothing, to the coupling an incident to a scattered beam. The beam is incident from the lower left at an angle of 45° and scattered to the lower right at an angle of 60°. The upper left panel is the real part of the incident/non-adjoint electric field, the upper right panel is the real part of the scattered/adjoint field, the lower left is the final permittivity of the medium. the lower right panel shows both field magnitudes, the red channel being the non-adjoint field and the green channel being the adjoint field. The animation of these iterations is shown in Visualization 3.

Download Full Size | PDF

The mutual reactions of the incident and scattered fields are plotted in Fig. 4 for all three simulations. Note that the blue and red curves, which represent the reaction of the incident field source on the scattered field, and the scattered field source on the incident field, respectively, are both nearly equal, showing that Lorentz reciprocity is maintained during the solution. Furthermore, the black curve, which is the reward function P and is the sum of these two reactions, is optimized by this process, but the iterations are not designed specifically to optimize P.

Fig. 4 The coupled power between the incident and the scattered fields at each iteration, for (a) the diffractive solution, (b) the refractive solution without smoothing, and (c) the refractive solution with smoothing. While the power is specified in relative units, these may be compared between the three solutions because the same fields are being coupled in all three simulations. The blue curves are the reaction of the incident field source on the scattered field, the red curves are the reaction of the scattered field source on the incident field, and the black curve is the sum of the two, or the reward function P. The blue and red curves follow closely, showing that Lorentz reciprocity is broadly satisfied. Furthermore, the method does optimize P, despite the fact it is not specifically designed to do so.

Download Full Size | PDF

As a second demonstration, a medium was designed to couple two incident waves to their two respective scattered waves. The first pair is a wave incident from the lower left at 60° inclination to be scattered to a wave that is 0° inclination at the upper right. The second pair is a wave that is −20° inclined from the upper left, scattered to a wave that is −40° at the lower right. The waves pass through each other, so the medium through which the waves intersect must simultaneously accommodate both beams. A diffractive structure may easily operate on both beams as the gratings that diffract each beam largely do not interact due to the Bragg selection effect. However, a refractive structure is more difficult as the refractive structures affect both beams and must refract both incident waves to their respective scattered waves.

We first examine the diffractive solution of Fig. 5. Both incident waves can be efficiently coupled to their respective output waves, as expected. The visualizations shows the gradual formation and “capture” of the wavefronts so that the incident and scattered waves are coupled together. The refractive solution of Fig. 6 is more interesting because the beams must cross. The permittivity is shaped so that the fields are formed into filament-like waveguide structures with the crossing point between the filaments being a region of more uniform permittivity near the center of the slab. This solution was obtained completely automatically, initially seeding the slab with a uniform refractive index ∊_r = 3. The animation of these iterations shows the two fields “rippling” back and forth as they both try to find a location where the fields might pass each other unimpeded. The coupled power between the two field pairs is plotted at each iteration in Fig. 7. The blue curve is the total coupled power P reward function for the first field pair with the incident field at 60° inclination and the scattered field at 0°, and the red curve is P for the second field pair with the incident field at −20° and the scattered field at −40°. The diffractive solution monotonically approaches an efficient solution as volume gratings are an efficient way to divert a beam. The refractive solution shows the competition between the two wavefronts as these are at first attempting to alter the same regions of refractive index oppositely, so that the coupled power for the two wave pairs oscillates between them. Finally, these settle into a relationship where each does not obtain efficiency at the expense of the other, and both steadily improve.

Fig. 5 The diffractive solution to the coupling two incident beams to two scattered beams. The left column is the real part of the incident/non-adjoint electric fields, the second column is the real part of the scattered/adjoint fields. The third column is the final permittivity of the medium. The right column shows both field magnitudes, the red channel being the non-adjoint field and the green channel being the adjoint field. The animation of these iterations is shown in Visualization 4.

Download Full Size | PDF

Fig. 6 The refractive solution to the coupling two incident beams to two scattered beams. The left column is the real part of the incident/non-adjoint electric fields, the second column is the real part of the scattered/adjoint fields. The third column is the final permittivity of the medium. The right column shows both field magnitudes, the red channel being the non-adjoint field and the green channel being the adjoint field. The animation of these iterations is shown in Visualization 5.

Download Full Size | PDF

Fig. 7 The coupled power between the two pairs of incident and the scattered fields at each iteration, for (a) the diffractive solution and (b) the refractive solution with smoothing. While the power is specified in relative units, these may be compared between the solutions because the same fields are being coupled in both simulations. The blue curve is the total coupled power P for the first set of incident and scattered waves, and the red curve is the total coupled power for the second set of waves. While the diffractive solution monotonically approaches the solution, the efficiency oscillates between the two wave pairs for the refractive solution as a solution is found which is consistent with both being coupled efficiently.

Download Full Size | PDF

The computational performance is limited by the need to simultaneously satisfy the boundary conditions and the wave equation in the medium. The updates of Eq. (37) ensure consistency with the wave equation, however, the equivalent current is only present at the boundary, and therefore its effects must propagate throughout the medium to ensure consistency with any changes in the medium so that many iterations are required to converge on a solution. Inverse scattering methods based on Green’s functions, such as the Born Iterative Method [28] and the Modified Gradient Method [29], have the advantage that the Green’s function may be used to rapidly compute the effects of changing the medium on the field. However, if the perturbation in the medium is too great, these methods may not converge [30]. To improve the performance of this method, ideally a method could be devised that ensures consistency in a nonlocal way in a manner similar to the Green’s function solutions, unlike the strictly local updates of Eq. (37) between adjacent field samples, but would not require that the variation of the medium be a perturbation of a background medium. The Green’s function solution is available because of the perturbation assumption of the medium, and in general this is insufficient for many design problems, especially for designing graded refractive index structures which are poorly modeled by perturbative media.

6. Conclusion

By examining the non-self-adjoint variational formulation Maxwell’s equations, we have created and demonstrated an iterative method for finding an electromagnetic medium that couples one or more incident waves to their respective scattered waves. The functional itself served as a means of determining if the solution is consistent and it was not required to use a cost penalty to ensure a fit to measurements. The disadvantage of this approach is that in general the solution is a stationary point of the functional and not a true minimum with a lower bound so that the solution may be more difficult to find. While a simple means of solving the wave equation based on the conjugate gradient method on the normal equations was used here, more sophisticated methods are likely more efficient, as long as these methods can be formulated so that the problem may be solved iteratively and incrementally rather than completely at each iteration. The boundary element method, for example, may be solved incrementally and is a useful way to incorporate the Green’s function and constrain the geometry of the solution. As the efficiency of a calculation is often determined how specific it is to the problem domain, and the described method is perhaps overly general, it might be customized by incorporating better field solvers to attack larger problems. Nevertheless, the demonstrated method shows that the inverse scattering problem can be solved using very few assumptions and therefore is highly adaptable.

Funding

Air Force Office of Scientific Research (FA9550-12-1-0491).

References and links

1. J.-M. Jin, The Finite Element Method in Electromagnetics (Wiley, 2014).

2. J. L. Volakis, A. Chatterjee, and L. C. Kempel, Finite Element Method for Electromagnetics (Wiley, 1998). [CrossRef]

3. R.-E. Plessix, “A review of the adjoint-state method for computing the gradient of a functional with geophysical applications,” Geophys. J. Int. 167, 495–503 (2006). [CrossRef]

4. M. B. Giles and N. A. Pierce, “An introduction to the adjoint approach to design,” Flow. Turbul. Combust. 65, 393–415 (2000). [CrossRef]

5. P. R. McGillivray and D. W. Oldenburg, “Methods for calculating Frechet derivatives and sensitivities for the non-linear inverse problem: a comparative study,” Geophys. Prospect. 38, 499–524 (1990). [CrossRef]

6. S. J. Norton, “Iterative inverse scattering algorithms: Methods of computing Frechet derivatives,” J. Acoust. Soc. Am. 106, 2653–2660 (1999). [CrossRef]

7. F. Hettlich, “Frechet derivatives in inverse obstacle scattering,” Inverse Probl. 11, 371 (1995). [CrossRef]

8. I. T. Rekanos and T. D. Tsiboukis, “A combined finite element-nonlinear conjugate gradient spatial method for the reconstruction of unknown scatterer profiles,” IEEE T. Magn. 34, 2829–2832 (1998). [CrossRef]

9. I. T. Rekanos, T. V. Yioultsis, and T. D. Tsiboukis, “Inverse scattering using the finite-element method and a nonlinear optimization technique,” IEEE T. Microw. Theory 47, 336–344 (1999). [CrossRef]

10. I. T. Rekanos, “Inverse scattering in the time domain: an iterative method using an FDTD sensitivity analysis scheme,” IEEE T. Magn. 38, 1117–1120 (2002). [CrossRef]

11. D. Bertsimas, O. Nohadani, and K. M. Teo, “Robust optimization in electromagnetic scattering problems,” J. Appl. Phys. 101, 074507 (2007). [CrossRef]

12. Y. Cao, J. Xie, Y. Liu, and Z. Liu, “Modeling and optimization of photonic crystal devices based on transformation optics method,” Opt. Express 22, 2725–2734 (2014). [CrossRef] [PubMed]

13. H. Khalil, S. Bila, M. Aubourg, D. Baillargeat, S. Verdeyme, F. Jouve, C. Delage, and T. Chartier, “Shape optimized design of microwave dielectric resonators by level-set and topology gradient methods,” Int. J. RF Microw. C. E. 20, 33–41 (2010).

14. C. M. Lalau-Keraly, S. Bhargava, O. D. Miller, and E. Yablonovitch, “Adjoint shape optimization applied to electromagnetic design,” Opt. Express 21, 21693–21701 (2013). [CrossRef] [PubMed]

15. A. A. Oberai, N. H. Gokhale, and G. R. Feijoo, “Solution of inverse problems in elasticity imaging using the adjoint method,” Inverse Probl. 19, 297 (2003). [CrossRef]

16. G. Bao and P. Li, “Inverse medium scattering for three-dimensional time harmonic Maxwell equations,” Inverse Probl. 20, L1 (2004). [CrossRef]

17. G. Bao and P. Li, “Inverse medium scattering problems for electromagnetic waves,” SIAM J. Appl. Math. 65, 2049–2066 (2005). [CrossRef]

18. A. Berk, “Variational principles for electromagnetic resonators and waveguides,” IRE T. Antenn. Propag. 4, 104–111 (1956). [CrossRef]

19. D. Jones, “A critique of the variational method in scattering problems,” IRE T. Antenn. Propag. 4, 297–301 (1956). [CrossRef]

20. T. G. Hazel and A. Wexler, “Variational formulation of the Dirichlet boundary condition,” IEEE T. Microw. Theory 20, 385–390 (1972). [CrossRef]

21. B. H. McDonald, M. Friedman, and A. Wexler, “Variational solution of integral equations,” IEEE T. Microw. Theory 22, 237–248 (1974). [CrossRef]

22. A. Konrad, “Vector variational formulation of electromagnetic fields in anisotropic media,” IEEE T. Microw. Theory 24, 553–559 (1976). [CrossRef]

23. K. Morishita and N. Kumagai, “Unified approach to the derivation of variational expression for electromagnetic fields,” IEEE T. Microw. Theory 25, 34–40 (1977). [CrossRef]

24. A. Mohsen, “On the variational solution of electromagnetic problems in lossy anisotropic inhomogeneous media,” J. Phys. A-Math. Gen. 11, 1681 (1978). [CrossRef]

25. W. K. H. Panofsky and M. Phillips, Classical Electricity and Magnetism (Addison-Wesley, 1962).

26. C. H. Chen and C.-D. Lien, “The variational principle for non-self-adjoint electromagnetic problems,” IEEE T. Microw. Theory 28, 878–886 (1980). [CrossRef]

27. G. Jeng and A. Wexler, “Self-adjoint variational formulation of problems having non-self-adjoint operators,” IEEE T. Microw. Theory 26, 91–94 (1978). [CrossRef]

28. W. C. Chew, Waves and Fields in Inhomogeneous Media (IEEE Press, 1995).

29. R. Kleinman and P. den Berg, “A modified gradient method for two-dimensional problems in tomography,” J. Comput. Appl. Math. 42, 17–35 (1992). [CrossRef]

30. K. Kilgore, S. Moskow, and J. C. Schotland, “Convergence of the Born and inverse Born series for electromagnetic scattering,” Appl. Anal. 96, 1737–1748 (2017). [CrossRef]

Name	Description
Visualization 1	Diffractive solution of a single incident field coupled to a single scattered field.
Visualization 2	Refractive solution without smoothing of a single incident field coupled to a single scattered field.
Visualization 3	Refractive solution with smoothing of a single incident field coupled to a single scattered field.
Visualization 4	Diffractive solution of two incident fields coupled to their respective scattered fields.
Visualization 5	Refractive solution with smoothing of two incident fields coupled to their respective scattered fields.

Inverse scattering with a non self-adjoint variational formulation

Abstract

1. Introduction

2. Deriving the functional for the non-self-adjoint Maxwell’s equations

3. Posing the inverse scattering problem

4. Varying the medium to find a solution to the inverse scattering problem

5. A demonstration of medium design

6. Conclusion

Funding

References and links

Supplementary Material (5)

Cited By

Figures (7)

Equations (39)

Optics Express