Upper bounds on absorption and scattering

A general framework for determining fundamental bounds in nanophotonics is introduced in this paper. The theory is based on convex optimization of dual problems constructed from operators generated by electromagnetic integral equations. The optimized variable is a contrast current defined within a prescribed region of a given material constitutive relations. Two power conservation constraints analogous to optical theorem are utilized to tighten the bounds and to prescribe either losses or material properties. Thanks to the utilization of matrix rank-1 updates, modal decompositions, and model order reduction techniques, the optimization procedure is computationally efficient even for complicated scenarios. No dual gaps are observed. The method is well-suited to accommodate material anisotropy and inhomogeneity. To demonstrate the validity of the method, bounds on scattering, absorption, and extinction cross sections are derived first and evaluated for several canonical regions. The tightness of the bounds is verified by comparison to optimized spherical nanoparticles and shells. The next metric investigated is bi-directional scattering studied closely on a particular example of an electrically thin slab. Finally, the bounds are established for Purcell's factor and local field enhancement where a dimer is used as a practical example.


Introduction
As the field of nanophotonics becomes more mature, interest is shifting away from the analysis of simple systems (uniform waveguides, spheres, rods, etc.) and toward the synthesis of structures with engineered electromagnetic behavior, e.g., maximal absorption [1], directional emission [2,3], directed scattering [4,5], fluorescence diplexing [6], waveguide power division [7], field confinement [8], and waveguide diplexing (wavelength splitters) [9]. In order to generate novel geometries optimized for these particular design objectives, computational inverse design is often employed, see, e.g., [10,11] and the references therein. While such methods excel in the exploration of extremely broad design spaces and the discovery of nonintuitive solutions, there nevertheless exists a strong need for analytic results which inform, direct, and truncate their computationally intensive calculations [10].
Physical bounds on performance objectives constitute a particular class of analytic results that aid inverse design in this way. For example, rather than searching for a design which achieves a particular objective value, one may instead employ inverse design tools to search for an optimal design with superior performance over all other possible structures satisfying some fixed constraints. In searching for this optimal design, the determination of physical bounds is critical in comparing realized performance to the theoretical optimum. Additionally, secondary analytic results and features of optimal designs may often be obtained through the process of deriving physical bounds, as in the study of optimal antennas operating in the microwave regime. These secondary results add understanding to the behavior of optimal designs and may provide valuable input to specialized inverse design algorithms, e.g., starting points or feature definitions.
Several approaches have previously been applied to calculate physical bounds on the performance of nanophotonic structures, though each of these methods has unique limitations. One approach is to analyze the broadband characteristics of scatterers using passivity-based sum rules [12]. Sum rules of this kind lead to important conclusions regarding many practical applications (e.g., cloaking [13]), though they fail to offer insight into performance bounds at single frequencies. Attempts to combine sum rules with other methods of calculating bounds rely on assuming simple frequency domain responses, e.g., Lorentzian line shapes [14]. Alternatively, analysis may be restricted to canonical structures, such as layered spheres ("core-shell structures"), where the low number of degrees of freedom allows for tight brute force optimization [15]. Dipolar interaction for electrically small structures is analyzed in [16], spherical mode expansions in [17], and effects of lossy background media in [18].
Based primarily on the optical theorem, shape independent bounds exist [7] assuming uniform fields and neglecting effects of electric size and shape. These bounds depend only on the relative volume of an object and its material properties. Shape independent bounds of this kind are necessarily loose and, in general, can only be approached by inverse design procedures in special cases. In the interest of driving optimal design, there is a need to tighten this class of bounds by making them dependent on both material and the shape of the allocated design region. Shape dependent bounds adapted in this way for maximum radiative heat transfer are given in [1,19]. Shape dependent bounds for maximum absorption are given in [20] and for photonic design in [21].
Unlike in photonics, the development of shape dependent bounds is fairly mature in the area of antenna theory, where the task of optimizing many important physical parameters may be cast as convex optimization problems in an unknown current distribution. Once the operator framework relating source currents to various metrics (e.g., gain, Q-factor, efficiency) is established, the development of new bounds reduces to the task of manipulating and combining optimization problems of canonical forms. Because of the convex (or relaxed) nature of such optimization problems [22], their solution is deterministic, many times reducing to the solution of a single eigenvalue problem. Expressing unknown currents in a finite dimensional basis, the finite dimensional matrix operators governing quantities of interest may be readily calculated using the tools developed for solving integral equation problems [23], i.e., using the method of moments [24]. Problems in antenna theory solved in this way include maximization of directive gain [25,26,27], maximization of radiation efficiency [28,26,29], minimization of Q-factor (analogous to maximizing bandwidth) [30,31], and the trade-off between these parameters [32,27].
The purpose of this paper is to transfer this general approach to the field of optics, where the fundamental electromagnetic physics remain the same as in classical antenna theory but the metrics of interest shift from antenna parameters to quantities related to scattering, absorption, extinction, and local field enhancement. To carry this out, we first develop a framework in which to discuss quantities of interest (e.g., scattered power, local fields) in terms of operators acting on current densities confined to a region of interest. This general methodology is shared by previous work [7], though here we use notation and nomenclature based heavily on the numerical solution of integral equations in antenna theory [24] which has been used extensively used for developing bounds in that area [23]. We then proceed by formulating optimization problems giving rise to shape-specific bounds on absorption, scattering, and extinction cross sections as well as local field enhancement (Purcell's factor) and directive scattering. We formulate each bound using two sets of constraints governing the enforcement of material constitutive relations. In a few certain cases, the bounds derived here replicate previous results, however we include them here to highlight the generality and flexibility of the approach taken. Numerical examples are calculated for each derived bound and serve to demonstrate salient features.

Physical components of problems
The main goal of this paper is to find bounds on various scattering and absorption metrics by formulating optimization problems which are either convex or can be relaxed to a convex form. Because many of the components and features are shared between several problems, we begin with an overview of the objectives and constraints used throughout the rest of the paper. Note that throughout the entirety of the paper, time harmonic steady state is assumed and all quantities are implicitly functions of angular frequency ω with exp(−iωt) dependence. The wavenumber k = ω/c 0 is always that of free space, with c 0 being the speed of light in vacuum. Horizontal axes in figures show free-space wavelength λ = 2π/k, which is in many cases normalized with respect to the radius a of the smallest sphere circumscribing the scatterer.

Objectives
The physical quantities investigated in this paper may be represented using linear or sesquilinear (quadratic) forms in a current density distribution J (r), with all integrations taken over the current's entire spatial support Ω ⊂ R 3 , see Fig. 1. Without loss of generality, we assume an appropriate basis {ψ n (r)} for the current density, e.g., is chosen in which these forms reduce to linear or quadratic forms involving vectors and matrices. For clarity (and compatibility with finite dimensional optimization tools) we only consider the vector/matrix forms except when their equivalence to the forms in (1) and (2) yields additional insight.
Here we briefly describe the physical interpretation and mathematical properties of the vector and matrix operators describing quantities used in optimization problems throughout this paper. The exact formulation of these operators are not given here but are readily available in the literature [33,24,23]. The following discussion assumes that the matrix or vector forms of all operators strictly follow the analytic features (e.g., definiteness) of their continuous counterparts.

Radiated power, absorbed power, and reactance
The total cycle-mean power radiated P r and absorbed P a by an arbitrary current I may be cast as quadratic forms [23] and respectively. Notice that rigorously the power P r refers to a cycle-mean scattered power since the current density I represents an equivalent contrast source situated in vacuum, see Fig. 1 and App. A for more details. As currents should not radiate negative power we have P r ≥ 0 for all currents and equivalently R 0 0. Similarly for currents in lossy media we have P a ≥ 0 and thus R ρ 0. Both operators are real symmetric, and their sum P r + P a represents the total cycle mean real power P t released by the current I, i.e., Employing spherical wave decomposition, the symmetric positive definite radiation operator R 0 may be constructed [34] as where the inner dimension of the matrix S depends on the number of spherical harmonics used. The number of necessary spherical harmonics, and thus the rank of R 0 , may be approximated using the electrical size of the current support under consideration. In general, electrically smaller support regions correspond to lower rank.
The reactance of a current, or equivalently the cycle mean difference in stored magnetic and electric energies, W m and W e , may be calculated via the quadratic form [35] 2ω (W m − W e ) = 1 2 I H XI.
The reactance matrix X is real symmetric and, in most cases, indefinite. Like the operator R, the matrix X may be decomposed into two components, one material dependent and one material independent, see Sec. 2.3. Together, the matrices R and X constitute the impedance matrix Z = R − iX, see Sec. 2.2 and App. A.

Radiation intensity
The scattered field componentê · E s (r) produced by the current I in ther = r/|r| direction at distance r → ∞ with polarizationê is represented by the linear form such that the corresponding partial radiation intensity is [36] U (r,ê) = 1 2 From the above expression we observe the rank-1 nature of the operator FF H , a feature that will enable certain problems to be solved in a computationally efficient manner. Note that the operator FF H depends directly on the observation direction and polarization. Like the total absorbed and radiated power operators, the operator FF H is positive semi-definite and Hermitian symmetric.

Radiation enhancement and Purcell's factor
The cycle mean power radiated by an electric dipole moment p in the presence of a scatterer is [37] where E s (r p ) is an electric field generated by the scatterer at the point of the dipole, P a is the cycle mean power lost within the scatterer (7) and is the cycle mean power radiated by the dipole in free space, with 0 being permittivity of vacuum (the assumed background medium).
Analogously to the electric far field in (11), the following linear form can be defined from which a Purcell's factor [38] can be defined as which assumes that the back reaction of the dipole to the scatterer may be neglected [39,40]. This Purcell's ratio characterizes radiation enhancement provided by the scatterer, a metric of primary importance in many areas of applied optics. Note that unlike in (16), Purcell's factor might also be defined without the subtraction of the loss term [40]. Although in the experimental characterization this might be the only option, evaluating bounds on Purcell's factor in the presence of realistic scatterers without subtracting losses leads to bounds that are too optimistic, as the Purcell's factor may be increased through absorption rather than radiation.

Real and reactive power constraints
An incident electric field E inc (r) may be projected onto the chosen basis in (3) via Such an excitation V uniquely induces a current distribution I obeying where Z is the impedance matrix representing the underlying integral operator mapping currents to fields [24], see App. A for more details. For the purpose of formulating general bounds on the behavior of currents existing with a given support, we relax this expression through testing with the current itself, i.e., which is an algebraic representation of power conservation [37, §6.9], see App. A for more details. Considering the excitation field V is given, two power constraints are introduced from (19) as representing the conservation of real (cyclic mean) and reactive power, respectively. The expression in (20) may be interpreted as a form of the optical theorem relating the total extincted power to that removed from the incident field [37, §10.11]. For plane wave incidence the equivalence to the optical theorem may be further used to relate the extincted power to a forward scattering amplitude.
As discussed in greater detail in Sec. 2.4, use of the constraint in (20) within an optimization over the current I amounts to enforcing both loss and radiation properties of an object when illuminated by the incident field V. For this reason we denote the use of this constraint as using "prescribed losses" (indicated throughout the paper by the subscript R). Similarly, inclusion of the constraint in (21) enforces material reactance properties, so together the use of (20) and (21) is denoted as using "prescribed materials" (indicated throughout the paper by the subscript Z).

Material properties and contrast current
Throughout this paper we consider isotropic 1 , non-magnetic materials with the constitutive relation where is the frequency dependent permittivity and χ is the frequency dependent susceptibility.
In the presence of the electric field E, the contrast (polarization) current density J may be written in terms of the complex resistivity ρ as where ρ has real and imaginary parts The real part of resistivity ρ r is related to a material figure of merit ζ = |χ| 2 / Im χ introduced in [7,20] as ρ r /a = η 0 /(ζka). This last equality can be used in Sec. 3.1 to directly compare the results obtained in this paper with the results obtained in [20]. Under the discretization in (3), the complex resistivity governs the nature of certain physical quantities and their associated operators. In general the radiation operator R 0 does not depend on the resistivity, while the loss operator depends on its real component ρ r , i.e., The situation is analogous for the reactance operator For homogeneous bodies, the matrices R ρ and X ρ scale linearly with the real and imaginary components of the resistivity, respectively.

Formulating bounds using optimal currents
In the remainder of the paper, we formulate bounds maximizing absorption, scattering, and extinction cross sections, directed scattering, and Purcell's factor for an obstacle with resistivity ρ = ρ r + iρ i confined to a region Ω, excited by a prescribed incident field. To do this, we construct and maximize relevant quantities over all possible current density distributions J supported in Ω, subject to certain constraints and fixed material parameters. These constraints (see (20) and (21)) enforce conservation of particular quantities (e.g., real power) but do not mandate that the optimized current be directly excitable by a given incident field. In this way, the formulated optimization problems yield bounds for all possible structures Ω ⊆ Ω supported within the region Ω consisting of either the structure medium (permittivity ) or the background medium (background permittivity 0 ), as shown in Fig. 2. Notice that the exclusion of a given subregion is realized in this paradigm by setting the contrast current density equal zero in that subregion. Because the optimization method is allowed to set currents within any such subregion in Ω to zero, all possible exclusions are automatically considered. This concept of searching for optimal currents has extensively been used in the study of bounds on antenna performance [23].
As stated previously, we may interpret the application of the real power conservation bound in (20) as enforcing only the prescription of material losses ρ r . The imaginary (reactive) part of the resistivity ρ i does not enter into this constraint and may be considered to be chosen freely. This implies that, under only the constraint of (20), the conservation of reactive power in (21) may always be satisfied by a suitable a posteriori choice of bulk reactivity ρ i . When both the conservation of real power (20) and of reactive power (21) are used, an optimization problem with two constraints is formed which is always more restrictive then the former, since now the reactance part of resistivity ρ i is prescribed prior to the optimization.
Application of both constraints in (20) and (21) effectively allows for the calculation of bounds on structures synthesized from fully known and characterized materials (such as gold). Though looser, application of only the real power bound (20), allows for the examination of hypothetical future materials which may have fixed losses but tunable bulk reactivity.

Maximization of cross sections
In this section we study upper bounds on the absorption, scattering, and extinction cross sections, denoted as σ a , σ s , and σ t , respectively, achievable by an arbitrary object constructed of material with resistivity ρ confined to the support Ω. In all cases, a fixed incident field V is assumed which, in the case of cross sections, corresponds to a plane wave [41,42]. Two forms of constraints are considered. We begin with the simpler of the two, where only the real part ρ r of the resistivity is prescribed (corresponding to constraint (20)). We then extend the problem to include both real and imaginary components of the resistivity (corresponding to both (20) and (21)), see Tab. 3.1 for corresponding section numbers for each problem.

Absorption
Maximization of the absorbed power (7) with prescribed losses over all possible current densities is determined by the solution to the optimization problem which can, with the help of Sec. 2.1.1, be written as a primal QCQP This QCQP is not convex [22] but can be analyzed using the techniques in App. B. Particularly, the maximum absorbed power is determined from solution to a dual problem

Constraints
Relative Problem Quantity (20) (21) Complexity Section maximum absorbed power Table 1: Optimization problems solved in the paper.
where the condition ν > 1 follows from the low-rank representation of the matrix R 0 , see App. B for more details. This dual problem can directly be solved using a line search algorithm such as the Newton or bisection algorithms [22]. Like many of the other problems detailed in later sections, solution to (29) can advantageously be formulated in a basis Q which simultaneously diagonalize the operators R ρ and R 0 . Vectors of this basis (columns I n of matrix Q) are called radiation modes, see App. C, and are the solution to a generalized eigenvalue problem [43] R 0 I n = n R ρ I n .
The radiation matrix R 0 is known to be of low rank, and can advantageously be decomposed as R 0 = S H S, see [34]. The loss matrix R ρ is a sparse full rank positive definite matrix and can be decomposed via a Cholesky decompostion [43] as R ρ = Υ H Υ. Substituting into (30) and using a singular value decomposition (SVD) [43] SΥ −1 = U 1 ΣU H 2 , where matrices U 1 , U 2 are unitary and matrix Σ is diagonal, it can be shown that n = Σ 2 nn and Q = Υ −1 U 2 . Assuming I = QĨ, the minimization problem (29), normalized to incident power flux S 0 , can be rewritten as The second equality in (31) assumes that the excitation vector has been composed of spherical waves as V = S H a/(k √ η 0 ), where a are dimensionless spherical wave expansion coefficients of an incident plane wave with unit amplitude, see App. C. A projectionã = U H 1 a has also been defined to ease the notation. As shown in App. B, the form in (31) (as compared to (29)) allows for simple analytical evaluation of the first and second derivative with respect to the Lagrange multiplier ν being well prepared for a minimization via Newton's algorithm. Furthermore, the summation involved in (31) converges quickly in realistic scenarios, since for small electrical sizes (ka ≤ 1), the eigenvalues n decay rapidly in n and only a few terms are needed, see App. C. Notice that a similar decomposition was implicitly used in [20,Eq. 4,5] to solve a related optimization problem.

Scattering
Maximization of the scattered power is determined by interchanging R ρ and R 0 in (28) and (29) which gives maximize with the dual problem where ν 1 = 1 /(1 + 1 ) with 1 the largest eigenvalue (30) for the radiation modes. Observe that this value of ν 1 coincides with the maximum radiation efficiency for antennas restricted to the region Ω [32].
Analogously to the bound on the absorption cross section in (31), the bound on the scattering cross section can be written in terms of radiation modes as

Extinction
The maximum extincted power is determined from maximization of the total cycle mean power P t . In this case, the power balance in (19) can be used to greatly reduce the complexity of the problem into maximize Re{V H I} with the explicit solution This expression is also recognized from bounds on the maximum effective area A e of antennas restricted to region Ω [27]. Following the same methodology as applied to absorption and scattering, the maximum extinction cross section (36) can be written in terms of radiation modes as

Electrically small scatterers
Although the bounds in (31), (34), and (37) are easily determined for an arbitrary shaped geometry, their explicit approximations depending only on resistivity ρ r , volume V , and free-space wavenumber k are of great interest [7]. Such approximation is possible in the limit of small electric sizes ka = 2πa/λ 1. In this case it can be assumed that the summations in (31), (34) and (37) are dominated by the first term. If only this dominant term (n = 1) is kept, an analytical solution exists in all three cases and reads where the dominant radiation mode 1 ≈ η 0 k 2 V /(6πρ r ) and the projection of the plane wave on a dipole mode |ã 1 | 2 = 6π for ka 1 have been used, see App. C. Relation (38) shows that the absorption cross section is, in the limit of small electric size, bounded by σ a,ρr introduced in [7] and by the dipole bound σ a,d [44]. The material-only bound is valid for all electrical sizes [7] and is obtained by neglecting the radiation term R 0 in (29). The dipole bound is restricted to dipole interaction and can be obtained from the directivity, D, and gain G, of a lossless Hertzian dipole D = G = 3/2 having effective area λ 2 G/(4π) = 3π/(2k 2 ) [45]. These two bounds intersect at k = 3πρ r /(2η 0 V ). The bound σ a,R has a maximum for 1 = 1 but decreases monotonically with increasing wavenumber k for electrically small objects. The single term (dipole) approximation is not valid as ρ r → 0 since the dipole contribution to the absorption is negligible and therefore it is necessary to include higher order terms.
Results for bounds on scattering cross section (39) show three distinct regions at small electrical sizes as compared with the two for the absorption cross section (38). The scattering cross section first increases with k 2 reaching its maximum σ s,ρr for the wavenumber after which it deceases with k −2 . Small self-resonant objects are often minimum scattering [46,47] and hence σ a = σ s which implies 1 = 1 which is recognized as (41) and from the maximum of (39). The small electric size limit of the bound on the extinction cross section (40) is similar to the two previous cases.

Numerical examples
Bounds on absorption, scattering, and extinction cross sections are depicted in Fig. 3 for obstacles supported in a spherical region with radius a and resistivities ρ r /a ∈ {0.01, 0.1, 1, 10} Ω. The bounds are computed using an analytical prescription for radiation modes of a sphere (88), normalized with the geometrical cross section A cross = πa 2 , and shown for electric sizes 10 −3 ≤ ka ≤ 10 2 . Compared with the small-size approximations from Sec. 3.1.4 (dashed and dotted lines) it is observed that σ a,R follows σ a,ρr [7] up to k = k s /2, cf., (41), where the small size approximations (38) intersect. Then the bound follows the dipole approximation until the structure starts to support higher order modes. The maximum of σ s,R is reached at k = k s and is well described by the small size approximation (39). The dipole approximations are most useful for low-loss cases where they approximate the bounds over a large range of frequencies.
The approximate onset of different radiation modes can be found from the condition n ≈ 1, with the radiation modes for the sphere in (88) depicted  (31), scattering (34), and extinction (37) cross sections for spherical regions compared with realized cross sections of solid spheres and spherical shells homogeneously filled with resistivity (24) ρ r /a = 1 Ω and optimized reactance ρ i and shell thickness.
the limit ρ r → 0. This increase is, however, slow for k > k s and negligible in the limit of electrically large structures where σ a,R approaches the geometrical cross section A cross and σ s,R , σ t,R approach 4A cross . Notice that the bound σ a,R is similar to that derived in [20]. The results shown in Fig. 3 for σ a,R can therefore be compared with those in Fig. 2 from [20] taking into account the relation ρ r /a = η 0 /(ζka). Fig. 4, where the realized cross sections for solid spheres and spherical shells are compared with the bounds. The imaginary part ρ i = Im{ρ} and shell thickness are swept numerically for each electric size and the maximum realized cross sections are depicted in Fig. 4. Solid spheres reach the bounds for small electric sizes ka ≤ 0.1 but are slightly below the bounds for larger sizes. Optimization over the layer thickness increases the cross sections slightly but does not reach the bounds.

Tightness of the bounds is investigated in
Bounds on absorption and scattering cross sections for obstacles circumscribed by spheroidal regions with semi axes a r and a z and an incident plane wave from theẑ-direction are depicted in Fig. 5. The overall behavior of the bounds for circumscribing spheres in Fig. 3 and spheroids are similar. A decrease for electrically small spheroid as compared with the sphere a r = a z is explained by the reduced volume, see (38) and (39).

Prescribed materials
Here we adapt the physical bounds in Sec. 3.1 to a stricter form where both real and imaginary components of the resistivity ρ are fixed.

Scattering
The scattering cross section is similarly obtained by interchanging R 0 and R ρ giving the primal problem maximize I H R 0 I subject to I H RI − Re{I H V} = 0 and the dual problem where D = {(ν, µ) : −R 0 + νR + µX 0}.

Extinction
For the extinction cross section, the optimization problem with both constraints is yielding the dual function where D = {(ν, µ) : νR + µX 0}. Substituting µ = νµ 1 and carrying out the minimization over ν (now entirely outside of the matrix inversion) gives ν = ±1/ 1 + µ 2 1 , which simplifies (47) to P t,Z = min Since the optimization problem (48) only involves two matrices, it can still be, analogously to problems in Sec. 3.1, diagonalized using a basis Q with basis vectors (columns I n of matrix Q) generated by an eigenvalue problem These eigenmodes resemble characteristic modes [35], though here, unlike the classical formulation for characteristic modes on perfectly conducting bodies, the operator R contains terms related to both loss and radiation and thus alters the modal set's orthogonality properties [48].
Normalizing Q H RQ to be an identity matrix, Q H XQ is a diagonal matrix of eigennumbers λ n and the dual formulation (48) can be rewritten as P t,Z = min whereṼ = Q H V and the + sign corresponds to the domain D 1+ while thesign corresponds to the domain D 1− . These domains depend on the definiteness of the matrix X and are defined as In (50), the inner minimization is performed first for the multiplier µ 1 in the domain D 1+ and D 1− separately and from those two results, the minimum is selected by the outer minimization.  Bounds on the extinction cross section are depicted in Fig. 6 for obstacles composed of gold (Au), see App. D, and circumscribed by a sphere with radius a ∈ {10, 20, 50, 100} nm. The bounds (40), (37), and (50) are compared with the realized extinction cross sections of solid gold spheres and spherical gold shells where the shell thickness is optimized to maximize the extinction cross section σ t . For all radii it is clear that inclusion of the reactance constraint (21) has a large effect for longer wavelengths or, more accurately, electrically small sizes (ka = 2πa/λ ≤ 1). The effect diminishes as the electric size increases and is negligible for ka ≥ 1, cf., [27]. Starting with radius a = 10 nm, it is observed that σ t,ρr ≈ σ t,R . Inclusion of the reactance constraint (21) reduces the bound on σ t by orders of magnitude for wavelengths λ ≥ 1 μm. The differences are smaller for shorter wavelengths and minuscule for λ ≈ 0.2 μm and λ ≈ 0.5 μm. The longer wavelength 0.5 μm corresponds to the dipole plasmonic resonance [49,41] which occurs for a relative permittivity with a real part close to −2 (Fröhlich condition), see Fig. 6, but the high losses weaken this effect. The extinction cross section for homogeneous spheres is also close to the σ t,Z bound for λ ∈ [0.2, 0.5] μm but far from the bound for longer wavelengths. The performance improves for spherical shells with shell thickness optimized for maximum σ t which are seen to follow the σ t,Z bound for the considered wavelengths. This is partly explained by the added degree of freedom from the shell thickness which is used to tune the plasmonic resonance at each wavelength. The conclusions are similar for the larger radii but σ t,ρr and σ t,R starts to differ and the difference is large for a ≥ 50 nm. Performance of the solid spheres and spherical shells relative to the bounds worsen with increasing radius.

Numerical examples
Dielectric materials (such as silicon, silica, or germanium) offer much lower losses than gold and could hence provide larger cross sections as suggested from the bounds in Fig. 3. In Fig. 7, bounds (37) and (50) on the extinction cross section are depicted for obstacles supported in a spherical region with radius a and composed of a material with relative permittivity Here, we observe that σ t,Z follows σ t,R for electrically larger ka ≥ 1 cases but has a sharp drop around ka = 1. Effects of the losses Im{ r } are also relatively small for ka ≥ 1 and negligible for ka < 1.
Contour plots depicting the upper bound (50) on the normalized extinction cross section for obstacles supported within a sphere with radius a ∈ [0.05, 0.2] μm and made of gold (Au), silver (Ag), and a low loss dielectric ( r = 11 + i10 −5 ) are shown in Fig. 8. The results show which

Directional scattering
Here we study bounds on scattering into a particular direction and polarization by a scatterer illuminated by an arbitrary incident wave. More specifically, let the excitation be described by vector V and let the scattering be measured by the radiation intensity U (r,ê) in the directionr with polarizationê, see Sec. 2.1.2. An example of this scenario where the incident field is a plane wave is depicted in Fig. 9. In a complete analogy to Sec. 3.1 and Sec. 3.2 an optimization problem is formulated for upper bounds on radiation intensity U (r,ê) using the power constraints (20) and (21). The first constraint would be used when only real part of material resistivity is prescribed, while both constraints should be used to enforce both real and imaginary components of the complex resistivity. Like absorption, scattering, and extinction studied in Sec. 3, under either selection of constraints, the upper bound on radiation intensity U (r,ê) may be interpreted in terms of an upper bound on bistatic radar cross section [36] σ bis = 4πU (r,ê) x y Ω Figure 9: Schematic of a directed scattering optimization problem. Currents within the structure Ω are optimized, subject to material-related power balance constraints, to maximize the scattered fields in a particular polarization and direction.
where the incident field is a plane wave with power flux S 0 .

Prescribed losses
In the case of resistivity ρ i being free, maximization of radiation intensity is determined by the solution to a QCQP maximize I H FF H I subject to I H RI − Re{I H V} = 0.
The associated dual problem reads where ν 1 = F H R −1 F. Since the matrix FF H is of rank 1, an analytical solution to (54) exists, i.e., the Sherman-Morrison-Woodbury identity [43] may be used to evaluate the matrix inverse analytically and the simplified scalar optimization problem may be solved in closed form. The result of this procedure is given by where

Prescribed materials
If the constraint on reactive power conservation (21) with the accompanying dual problem where D = {(ν, µ) : −FF H + νR + µX 0}. Analogously to Sec. 3.2.3, substitution of µ = νµ 1 can be used to eliminate the Lagrange multiplier ν. This, together with application of the Sherman-Morrison-Woodbury identity, yields the simplified optimization problem U Z = min where α, β, and γ are the same as in (56) except for the added dependence on the Lagrange parameter µ 1 via the altered form of the matrix G = (R + µ 1 X) −1 . The domain D 1 is given by the union of domains D 1+ and D 1− from (51), i.e., The inverse in G may be conveniently factored using the characteristic-mode-like eigenvalue problem (49) such that each term in (56) may be calculated efficiently using where {λ n } are again the eigenvalues of the characteristic-mode-like problem in (49) and the tilde denotes the projectionã = Q H a.

Numerical examples
As a simple example, consider the directive scattering bounds for the rectangular prism (slab) in Fig. 9. The object's aspect ratio x : y : z is 10 : 5 : 1. Illumination is aŷ-polarized plane wave incident from theẑ direction. Observation directions are restricted to the xz plane and described by the angle θ. Two observation polarizations are selected: parallel ( , electric field y-directed) and perpendicular (⊥, electric field in the xz plane). First, we calculate the bounds for directive scattering with prescribed losses for slabs of three electrical sizes ka = {0.01, 0.1, 1} and two real resistivities ρ r /a = {0.01, 1} Ω with only the real power constraint enforced. Results in Fig. 10, presented in terms of bistatic radar cross sections, show similar trends as observed for scattering cross sections in Fig. 3, with nonmonotonic shifting of the overall maximal scattering response over all observation angles as a function of electrical size. We observe that both endfire (θ = π/2) and broadside (θ = 0, 2π)  Figure 10: Upper bounds on directed scattering, plotted in terms of normalized bistatic radar cross section σ bis /(πa 2 ), for the scattering geometry in Fig. 9 with prescribed losses. may serve as the direction with highest scattering potential, though for small electrical sizes the angular variation in maximal scattering disappears as the optimal current becomes a uniform dipole moment. Using (59), we calculate the directive scattering bounds for the same slab with prescribed material properties corresponding to gold. Absolute units are required in this case, so the thickness of the slab is set at 40 nm and three wavelengths λ = {470, 550, 665} nm within the optical regime are examined. The difference in wavelength in this case is small, nonetheless the dispersive properties of gold lead to interesting variation in the maximal scattering properties, shown in Fig. 11. We observe that over this range, the parallel polarization ( ) transitions from a maximum in the endfire direction (θ = 0, λ = 665 nm) to a maximum in the broadside direction (θ = π/2, λ = 470 nm). Additionally, in the parallel polarization the effect of the loss versus material constraints is minimal, as also observed in the trends for extinction cross section on large spherical scatterers in this range of frequencies, cf., Fig. 6. This is not the case, however, for the perpendicular polarization, where the inclusion of the reactance constraint impedes endfire scattering at all three wavelengths.

Purcell's factor
Upper bounds on enhancement of radiation from a point electric dipole (Purcell's factor, see Sec. 2.1.3) are studied in this section for the case of prescribed losses and for the case of prescribed materials.

Prescribed losses
Assuming a fixed excitation vector V which is the projection (17) of the electric field generated by the dipole p in free space onto the current basis (3), a maximum Purcell's factor (16) with prescribed material losses is generated by with the corresponding dual problem where ν 1 = −1/(1 + 1 ) with 1 being the largest radiation mode eigenvalue generated by (30). Employing diagonalization of the underlying matrices, the dual problem can be, similarly to Sec. 3.1.1, written as A sweep of the upper bound to Purcell's factor F R over frequency is shown in Fig. 12 for a gold spherical dimer together with the realized Purcell's factor F for the same structure. A  Figure 12: Frequency dependence of Purcell's factor (16) and its upper bounds (64), (66) for a spherical dimer with sphere radius r = 30 nm, separation distance between the spheres d = 0.4r and an electric dipole moment centered between the spheres and oriented along the axis of the dimer. Gold has been used for the dimer, see App. D. spherical geometry is used as a prototype of a core-shell particle which is commonly used for the metal-enhanced fluorescence [50,51]. In these cases, either the core or the shell are made of a plasmonic material (e.g., gold). Within the optimization paradigm used here, the contrast current chooses the preferred option.
Comparison of both lines in Fig. 12 shows that in the real scenario, the dimer structure is effectively excited only in the vicinity of plasmonic resonance of the dimer and even there the realized Purcell's factor is one order in magnitude lower than upper bound F R . Away from the plasmonic resonance a simple spherical dimer is far from acting as an optimal structure for radiation enhancement, as evidenced by a considerable gap between the realized Purcell's factor and its upper bound.

Prescribed material
The addition of a reactive power constraint to (62) results in a maximal Purcell's factor under prescribed material given by the optimization problem and dual problem where D = {(ν, µ) : R ρ + νR + µX 0}. The solution to this optimization problem must be performed by general purpose solvers, see App. B. The frequency dependence of the upper bound F Z is depicted in Fig. 12. As compared to upper bound F R , the additional constraint on the conservation of reactive power considerably tightens the bound at long wavelengths (small electrical sizes), pushing it closer to the realized Purcell's factor F . The noticeable gap between the realized Purcell's factor and its bound nevertheless still exists away from the plasmonic resonance. Thus, away from plasmonic resonance the realized structure would have to be modified in topology in order to achieve resonance coupling to the exciting dipole.

Conclusions
In this paper we have laid out a general optimization framework for determining bounds on several metrics related to electromagnetic scattering, absorption, and field enhancement. In all cases, bounds were formulated using contrast current density representing the optimization domain and matrix operators mapping current distributions onto physical quantities. Techniques such as modal decomposition and rank-1 inverse updates were applied to simplify each problem based on its distinguishing features. Through that process, many of the problems studied here may be solved in computationally efficient ways.
The framework discussed here is not limited to the problems studied in this paper. Further bounds on both near-and far-field metrics may be derived given the ability to calculate the necessary matrix operators. Though in this paper all scattering objects were considered to exist in vacuum, no substantial modification is necessary for the case of objects suspended within any other lossless dielectric background. Similarly, material inhomogeneity and anisotropy may be introduced with minimal technical changes.
The bounds presented in this paper represent the absolute optimal behavior achievable by objects of prescribed material and bounding support. In addition to providing insight into the fundamental physical limitations of specific electromagnetic processes, these bounds will serve as benchmarks for future topology optimization implementations. There, the associated optimal current solutions may also find utility in accelerating and informing the chosen topology optimization algorithm.
A complex power balance in this scenario is typically formulated as [37] − see Fig. 1 for definition of the region V and its boundary ∂V . Using a volume equivalence principle [52], the same scattering problem can be modeled by a contrast current density J which replaces the scatterer and is the source of scattered fields. In this point of view, the complex power balance is more naturally stated as where sources J + J i radiate in free space.
Comparing (67) and (68) identifies the term with a complex power flow within the material of the scatterer. The last equality results from a current expansion (3) and definitions introduced in Sec. 2.3 and defines the matrix Within the volume equivalence, the scattered field is itself a valid solution to Maxwell's equations in free space and therefore generates a complex power balance relation where defines matrix R 0 − iX 0 introduced in Sec. 2.2 as a projection of scattered field operator onto a current basis (3) see [24] for more details. Addition of (69) and (72) together with definition (17) generates the complex power constraint (19).

B Quadratically constrained quadratic programs
Assume an optimization problem of the form ] . The dual function [22,53] for this problem is where the Hessian matrix of the Lagrangian (strictly speaking, twice the Hessian matrix)

B.1 Single quadratic constraint
If C, c, c 0 = 0 in the optimization problem (75), considerable simplifications may be made. For example, assume that a matrix Q exists that simultaneously diagonalizes the matrices A, B.
The optimization problem (75) can then be recast into an alternative ("Q") basis as which in many cases allows for an analytical solution to the problem (84) that is given by the root of (79).
For two important cases considered in this paper, an analytical solution to (84) can be found. The first case occurs when A = 0 and the matrix B is of full rank and not indefinite. The Qbasis is in that case formed by an eigenvalue decomposition BQ = QD and the Lagrange multiplier which minimizes the dual function is given by under the condition that ν 2 opt > 0. The positive root of (86) is chosen whenB nn > 0, while the negative root of (86) is chosen whenB nn < 0.
The second case is realized when the matrix A = αL H 1 L 1 is a rank-1 Hermitian matrix, the matrix B = βL H 2 L 2 is Hermitian and of full rank, the vector a = 0 and b 0 = 0. The Q-basis is most easily obtained via Q = L −1 2 U 1 , where U 1 is a unitary matrix coming from the singular value decomposition L −H 2 L H 1 = U 1 ΣU H 2 . Since the matrix A is rank-1, only one element of the diagonal matrixÃ is non-zero. Denoting its index as i = 1, the solution to this special case reads where the correct sign in (87) is chosen according to the negative definiteness of the Hessian matrix (77).

C Radiation modes
Radiation modes are determined from the eigenvalue problem (30). Although radiation modes are easily determined numerically for arbitrary shapes using the S matrix [34], analytic expressions are valuable. The radiation modes for a solid sphere with homogeneous resistivity ρ r can where u (1) υ denote regular spherical waves [56,42], R τ,l = R (1) τ,l (ka) radial functions [56,34], and υ = (τ, s, m, l) a multi-index with s, l ∈ {1, 2}, m ∈ {0, . . . , l}, and l ∈ {1, . . . } ordered as n = 2(l 2 + l − 1 + (−1) s m + τ . Amplitudes of radiation modes for a homogeneous sphere with radius a are depicted in Fig. 13. The radiation modes are dominated by the three transverse magnetic (TM) (τ = 2) dipole modes (l = 1) for electrically small spheres followed by the three transverse electric (TE) (τ = 1) dipole and five TM quadrupole (l = 2) modes. They increase monotonically with the electrical size ka. The sum of radiation modes for a homogeneous object with volume V and resistivity ρ r are determined from the trace of matrix R 0 to n = k 2 η 0 V /(2πρ r ). Radiation modes can also be determined analytically in the limit of electrically small objects. In this limit, electric dipoles (τ = 2 and l = 1) dominate radiation and the Rayleigh quotient related to the eigenvalue problem (30) simplifies to where we used the low-frequency expansion of the regular spherical waves u (1) υ ≈ê/ √ 6π for electric dipoles in theê ∈ {x,ŷ,ẑ}-directions.
For inhomogeneous objects with resistivity ρ r (r), we use variational calculus and set J = J 0 +δJ 1 and evaluate for small perturbations δJ 1 with δ → 0 giving The functional is stationary for all perturbations J 1 if Ωê · J 1 dV = 0 and Ω ρ r J 1 · J 0 dV = 0 (91) for all J 1 which imply a current density J 0 ∼ê/ρ r . Within this approximation, the first three dominant radiation modes have eigenvalues for ka 1 and the last equality is for homogeneous objects with volume V . Radiation modes for a spherical and oblate and prolate spheroidal regions with semi-axes a r and a z are depicted in Fig. 14. The first three modes radiates as electric dipoles and have eigenvalues approximately given by the volume V = 4πa 2 r a z /3. This explains the decrease of the first modes compared to the sphere and ordering of the curves. Eigenvalues for the sphere appear in groups with 2l + 1 elements, where l denote the order of the spherical mode, see also (88) and Fig. 13. The eigenvalues decrease rapidly with the mode index such that higher order become negligible for l ka with n = 2l(l + 2) and the order, l, can often be estimated from l = ka + 7 3 √ ka + 3 [57]. Radiation modes together with an expansion of the incident plane wave in spherical waves [42] a υ = 4πi l−τ +1ê · Y υ (k), are used in (31), (34), and (37) to express the bounds on the cross sections in a form that is suitable for small size approximations. Here, Y υ denotes spherical harmonics [42] and the explicit value |ê · Y υ | 2 = 3/(8π) for the dipole terms (l = 1) is used in (38), (39), and (40).  Figure 14: Radiation modes n for spheroids with semi-axes a z and a r for ka = 1 with a = max{a z , a r } and resistivity ρ r normalized with η 0 /(ρ r k).

D Material models
A MATLAB implementation [58] of a Drude-Lorentz model of gold and silver [59] is used for the permittivity and resistivity ρ throughout the paper. The frequency dependence of both parameters in the interval λ ∈ [0.2, 12] μm is depicted in Fig. 15 and Fig. 16. This permittivity model also agrees well with older reference data for gold from [60]. Bulk material parameters are assumed throughout the paper. Surface effects resulting in a non-local material response are neglected.