PHYSICAL REVIEW RESEARCH 2, 033172 (2020) Global T operator bounds on electromagnetic scattering: Upper bounds on far-field cross sections

We present a method based on the scattering T operator, and conservation of net real and reactive power, to provide physical bounds on any electromagnetic design objective that can be framed as a net radiative emission, scattering or absorption process. Application of this approach to plane-wave scattering from an arbitrarily shaped, compact body of homogeneous electric susceptibility χ is found to predictively quantify and differentiate the relative performance of dielectric and metallic materials across all optical length scales. When the size of a device is restricted to be much smaller than the wavelength (a subwavelength cavity, antenna, nanoparticle, etc.), the maximum cross-section enhancement that may be achieved via material structuring is found to be much weaker than prior predictions: the response of strong metals (Re [χ ] 0) exhibits a diluted (homogenized) effective medium scaling ∝|χ |/ Im [χ ]; below a threshold size inversely proportional to the index of refraction (consistent with the half-wavelength resonance condition), the maximum cross-section enhancement possible with dielectrics (Re [χ ] > 0) shows the same material dependence as Rayleigh scattering. In the limit of a bounding volume much larger than the wavelength in all dimensions, achievable scattering interactions asymptote to the geometric area, as predicted by ray optics. For representative metal and dielectric materials, geometries capable of scattering power from an incident plane wave within an order of magnitude (typically a factor of two) of the bound are discovered by inverse design. The basis of the method rests entirely on scattering theory and can thus likely be applied to acoustics, quantum mechanics, and other wave physics.

Nevertheless, despite these apparent challenges, computational (inverse) design techniques based on local gradients have proven to be impressively successful [43][44][45], offering substantial improvements for applications such as on-chip optical routing [45][46][47], meta-optics [48][49][50], nonlinear frequency conversion [51,52], and engineered bandgaps [53,54]. The widespread success of these techniques, and their increasing prevalence, begs a number of questions. Namely, how far can this progress continue, and, if salient limits do exist, can this information be leveraged to either facilitate future inverse approaches or define better computational objectives. Absent benchmarks of what is possible, precise evaluation of the merits of any design methodology or algorithm is difficult: failure to meet desired application metrics may be caused by issues in the selected objectives, the range and type of parameters investigated, or the formulation itself.
Building from similar pragmatic motivations, and basic curiosity, prior efforts to elucidate bounds on optical interactions, surveyed briefly in Sec. III, have provided insights into a diverse collection of topics, including antennas [9,10,55], light trapping [6,[56][57][58][59][60], and optoelectronic [61][62][63][64] devices, and have resulted in improved design tools for a range of applications [39,[65][66][67]]. Yet, their domain of meaningful applicability remains fragmented. Barring recent computational bounds established by Angeris, Vučković and Boyd [39], which are limits of a different sort, applicability is highly context dependent. Relevant arguments largely shift with circumstance [68,69], and even within any single setting, attributes To what extent does the specification of an electric susceptibility (χ ) and a confining spatial domain for the design of a structured body (optical device), determine the maximum absorbed P abs , scattered P sct , and extinguished (total) P ext = P sct + P abs power that can be extracted from a known incident field? widely recognized to affect performance (e.g., differences between metallic and dielectric materials, the necessity of conserved quantities, required boundaries, minimum feature sizes) are frequently unaccounted for, leading to unphysical asymptotics and/or bounds many orders of magnitude larger than those achieved by state-of-the-art photonic structures, Fig. 1.
In this paper, we exploit the requirement of global (net) power conservation as constructed from the defining relation of the scattering T operator, in conjunction with Lagrange duality, to derive bounds on any electromagnetic objective equivalent to a net extinction, scattering, or absorption process. With minor modifications for cases where one is interested in only a portion of the output field [70], such phenomena encompass nearly every application mentioned above. The scheme, capturing the potential of any and all possible structures under the fundamental wave limitations contained in Maxwell's equations, requires only three specified attributes: the material the device will be made of, the volume it can occupy, and the source (current or electromagnetic field) that it will interact with. Directly, the conservation of real power is seen to set an upper bound on the magnitude of a system's net polarization response, while the analog of the optical theorem for reactive power, introducing the polarization phase, is shown to severely restrict the conditions under which resonant response is attainable, particularly in weak dielectrics (glasses), leading to significantly tighter limits compared to related works [71,72].
The utility of a more robust, methodic, approach for treating electromagnetic performance limits has been recognized by several other researchers (especially in the field of radio frequency antennas [70,73,74]), and in particular, concurrent works by Kuang et al. [72] and Gustafsson et al. [75] have converged on the same basic optimization approach given in Sec. II. Although Ref. [72] considers only the conservation of real power, and thus overestimates achievable performance for certain parameter combinations (Figs. 2 and 3), the findings presented in these articles are excellent complements to our results, further highlighting the adaptability and utility of Lagrange duality and scattering theory for predicting possible performance in photonics. Moreover, during review of this paper, another program for calculating bounds via the selfconsistency of the total scattering field has also been put forward by Trivedi et al. [76]. A brief description of this work, as well as those cited above, is given in Sec. III.
Application of the technique to compute limits on far-field scattering cross sections for any object of electric susceptibility χ that can be bounded by (contained in) either a ball of radius R or a periodic film of thickness t interacting with a plane wave of wavelength λ, codifies a substantial amount of intuition pertaining to optical devices. As R/λ → 0, the requirement of reactive power conservation means that the energy transferred between a generated polarization current and its exciting (incident) field, averaged over the confining volume, can never scale as the material enhancement factor of ζ mat = |χ | 2 / Im [χ ] introduced and broadly discussed in prior works [67][68][69]71,[77][78][79]. Instead, for metals (Re [χ ] −3), with Re [χ ] = −3 corresponding to the localized plasmon-polariton resonance [80] of a spherical nanoparticle, the relative strength of such interactions cannot surpass, as compared to ζ mat , the reduced form factor 3|χ |/ Im [χ ], consistent with a "dilution" of metallic response implied by homogenized or effective medium perspectives [81][82][83]. When Re [χ ] > −3, even this level of enhancement is not possible. For dielectrics (Re [χ ] > 0), enhancement is limited by the same material dependence that appears in Rayleigh scattering [84], approximately 3|χ |/|χ + 3|. (In either case, comparison with cross section limits based on shape dependent coupled mode or effective medium models [85,86] reveals notable inexactness.) After surpassing a wavelength condition inversely proportional to the index of refraction, the importance of reactive power is superseded by the conservation of real power, causing Re [χ ] to have less drastic consequences on the magnitude of achievable cross sections; for R ≈ λ/2, limits for dielectrics meet or surpass those of metals for realistic material loss values, Im [χ ]. In the macroscopic limit of R λ or t λ/2, the selected material plays almost no role in setting achievable scattering cross sections, other than determining the level of structuring that would be required, and the predictions of ray optics (geometric cross sections) emerge. Additionally, beyond these asymptotic features, for representative metallic and dielectric material selections, geometries discovered through topology optimization [43,44] are shown to come within an order of magnitude of the bounds for domain sizes varying between R/λ = 10 −3 and R/λ = 1, with connate agreement for periodic films, providing supporting evidence that the bounds are nearly tight and predictive.
These same findings also shed light on a range of fundamental questions, such as limitations for light extraction and trapping efficiency, and the relative merits of different materials for particular applications [64,83,87], and provide a much more quantitative perspective on which aspects of a design are most critical to device performance than prior analyses. Shortly, we foresee extensions of the framework to embedded sources and user defined design (containing) volumes as providing a means of formalizing, comparing, and contrasting existing paradigms within photonics, revealing 033172-2 FIG. 2. Scattering cross-section bounds for a plane wave incident on a compact object. The four panels show different aspects of scattering cross section bounds for any structure of electric susceptibility χ that can be contained within a spherical boundary of R/λ (e.g., cavities, nanoparticle, etc.). Dashed lines result from imposing the conservation of real power, as in Ref. [72]. Full lines result from additionally requiring the conservation of reactive power, as in Ref. [75]. As R → 0, limit values agree with (38) in all cases, and with (40) so long as Im [χ ]/| Re [χ ]| 10 −4 . Dots in panels (a) and (b) mark scattering cross sections achieved in actual geometries discovered by inverse design, for χ = −10 + i 10 −1 and χ = 16 + i 10 −6 , respectively. For the metal structure in (a), vertically aligned cross-hatched dots result from enforcing binarized permittivity profiles. (Binarization does not alter the dielectric results.) Two sample structures are shown as insets, with the plane wave incident from the more solid side of both designs, left side in (a), right side in (b), and aligned along the left-right symmetry axis. Unbounded cross sections are encountered only for fictitious metals with vanishing material loss, with the bounds exhibiting a weak, sublogarithmic, divergence as Im [χ ] → 0. More practically, cross-section enhancements surpassing ≈200 should not be expected. Descriptions of other major features are given in the accompanying text.
The paper is divided into four main sections. Section II begins with an overview of the T operator relations governing absorption, scattering and radiative processes, followed by a statement of the wave constraints and relaxations that are 033172-3 FIG. 3. Absorption cross-section bounds for a planewave incident on a compact object. The four panels provide analogous information as Fig. 2, but for absorbed rather than scattered power. Again, the dashed lines represent bounds determined by insisting only that real power is conserved, the solid lines result from additionally requiring the conservation of reactive power, and dots correspond to actual structures discovered by inverse design. For the example metal structure shown in (a), the planewave is incident from the solid side of the half shell (from the left in the first view and from the back left in the second). For the example dielectric structure, the planewave is incident from the left along the axis of symmetry of the mushroom cap in the full view, and from the back left in the cross cut. Small radii features are found to be in good agreement with the asymptotic predictions of (38) and (40) for Im χ/|Reχ | 10 −4 . When the confining volume is small, R/λ 1/20, the ability of a structured object to absorb radiation is found to be substantially weaker than past predictions [68,71,77]. As R → ∞, the geometric cross-section limit predicted by ray optics is recovered regardless of the material considered. Like scattering, these results indicate that absorption cross section enhancements larger than ≈ 200 should not be expected. explored in the rest of the manuscript. From these preliminaries, the calculations of limits is then cast in the language of optimization theory, and a solution in terms of the Lagrangian dual is given. In Sec. III, a brief examination of the similarities and differences of this approach with current art is supplied. Next, in Sec. IV, basic computational mechanics are examined, and the quasistatic R/λ → 0 scaling factors quoted elsewhere are derived based on analytically tractable single-channel asymptotics. Finally, Sec. V provides sample applications of the method as described above. This is likely 033172-4 the section of the text of most relevance to the majority of readers. Although only single frequency examples are given, broad bandwidth objectives should present no major hurdles [98,99]. Further technical details, necessary only for reproducing our results, appear in Appendix. Because the approach relies exclusively on the validity of scattering theory, it is likely that counterparts of the presented findings exist in acoustics, quantum mechanics, and any other wave physics.

II. FORMALISM
The key to the approach of this paper rests on the use of partial relaxations [100]. Past electromagnetic limits, including our own work, have been predominately formulated by placing bounds on the individual physical quantities entering an objective and then deducing a total bound by composing the component limits [101,102]. We begin, alternatively, with the fundamental scattering relation that any physical system must satisfy, derive algebraic consequences of this relations (e.g., energy conservation) and then suppose a subset of theses conclusions as constraints on an otherwise abstract optimization. This general procedure is detailed below. The formulas given in the Power objectives and Scattering constraints sections are exact. Relaxations (omissions of certain physical requirements) begin only after the Relaxation and optimization heading.

A. Power objectives
Considerations of power transfer in electromagnetics typically belong to one of two categories: initial flux problems, wherein power is drawn from an incident electromagnetic field, and initial source problems, wherein power is drawn from a predefined current excitation. Initial flux problems are typical in scattering theory, and as such, our nomenclature follows essentially from this area [103]. That is, we will denote the initial (incident, given, or bare) field with an i superscript (either |E i or |J i ) and the total (or dressed) field with a t superscript. For a pair of initial and total quantities referring to the same underlying field, the scattered field, s superscript, is defined as the difference |F s = |F t − |F i . There is a certain appeal to transforming one of these two classes of problem into the other via equivalent fields. However, due to the additional back action that can occur in initial source problems, in our experience a unified framework promotes logical slips. For this reason, the total polarization field of an initial flux problem (or total electromagnetic field of an initial source problem) will be referred to as a generated field, g superscript. With this notation, scattering theory for initial flux and source problems consists of the following relations: Here and throughout, G 0 marks the background or environmental Green's function, which may or may not be vacuum [80]. The V operator refers to the scattering potential (susceptibility) relative to this background (whatever material was not included when G 0 was computed), and |E i and |J i are similarly defined as initial fields (solutions) in the background. The remaining quantities in (1) and (2) are the impedance of free space Z, the wavenumber k = 2π/λ (with λ the wavelength), and the T operator, defined by the formal equality I = (V −1 − G 0 )T [69]. The three primary operator forms for energy transfer in an initial flux problem are the extracted power, the absorbed power, and the scattered power, with S E = |E i E i | denoting projection of the corresponding operators onto the incident fields. Reciprocally, the three principal forms characterizing power flow from an initial current excitation are the extracted power, the radiated power, and the material (loss) power, with S J = |J i J i | denoting projection of the corresponding operators onto the initial current sources. The naming of the final two forms, which appear less frequently than the other four, follows from the observation that once the total source |J t is determined the corresponding electromagnetic field is generated exclusively via the background Green function. Hence, the energy transfer dynamics of a total source are exactly those of a special "free" current distribution. Because the only pathway for power to flow from a current source in free space (or lossless background) is radiative emission, P rad src must be interpreted in this way-energy transfer into the source free solutions of the background-and similarly, P mat src must be equated with loss into the scatterer. This reversal of forms and physics (compared absorption and scattering) is sensible from the perspective of field conversion. Absorption is the conversion of an electromagnetic field into a current, and radiative emission the conversion of a current into an electromagnetic field. Scattering in an initial flux setting is the creation of a new field of the same type, as is material loss in an initial source setting.
Note, however, that there is a caveat to these interpretations. G 0 , as we have defined it, accounts for the entire electromagnetic background. As such, in non vacuum cases, any polarization current associated with a background solution will not appear in any |J field, generally comingling the corresponding physical interpretation of the various power quantities. For instance, as Asym [G 0 ] describes power flow into the full homogeneous solutions of Maxwell's equations, if the environment for which G 0 is determined contains absorptive material then Asym [G 0 ] does not represent radiation. Implied meaning can be restored by appropriate alterations; but, as this point will be treated in an upcoming work, for the moment we will simply accept it as a limitation for our study.
Setting such possibilities aside, the equivalence of (7) with radiative emission is also supported both by the analogy between its operator form and that of the scattered power, and by a direct calculation for thermal (randomly fluctuating) currents [104]. By the fluctuation-dissipation theorem |J i J i | th = 4k (ω, T ) Asym [V ]/(π Z ), and so The final line above is precisely what we have derived in Ref. [68] from the perspective of integrate absorption.

B. Scattering constraints
As supported by the previous subsection, any quantity in electromagnetics can be described by combinations of G 0 , V , T , and projection operators. The basis of this reality rest on the fact that a defining relation for T , supposing G 0 and V are known, is abstractly equivalent to complete knowledge of the electromagnetic system [103]. Thus, like Maxwell's equations, any facet of classical electromagnetics, beyond the definitions of G 0 and V , must be derivable from the definition of the T operator [68,103] (The s subscript explicitly marks the domain and codomain of each operator as begin restricted to the scattering object, as opposed to the background.) To derive constraints, we will focus on the equivalent relation As both T ss and its Hermitian conjugate T † ss have support only within the volume of the scatterer, in this form, the domain and codomains of V †−1 and G 0 † are unimportant. For any projection operator (∀n ∈ N )P = P n , and so the geometric description of the scatterer contained in V †−1 equivalent to (11), with the scattering potential V †−1 and Green's function G 0 † filling whatever volume we would like to consider. Taking so that Asym [U ] is positive definite, treating the symmetric (Hermitian) and anti-symmetric (skew Hermitian) parts of (12) separately then gives Asym The constraints used to generate the cross-section bounds shown in Sec. V follow directly from (14) and (15) under the relaxation described in the next section. As recently noted in Refs. [68,69,75,79], (14) and (15) contain a surprising amount of physics. Taken together, these relations give a full algebraic characterization of power conservation [75,105], with (14) representing the conservation of reactive power and (15) the conservation of real power. (The T † ss Sym [U ]T ss piece of (14) produces the difference of the magnetic and electric energies for any incident electric field [106].) Because both real and imaginary response are captured, when these two constraints are employed there are requirements that must be satisfied on both the magnitude and phase of any potential resonance.

C. Relaxations and optimization
For the single source problems of concern to this paper, it is simplest to work from the perspective of the image field resulting from the action of T ss on a given source |S (1) , T ss |S (1) → |T . A bound in this setting amounts to a global maximization of one of the six power transfer objectives, (3)-(8), taking |T and a known linear functional S (2) | as arguments, subject to satisfaction of (14) and (15) as applied to the source and its image. So long as the known fields are not altered at previously included locations by expanding the domain, this procedure leads to domain monotonic growth: if |S (1) and |T satisfy all constraints on some subdomain, then these same vectors will also satisfy the constraints if they are embedded (included without alteration) into a larger domain. Because the value of any power objective is similarly unaffected by inclusion, the global maximum of a larger domain will thus always be larger than the global maximum of a smaller domain.
The above view also underlies the central relaxation, persisting throughout the remainder of the paper, that makes global optimization over all structuring alternatives possible. For any true T operator, nonzero polarization currents can exist only at spatial points lying within the scattering object. This fact will never be truly enforced on the image of the source resulting from the action of T (|T ), alleviating the need for an explicit geometric description of the scatterer. Rather, |T will be considered simply as an unknown vector field confined to some predefined design domain. (For instance, in the examples of Secs. IV and V the encompassing design domain is taken to be ball of radius R.) Therefore, when a bound is found, it must necessarily apply to any possible structure that can be contained in the given region, as the freedom of choosing different device geometries is already explored by optimizing over the |T field. As illustration, through this relaxation of structural information and the monotonicity property, a bound for a cuboid is both a bound for any device that could fit inside the cuboid, no matter how exotic, and for any subdesign region that could be included inside the cuboid, be it a needle, bounded fractal, or a disconnected collection of Mie scatterers.
With this easing of true physical requirements in mind, scattering operator bounds for any of (3)-(8) are equated with an optimization problem on |T : The corresponding Lagrangian is given by To help decompose later analysis, in (17) and (18) the subscript has been introduced to stand for a family of j indices coupled together by U (U = V †−1 − G 0 † ) in some complete basis {|G , j } for vector fields in domain, as would likely occur in any numerical solution approach. For the spherically bounded examples examined in Secs. IV and V is the familiar angular the momentum number, for a periodic films is an in-plane wave vector. (Generally, can be thought of as the th block, invariant subspace, in the representation G ,i |U |G , j , and j is used as a subindex.) Due to the further equivalence of the angular index and the radiation modes of a spherical domain, as well as relations with existing literature [68,101,107,108], we will also occasionally refer to families as channels (as a shorthand for radiation channels). The constraints C ζ and C γ are determined by applying (14) and (15) to { S (1) |, |S (1) }, and forgetting any information related to the geometry of the scatterer. O is either Asym [G 0 ], Asym [V †−1 ], or 0, depending on whether the problem is absorption/material loss, scattering/radiation or extracted power from a field. As exemplified in Sec. IV and illustrated in Sec. V, the necessity of conserving reactive power imparted by the symmetric γ constraint is crucial for accurately anticipating how a particular choice of material and domain influences whether or not a family can achieve resonant response.
For all cases except extracted and radiated power from an external unpolarizable source, S (2) | = S (1) |. In these instances, although (6) and (7) show that the power quantities of interest can be cast in a form similar to the corresponding initial flux problems, the inclusion of the second source image is necessary. If an unpolarizable source is taken to lie outside the domain being optimized, V −1 is defined only as a limit (tending to infinity). Once this limit is taken, the G 0 based expressions (final forms) for the extracted and radiated power result, which include the introduction of the field |S (2) = G 0 * de |J i to the objective due to the Asym [G 0 T G 0 ] term. (With e denoting the external space of the emitter and d the optimization domain.) These differences amount to the introduction of cross terms describing the interference of the fields generated by the bare and induced currents that are no longer inherently accounted for by the scattered currents at the location of the source (multiple scattering and back action). Still, the form of these problems remains like (18) up to the addition of an unalterable background contribution of Tr [S J Asym [G 0 ]].

D. Solution via duality
To solve (17), we make use of the following lemma, commonly referred to as Lagrange duality [38]. (Lagrange duality is closely associated with the alternating direction method of multipliers [109][110][111] often used for solving multiply constrained optimization problems.) Lagrange duality. Take O, (∀ j ∈ J ) E j and (∀k ∈ K ) I k to be differentiable real valued functions on subsets of R n defining a well-posed optimization problem Let m * be the corresponding maximum value and let D be the domain on which all functions are simultaneously defined, which includes the subset D * on which all constraints are satisfied.
Additionally, the function max D L is convex, and, if a set Proof. ∀x * ∈ D * , points satisfying the constraints of the original (primal) optimization, (∀k ∈ K )I k (x * ) 0, and max D L is convex as it is a sum of compositions of convex functions, max and affine functions of Whenever the operator appearing between a field and its associated linear functional in a sesquilinear relation is positive definite, the constraint describes a compact manifold. This is always true of (15), and so, as both constraints are closed sets, the domain of (17) is compact. Moreover, by the validity of the |T = 0 solution, the domain is nonempty. As such, (17) is assured to have a unique (non trivial) maximum value occurring at some stationary point (or points), and it is meaningful to consider the Lagrangian dual where the domain F is set by the criterion that max L is finite. Under this assumption, taking partial derivatives over |T , a stationary point of L requires (∀ ), is inherently both defined and positive definite. (If any A −1 is not positive definite, there is then a field |T such that G → ∞.) Letting |S (3) = (ζ − iγ )|S (1) + |S (2) , and |S (4) = A |S (3) , it follows that the unique stationary point of the dual occurs when |T = i 2 |S (4) . Hence, within F (using the above replacement), The gradients of (22) exactly reproduce the constraint equations ∂G ∂ζ = C ζ and ∂G ∂γ = C γ , with (4) . (23) Therefore, if a stationary point within the feasibility region is found, strong duality holds. The function is convex, and thus has a single stationary point, either on the boundary or within the domain of feasibility. If the point is within the domain of feasibility, then by (23) the corresponding field satisfies the constraints, and by the lemma, the corresponding point is strongly dual. In this case, the solution of (17) is with {ζ , γ } set by the simultaneous zero point of (23). If no such point exists in F, then the unique minimum value of G in the domain of feasibility, attained on the boundary of some A becoming semidefinite, remains a bound on O in (17). That is, in all cases. Comments on solving (23) are given in Sec. IV.

III. RELATIONS TO PRIOR ART
Previous work in the area of electromagnetic performance limits can be loosely classified into three overarching strategies: modal limits, shape-independent conservations limits, and scattering operator limits. Each approach presents its own relative strengths and weaknesses. Below, we provide a rough sketch of these contributions as they relate to this work, particularly the use and interpretation of constraints (14) and (15).
Arguments for response bounds based on modal decompositions, exploiting quasinormal, spectral, characteristic, Fourier and/or multipole expansions [85,[112][113][114][115][116][117][118][119][120][121][122][123][124][125], have been widely considered for many decades. And like the classical diffraction and blackbody limits of ray optics, they have proven to be of great practical value for describing possible interactions between large objects and propagating waves [60,126,127]. Yet, in complement, the need to enumerate and characterize what modes may possibly participate has also long proved problematic. Small sources, separations, and domains typically require many elements to be properly represented in any basis well suited to analysis of Maxwell's equations, and so, especially in the near-field and without knowledge of the geometric characteristics of the scattering object, or inclusion of material considerations, there is no systematically effective approach to bound modal sums (without introducing additional aspects as is done in scattering operator approaches). While a variety of considerations (transparency, size, etc.) have been heuristically employed in an attempt to introduce reasonable cutoffs [58,67,107,[128][129][130], the values obtained by modal methods in such settings are consistently orders of magnitude too large [71,79,131]. Nevertheless, the idea that modes often separate otherwise muddled aspects of photonics remains a key insight.
Shape-independent conservation limits, utilizing energy [59,71,78] and/or spectral sum rules [99,101,[132][133][134][135][136][137] to set local limits based on physical laws, generally display the opposite behavior, and are known to give highly accurate estimates of maximal far-field scattering interactions in the limit of vanishingly small sizes (quasistatics) for certain (near resonant) metallic materials [9,10,71,77,136]. Notwithstanding, as we have found in our work on bounds for radiative heat transfer [69,79] and angle-integrated radiative emission [68], they are not sufficient in and of themselves to properly capture various relevant, performance limiting, wave effects. Intuitively, without any geometric information, a conservation argument must apply on a per volume basis, which is at odds with the area scaling of ray optics.

033172-8
As a relevant example, consider the fact that the global power quantities given in Sec. II must be non-negative. Two of these turn out to be unique and thus set physically motivated algebraic constraints on the T operator. For any vector field |E , the non-negativity of scattering (known as passivity [138]) imposes that while the non-negativity of absorption imposes that Both conditions are concurrently captured in (15), which is equivalent to the optical theorem [105]: physically, the sum of the absorbed power and scattered power, (4) and (5) If no additional information pertaining to possible geometries or the characteristics of |E is given, then the most that can be said from (15) is that no singular value of T can be larger than the inverse of the smallest singular value of Asym [V −1 † ]. For a single material with response characterized by a local electric susceptibility χ , this logic yields the material loss figure of merit originally derived in Ref. [71] directly using the implications of passivity for polarization fields. The universal applicability of this largest possible response has profound consequences for the design of many photonic devices relying on weakly metallic response (−1 Re [χ ] −10) and small interaction volumes [67,77]. In these cases, it is often fair to assume that a resonance can be created and that Asym [G 0 ] ≈ 0. Little structuring is required to achieve a plasmonic resonance, and the maximum achievable polarization current is indeed dominated by material losses. But, for single material devices where light-matter interactions occur on length scales comparable to (or greater than) the wavelength ( λ/10), or the real part of χ is outside the range stated above (e.g., strong metals or dielectrics), such estimates are overly optimistic for single material devices, Sec. V. Over a large enough domain, the generation of polarization currents capable of interacting with propagating fields leads to radiative losses which have been neglected by supposing Asym [G 0 ] ≈ 0. To create a bright (active) far-field resonance, it must be possible to couple to radiation modes and then confine the resulting generated field within the domain. This is not always possible for any prescribed choice of material and device size. Scattering operator approaches aim to eliminate the weaknesses of modal and shape-independent conservation arguments by combining their strengths [68,69,102,107,[139][140][141][142][143][144][145]. Innately, the Green function of an encompassing domain (through its link to Maxwell's equations) provides both a modal basis for, and constraints on, modal sums. In concert, restrictions on the possible characteristics of the T operator can be used to ensure that physical laws and scaling behavior are properly observed. A number of encouraging conclusions have been derived in this manner. Drawing from our own work, in Ref. [69], it was shown that imposing (25) on the operator expression for angle-integrated absorption and thermal emission, Eq. (9), is sufficient to generate bounds smoothly transitioning from the absorption cross-section limit of a resonant metallic nanoparticles (the product of the volume and ζ mat ) to the macroscopic blackbody limits of ray optics. Similar methods were used in Ref. [79] to prove that, for equal values of ζ mat , nanostructuring cannot appreciably improve near-field thermal radiative heat transfer compared to a resonant planar system.
Nevertheless, careful investigation of the previous situations where scattering operator amalgamations have been successfully applied reveals a consistent use of niceness properties that are not generally valid. In the examples given above, we were aided by the fact that thermal sources are completely uncorrelated and, for thermal emission and integrated absorption, that only propagating fields needed to be treated. Without these helpful aspects, there are situations where past scattering operator approaches, which focused exclusively on conservation of real power (15), would add complexity without tightening the asymptotics provided by shape-independent conservation arguments (dashed lines in the figures of Sec. V). Furthermore, using the currently practiced technique of translating established physical principles back to implied operator properties and then using inequality compositions to produce limits, it is difficult to see how the interaction of more than one or two additional constraints could ever possibly be accounted for. The switch to exploiting algebraic deductions beginning from (10) in combination with standard optimization theory may seem a subtle distinction; however, the flexibility offered by Lagrange duality suggests that this view may be of substantial benefit going forward.
A related shift towards systemization has been realized in the recent report of computational bounds by Angeris, Vučković, and Boyd [39], treating linear electromagnetics as an optimization problem with respect to a target field (or a collection of target fields). The result, also making of use Lagrange duality, has immediate consequences for qualitatively understanding and improving inverse design. Yet, it does not allow one to make conclusive statements about feasibility and relative merits, as is true of traditional limits. More precisely, what is found is a "computational certificate": given a target field and an evaluation metric, the algorithm returns a number; any vector satisfying Maxwell's equation will have a metric disagreement with the target at least as large as the number. That is, the algorithm does not find physical limits, but instead a minimum bound on distance, in a certain user determined measure, between a particular field and the set of physically possible fields. There may be situations where this difference is of little consequence, or provably zero, but a priori there are no guarantees. There need not be any relation between the value taken by a function at a point and how near that point is to some set.
Finally, while concluding the writeup and peer review of this paper, we have become aware of contemporary works by Gustafsson et al. [75], extending developed methods for bounding the performance of radio frequency antennas [9,10,70,146], Kuang et al. [72], and Trivedi et al. [76]. Independently and simultaneously developed, the first two formulations are in many respects quite like the method presented here. Working from the perspective of polarization currents, both articles bound the optimization of objectives 033172-9 equivalent to (4)-(8) subject to power constraints via Lagrangian duality. Ref. [75] incorporates both (14) and (15) while Ref. [72] uses only (15), leading to the discrepancies between the dashed and solid lines depicted in Sec. V. Taken together the two articles present a rich collection of findings reinforcing some of the observation of Sec. V. Nonetheless, tangible differences do exist. Although largely unexploited in this initial investigation, there are many advantages offered by working with operator as opposed to vector relations for further generalizations. In contrast, the more recent work of Trivedi et al. [76] focuses on how the physically required self-consistency of scattering theory sets constraint on the possible properties of the scattered electric field. Branching off from this alternative starting point, the optimization problem and subsequent application of Lagrange duality acquires characteristics quite different from those encountered here. The closest comparison would appear to be something akin to imposing only the conservation of reactive power, as the approach described in Ref. [76] appears to produce nontrivial bounds only when resonant response is not possible.

IV. COMPUTATIONAL MECHANICS AND SINGLE CHANNEL ASYMPTOTICS
To elucidate the mechanics of (17), we now describe the basics of our computational procedure for the specific example of a spherical confining boundary. The discussion is broken into two subsections. The first sketches an outline of the approach by which the results of Sec. V for compact bodies were obtained. The second considers a simplified singlechannel (family) version of (17) that becomes exact as R → 0, predicting many of the trends seen in Figs. 2 and 3 of Sec. V. Namely, the largest possible interaction enhancements are found to obey either an effective medium "dilution" response, or the material dependence encountered in Rayleigh scattering [84]. For low-loss dielectrics (Re [χ ] > 0) and strong metals (Re [χ ] −1), this can lead to large discrepancies with respect to previously established per volume bounds based on the material loss figure of merit ζ mat [67][68][69]71,77,78].

A. Computational mechanics
Recall that, once an origin has been specified, the Green function can always be expanded in terms of the regular (finite), RN and RM, and outgoing, N and M, spherical wave solutions to Maxwell's equations as [103,147] In (26), x and y are used to denote the wave vector normalized radial vectors of the domain and codomain, i.e., x = 2π r/λ, θ, φ , with x and y used for the corresponding radial parts. The integral over Y is taken to mean integration over the y coordinate. [Note that there is no complex conjugation in these integrals, and that our notation for the Green function is unconventional in that an additional factor of k 2 = (2π/λ) 2 is included as part of the definition, as opposed to a separate conversion factor, allowing all spatial distances to be normalized in terms of the wavelength.] So long as the current source is not located within the domain in question, any incident field can be expanded in term of the regular waves [103,147,148]. As such, the spectral basis of the asymmetric part of (26) the unit normalizedR M ,m andRN ,m waves, serves as convenient choice for starting point for generating the basis vector families required to properly represent U [149]. That is, given the form of the regular solutions the orthonormality of the vector spherical harmonics (A (1) ,m , A (2) ,m , and A (3) ,m , see Ref. [150] for details) means that the Green function (26) does not couple the , m or RN and RM labels. Hence, the individual radiation channels act as an effective partitioning, and by taking these vectors as the "seeds" or "family heads," a complete (simplifying) basis for (17) can be generated through the Arnoldi (Krylov subspace) procedure [151]. Briefly, starting with a given unit normalized regular wave,RN ,m , one generates U |RN ,m = (V †−1 − G 0 † )|RN ,m . Projecting out theRN ,m component of this image and normalizing, one obtains a new vector |PN (2) ,m . |PN (2) ,m then serves as the input for the next iteration, and in this way the family (block), more properly thê RN ,m family, of the matrix representation of the U operator (U ) is computed, i.e., U is represented in the basis {|RN ,m , |PN (2) ,m , |PN ,m , . . .}. Technically the above process does not terminate, but regardless, two practical consideration lead to workable numerical characteristics [152]. First, due to the fact that each vector is orthogonal to all others, the off-diagonal coupling components of U in every family originate entirely due to the volume integrals in (26). Therefore, in the limit of vanishing volume, or high , each U is effectively 2 × 2. Second, by the Arnoldi construction, all upper diagonals beyond diag 1 , with diag 0 standing for the main, are zero. Thus, because at each step the generation of new basis components is driven entirely by Sym [G 0 ], and Sym [G 0 ] = Sym [G 0 ] † , the matrix representation of every U is tridiagonal. Due to this banded nature, each A −1 gives a simple, conclusive, estimate of the error for images generated by A . All that is required is to pad the current solution with zeros and calculate its image under A −1 in a basis augmented by three additional elements. The magnitude of the error of the image compared to the source is exactly the same as would be found in any larger (even infinite) basis, see Appendix for additional details. Coupled with (22), the determination of bounds for any electromagnetic process in a compact body (an object of bounded extent) that can be described as a total absorption, scattering or extinction process, is thus mapped to the 033172-10 numerical determination of the minima of a constrained convex function. Many efficient algorithms exist to solve such problems [153][154][155], along with a variety of nice introductions [38,151,156].

B. Single channel asymptotics
Consider the simplified optimization problem max T (1) P 1 T (1) such that where it has been assumed that that (∀ f i = ) S f i | = 0, and P 1 represents the projection of |T onto the = 1 "family head." So stated, (29) represents the maximum possible interaction that can occur between a generated polarization current and an exciting field for a single radiation channel, respecting the conservation of total power. Based on the Taylor series representation of the spherical Bessel functions, there are two situations in which this problem is fairly simple to treat analytically. If either the radius kR 1, with k = 2π/λ, or is large compared to kR, then both the regular and outgoing waves appearing in the Green's function, (26), are well approximated by two-term expansions. This feature causes the above Arnoldi procedure for basis generation to effectively terminate after constructing a single image vector. As high contributions (for propagating waves) occur only when large low contributions are also present, no meaningful asymptotic behavior can be extracted from a single channel analysis of the second possibility. Accordingly, we will focus on the assumption that kR 1. Symbolically carrying out the required image generation and orthogonalization steps, the representation of U f 1 in this quasistatic (R → 0) regime is Denoting the symmetric and anti symmetric components of the representation of U 1 as u (·,·) s and u (·,·) a , solving (17) amounts to determining the {t 1 , t 2 } component pair producing the largest magnitude t 1 such that Here, the t 1 and t 2 variables are the (positive) magnitude coefficients of |T 1 in the first and second Arnoldi vectors of the = 1 family, s 1 is the coefficient of the source, θ is the relative phase difference between the source and first coefficient of |T 1 , and φ is the relative phase difference within the two coefficients of |T 1 . As a response operator, Asym [T ] must be positive semidefinite and so θ ∈ [0, π] [157]. Using the symmetric constraint to solve for t 2 in terms of t 1 , forgetting off-diagonal terms when they appear as sums against diagonal terms in the resulting quadratic equation, the asymmetric constraint determines that the maximum polarization current that can couple to the source is subject to the condition, resulting from the requirement that t 1 and t 2 are real, that Neglecting all higher-order corrections but the radiative efficacy ρ GN 1 = Asym [G 0 ] (1,1) , the relative magnitude t 1 /s 1 = ζ eff is therefore limited by where δ GN 1 is the domain dependent, material independent, portion of u (1,1) s − u (2,2) s [158,159]. In (33), the selected symbols are motivated by three factors: (1) analogs to the analysis given above likely exist whenever all dimension of the design domain are small compared to the wavelength, (2) δ GN 1 is the delta function portion of (26) for the regular RN 1,m wave of a spherically bounded domain, δ GN = 1/3, and (3) we have previously used ρ GN 1 = 2 9 ( 2πR λ ) 3 to denote "radiative efficacy" in studying the existence of bounds on radiative thermal power transfer [69,79]. Utilizing |T = ζ eff |S 1 in the power forms provided in Sec. II produces the single channel cross-section limits seen in Sec. V, (38) and (40). Full solutions of (17) for an incident planewave are found to be accurately predicted by (33) as R → 0, Figs. 2 and 3, outside the −3 Re [χ ] −1 region where the assumption that the u (1,2) s terms can be neglected does not hold. While the monotonicity property of these bounds means that (properly scaled) they are accurate for any compact domain geometry, and one may reasonably guess that the characteristics of the Arnoldi process on which the above arguments rest are similar in any small volume limit, it should be kept in mind that other domain geometries (e.g., ellipsoids [71,136]) may well display stronger per volume response, if the resonance condition shifts to more negative values of Re [χ ] (δ GN 1 > ⅓). Still, while the generated polarization current and exciting plane wave may exhibit a larger average interaction within the volume of a structured body, the net enhancement will be weaker than ζ eff once the ratio of its volume to an encompassing ball is accounted for.
Given the revealing nature of this simple symbolic analysis for spherically bounded domains, it is reasonable to ask why we have not included corresponding results for periodic films. Briefly, the reason for this omission has to do with the altered proprieties of G 0 that occur in the presence of two periodic boundary conditions. Unlike the case of the compact domain, in the limit of vanishing film thickness t → 0, important δ function contributions of G 0 are simultaneously spread across many in-plane wave vectors, with no one channel containing all key characteristics. The only case where single channel asymptotics reproduce observed behavior is for extremely thin (nonresonant) dielectric films.

V. APPLICATIONS
In this section, the program developed above is exemplified for two canonical scattering processes: limits on absorbed and scattered power for a plane wave incident on any structure that could be confined within a ball of radius R, characterizing the maximum cross-section enhancement a body may exhibit [103] subject to the constraints that net real and reactive power are conserved; and limits on the absorptivity of a periodic film for a normally incident plane wave as function of the film thickness, t. Because of the complete exploration of structural possibilities conducted in calculating these bounds, and the domain monotonicity property explained in Sec. II, all results are equally applicable to any subdomain, or disconnected collection of subdomains, that fit inside any given ball or periodic film; therefore, in addition to any possible connected structures, the bounds pertain to geometries like arrays of Mie scatterers [145,160] or plasmonic building blocks [161,162]. Throughout the section, R and t are unnormalized unless otherwise stated, and stands for the angular momentum number (spherical harmonics A (−) ), originating from representing (17) using the Arnoldi construction described in Sec. IV.
For both spherically bounded examples, the distribution of the power density within the domain between different angular momentum numbers is strongly tethered to the radius of the boundary. Specifically, the coefficients of the electric field, for a unit normalized electric field amplitude, in terms of the regular (finite at the origin) RM ,m and RN ,m solutions of Maxwell's equations in spherical coordinates, (28), arê with r standing for the wave vector normalized radius (the product of the true radius and k = 2π/λ). Through Asym [U ] (34) establishes a link between the magnitude of the radiative efficacy of each channel [69], or number, and its potential for enhancing scattering cross sections: if the plane-wave expansion coefficient of a given channel is large, then so is Asym [G 0 ], tightening the constraint that real power must be conserved. This leads to a complementary action of the power constraints. For any given combination of material and radius, save Re [χ ] = −3 in the R → 0 limit, either real or reactive power conservation limits induced polarization currents in the medium more severely than what would be expected based solely on the material loss figure of merit widely considered in past work on electromagnetic bounds for arbitrary materials and structures [67,69,71,77]. (An explanation of the origin and usefulness of this quantity is given in Sec. III.) In Figs. 2 and 3, dashed curves depict cross-section limits attained when only the conservation of real power, (15), is imposed, as in Ref. [72], and solid curves result when reactive power conservation is additionally included. All results are found using the Lagrange duality approach described in Sec. II and are strongly dual [38]. Conversely, the example bounds for periodic films shown in Fig. 4 are strongly dual only when absorption is near zero. Outside this regime, and the associated sharp transition to resonant response, all major features, including the half absorption plateaus, are seen to be accounted for by the conservation of real power, and are thus described by the growth of the radiative efficacies (singular values of Asym [G 0 ]) of the two channels with zero in-plane wave vector [68,72].

A. Quasistatic regime (R/λ → 0, t/λ → 0)
Recalling the conclusions of Sec. IV B, the simultaneous conservation of real and reactive power has far-reaching implications for electromagnetic power transfer when all dimensions of the confining domain are small. The analog of the optical theorem for reactive power, (14), adds phase information on top of the maximum polarization magnitude set by the conservation of real power, (15); and so, when both constraints are taken into account, (17) captures the fact it is not always possible to produce a resonant structure given any single material, of electric susceptibility χ , and a maximal characteristic size, R. In fact, there are rather strict requirements that must be met for resonant response to be achievable. Crucially, it must be possible to effectively confine the scattered electromagnetic field, resulting from the polarization currents created in the structure by the incident (source) wave, within the spherical volume. As validated by the cross-section bounds of Figs. 2 and 3, the only mechanism by which such confinement can be achieved while also allowing interactions with propagating waves as R/λ → 0, is the excitation of localized plasmon-polaritons, occurring for Re [χ ] = −3 if the domain is completely filled with material [80].
If Re [χ ] is larger than this value, excluding small deviations that appear for weak metals between −2 Re [χ ] −3, then, as confirmed by the tiny achievable cross-section values, resonant response with a propagating wave is not possible. Therefore the largest allowed power transfer happens when the material simply fills the entire domain, and as such, maximal scattering cross sections exhibit the same susceptibility dependence encountered in Rayleigh scattering [103]. Applying (33), the maximal magnitude of the interaction that can occur in a dielectric structure between the (normalized) 033172-12 FIG. 4. Absorption cross-section bounds for a plane wave incident on a periodic film. Absorption cross-section bounds for any periodic structure (e.g., grating, metasurface) of period and electric susceptibility (a) χ = 13.6 + i 5 10 −2 (silicon at 0.8μm) and (b) χ = −87 + i 10 (gold at 1.5μm) that can be contained within a planar boundary of thickness t/λ. Dots correspond to actual absorption cross-section values discovered using inverse design following the algorithm described in Ref. [163]. A number of visualizations of the associated structures, along with a schematic of the absorption process investigated, are included as insets. All discovered designs with t/λ < 10 −2 consist of a single layer, uniform in the normal direction. For larger thickness, optimization of three equally thick uniform layers is needed to achieve the reported absorption cross sections (eventually reaching absorption saturation). Like prior figures, solid lines are found by solving (17) under the constraint that total power is conserved, while dashed lines, coinciding with the solid line in (b), are found by only imposing the constraint of real power conservation. As shown in (a), total power conservation bounds for dielectric media depend on the square period supposed. Moving from left to right, the periods used in the inverse design optimizations are = {5λ, 5λ, 5λ, 5λ, 4λ, 2λ, 3λ/2, 3λ/2, λ, λ}. Results for metals (Re [χ ] < 0) do not depend on the periodicity of the system. Discussion of additional bound features relating solely to imposing the conservation of real can be found in Refs. [68,72]. incident field and the polarization current it excites is with denoting the radiative efficacy of the = 1, m = ±1, 0N polarized channel. Employing the power forms given in Sec. II, the scattering cross section of any structure, under the above assumptions, must obey the relations In stark contrast to ζ mat , ζ Ray decreases for increasing Re [χ ] and has a negligible dependence on material absorption, Im [χ ]. Comparing the dashed and full lines of Figs. 2 and 3, particularly Figs. 2(b) and 3(b), the resonance gap between these two forms can be quite extreme for realistic dielectrics. For example, ζ mat ≈ 10 7 for silicon at λ = 1.5 μm, with χ ≈ 11 + i 10 −5 [164].
As Re [χ ] shifts to increasingly negative values, geometries supporting localized plasmon-polariton resonances become possible, and past Re [χ ] = −3 cross-section limits display resonant response characteristics. The power exchange between the incident field and the generated polarization currents is then, asymptotically, restricted to be smaller than the "diluted" material figure of merit leading to cross-section limits of The "dilution" modifier is chosen as the form of (39), disregarding the radiative efficacy ρ GN 1 , is equivalent to the material loss figure of merit ζ mat if a "dilution factor" is introduced to shift Re [χ ] to −3. That is, considering the low-loss limit Im [χ ] |χ |, if it is supposed that the magnitude of Re [χ ] is rescaled to match the localized resonance condition of a spherical nanoparticle, then ζ mat → 3|χ |/ Im [χ ], which is equal to the expression for ζ dil in the limit ρ GN 1 → 0. (Due to its connection with the localized plasmon resonance of a spherical nanoparticle, the ratio |χ |/ Im [χ ] is commonly 033172-13 encountered in discussing the potential of different material options for plasmonic applications [165,166].) The validity of (39) internally rests on the assumption that the wavelength is much larger than any structural feature. Since this is also the central criterion for most homogenization descriptions of electromagnetic response to be applicable [167], it is sensible that material structuring limited to tiny domains can, at best, alter effective medium parameters [81][82][83]. Equation (39) proves that the implications of this picture are universally valid for both scattering and absorption in strong metals. However, many commonly stated effective medium models also predict that there are structures capable of creating effective susceptibility responses more negative than either of the constituent materials [168][169][170][171][172]. For example, the Maxwell-Garnett formula for monodispersed spherical vacuum inclusions in a background host is where f is the volume filling fraction of the inclusions and χ h is electric susceptibility of the host [86]. Based on (41), using the iterative argument given in Ref. [86], it should be anticipated that low-loss resonant response would be achievable shortly after Re [χ ] drops below −1. While cross sections do begin to grow before Re [χ ] = −3, it is clear that this Maxwell-Garnett condition is not sufficient. Dilution via (41) also yields a different dependence on material loss-values than that predicted by (39), with (41) consistently giving slightly larger effective losses. It is also interesting to compare (38) and (40) with coupled mode descriptions of scattering phenomena in the single channel limit [85], where γ rad and γ abs are the geometry-specific radiative and absorptive decay rates associated with a given resonant mode of frequency ω res . Up to a missing factor of 4, which is accounted for by the facts that (38) and (40) represent maximum quantities [116,142] and that scattering cannot occur without absorption [71], there is a clear symmetry of form between (38) and (40) and (42) when the system is off resonance, and when the system is on resonance; in both situations, γ rad → ρ GN 1 . Since (38) and (40) are bounds, and not descriptions of any particular mode, these associations may be understood as "best case" parameters for what could be achieved in any geometry supporting a single mode, and are thus closely linked to prior limits based on coupled mode theory [85,117,119,126,129,130]. Notably, the comparison precludes any resonant geometry from achieving the rate-matching condition of γ rad = γ abs if it is confined to a small ball. Precisely, the only candidate materials are fictitious metals with −3 Re [χ ] and Im [χ ] → 0, since the radiative efficacy ρ GN 1 → 0 with vanishing object size.
This situation, a fictitiously low loss metallic nanoparticle, is also the most relevant condition under which the bounds asymptotically reach arbitrarily large values. However, as we have discussed in Ref. [68] in the context of angleintegrated planewave absorption, unbounded growth requires saturation of an unbounded number of angular momentum indices (radiation channels). For any particular , saturation is approximately achieved as R → 0 when ρ GN , ρ GM Im [χ ]/|3 Re [χ ]|. Therefore the relation between the radiative efficacies and the angular momentum number , with imparted through the cylindrical Bessel functions, with , implies that so long as real power is conserved, bounds on cross-section enhancement exhibit sublogarithmic growth with vanishing material loss, Im [χ ] → 0. (Proof of this statement, in all important regards, follows from the derivation given in Ref. [68]. The stated inequality follows from the power series of the cylindrical Bessel functions [173].) As seen in the supporting figures, Figs. 2 and 3, in practice this scaling behavior is of little consequence.
For periodic films (e.g., gratings, photonic crystals, and metasurfaces), the central feature of (17) absent from the models of Refs. [68,72] is the initial suppression and sharp onset of resonant absorption for thin dielectrics. The above findings for compact domains, and practical experience, both suggest that the existence of such a critical thickness (depending of the magnitude of Re [χ ]) for dielectric materials is reasonable. However, the dependence of this thickness threshold on the period of the system, physically associated with the presence of leaky-mode resonances [174], is perhaps less expected. The origin of the relation traces to the properties of G 0 for an extended (infinite) system. Crossing over the light line boundary between propagating and evanescent waves, there are vectors within the basis described in Appendix, |X = |RX (−) (k p ) , that allow X| Sym [G 0 ]|X to be arbitrarily negative and X| Asym [G 0 ]|X = 0. These characteristics allow reactive power conservation, (14), to be trivially satisfied for any possible polarization current. However, when a finite period is imposed on the system, modes arbitrarily near the light line are not allowed, and in turn, the necessity of conserving reactive power may imply that resonant response is not possible  for particular values of t and χ . As exemplified in Fig. 4, knowledge of the critical thickness at which such leaky modes can be supported for a given period and material may be of substantial benefit to the design of large scale metasurfaces [163,175,176].
For boundary radii approaching wavelength size, the applicability of the quasistatic results for spherical confining region quoted under the previous subheading becomes increasingly tenuous. The growth of planewave amplitude coefficients into angular momentum numbers (channels) beyond = 1 opens the possibility of utilizing a wider range of wave physics (e.g., leaky and guided resonances [177,178]), and correspondingly, reactive power conservation (resonance creation) becomes a weaker requirement.
These factors lead to a more intricate interplay between the two power constraints, causing the sharp jumps observed for dielectrics in Figs. 2 and 3, which manifest, mechanically, in rapid changes to the properties of the scattering T operator constraint relations, especially (14) applied to dielectric materials. The behavior is first observed in the = 1 channel, with the initial peaks in Figs. 2(b) and 3(b) following closely after the half wavelength condition 2π r λ ∂r = 0, and the second peaks occurring near the full wavelength condition, This second criterion is also the approximate resonance location for a homogeneous dielectric sphere of index √ χ [179], the spherical analogs of the Fabry-Perot condition [180], making its appearance consistent with the Rayleigh response predictions of (36). At the same time, as previously remarked, the inflation of the boundary also increases the radiative efficacy of each channel as described by (43)  For periodic films, all features seen at both wavelength and ray optic thickness scales are fully accounted for by the conservation of real power and the associated radiative efficacies of Asym [G 0 ] at normal incidence. Detailed discussions of these quantities are given in Ref. [68] and Ref. [72]. The most interesting result, that the bounds plateau for film thickness between roughly t/λ = 10 −4 and 10 −2 , is caused by the fact that at these thicknesses only a single symmetric mode is bright. Accordingly, the power radiated to the far-field by the excited polarization current can not be canceled, and only half of the total power of the incident wave can be extracted through absorption [68].

C. Ray optics regime (R/λ 1, t/λ 10 −1 )
In the large boundary limit, achievable scattering interactions in any given channel are increasingly dominated by the conservation of real power through the growth of radiative losses. Correspondingly, the dash bounds, calculated by asserting only that the sum of the scattered and absorbed power must not exceed the power drawn from the incident beam, coincide with those arising from total power conservation to increasingly good accuracy. Making this reduction, limits for either cross-section enhancement quantity become largely congruous to the angle-integrated absorption bounds given in Ref. [68]. The planewave expansion coefficients of (34) exhibit exactly the same per-channel characteristics considered in that paper, and so, the same asymptotic behavior is encountered. Regardless of the selected susceptibility χ , for a sufficiently large radius, each of the power objectives described in Sec. II begins to scale as the geometric cross section of the bounding sphere. For absorption, this leads to a value equal to the geometric cross section of the confining ball, π R 2 . For extinction and scattering, a value of 4π R 2 is found, two times larger than what would be expected based on the extinction paradox [181,182]. The genesis of this additional factor is presently unknown, and investigation of the properties of the optimal polarization current of these curious results merits further study.

VI. SUMMARY REMARKS
The ability of metals and polaritonic materials to confine light in subwavelength volumes without the need for any other surrounding structure (plasmon-polaritons [183,184]), coupled with the variety of geometric wave effects achievable in dielectric media (band gaps [185,186], index guiding [187,188], topological states [189,190]), rest as the bedrock of contemporary photonic design. Yet, the relative abilities of these two overarching approaches for controlling light-matter interactions remains a widely studied topic [191][192][193]. The broad strokes are well established. The possibility of subwavelength confinement and large field enhancements offered by metals is offset by the fact these effects are fundamentally linked to substantial material loss [192]. Through interference, dielectric architectures may also confine and intensify electromagnetic fields, and can do so without large accompanying material absorption [51]; but, accessing this potential invariably requires larger domains and more complex structures. While comparisons within rigidly defined subclasses have been made [64], the merit of a particular method for a particular design challenge is almost always an open question. As with the rising need for limits in computational approaches highlighted in the introduction, a central driver of debate is the lack of concrete (pertinent) knowledge of what is possible, beyond qualitative arguments.
We believe that the simple instructive cross-section examples shown in Sec. V are compelling evidence that the generation of bounds based on constraints derived from the T operator and Lagrange duality offers a path towards progress;  and that by translating this method beyond the spectral basis employed here, onto a completely geometry agnostic numerical algorithm, it will be possible to analyze the relative trade offs associated with various kinds of optical devices. Through bound calculations varying material and domain parameters, the significance of different design elements from the perspective of device performance should be ascertainable in a number of technologically relevant areas. The basic scattering interaction quantities given in Sec. II lie at the core of engineering the radiative efficacy of quantum emitters [88,89], resonant response of cavities [90][91][92], design characteristics of metasurfaces [160,194,195], and the efficacy of light trapping [7,126] devices and luminescent [94,95] and fluorescent [96,97] sources. They are also central building blocks of quantum and nonlinear phenomena like Förster energy transfer [196], Raman scattering [78], and frequency conversion [51].
As seen in Sec. V, relations (14) and (15) are amenable to numerical evaluation under realistic photonic settings (for practical domain sizes and materials) and sufficiently broad to provide both quantitative guidance and physical insights: as the size of an object interacting with a planewave grows, there is a transition from the volumetric (or super volumetric) scaling characteristic of subwavelength objects to the geometric cross-section dependence characteristic of ray optics; critical sizes exist below which it is impossible to create dielectric resonances; material loss dictates achievable interactions strengths only once it becomes feasible to achieve resonant response and significant coupling to the incident field.
Several generalizations of the formalism should be possible. First, there is an apparent synergy with the work of Angeris, Vučković, and Boyd [39] for inverse design applications. The optimal vectors found using (17) provide intuitive target fields. Second, following the arguments given in the work of Shim et al. [99] it would seem that (17) can be further enlarged to include finite bandwidth dispersion information, accounting for the full analytic features of the electric susceptibility χ (ω). Finally, by combining the respective strengths of both classes of materials, hybrid metal-dielectric structures offer the potential of realizing more performant devices. The generalization of (17) to incorporate multiple material regions (multiregion scattering [69]) as an aid to these efforts stands as an important direction of ongoing study. As we have stated earlier, as the method relies only on scattering theory, almost all lines of reasoning we have presented apply equally to acoustics, quantum mechanics, and other wave physics.  With perfect numerical accuracy, the convergence of A is guaranteed in a finite number of iterations. The strictly diagonal elements of each U matrix, V †−1 , remain constant while the off diagonal coupling coefficients introduced by Sym [G 0 † ] gradually decay with every iterations. Thus, at a certain point, the diagonal V †−1 entries eventually overwhelm all other contributions, terminating the U matrix. (As the magnitude of the susceptibility considered increases, V †−1 shrinks and more Arnoldi iterations are required.) Still, there are pitfalls that must be avoided when numerically implementing an Arnoldi iteration, caused by the singularity of the outgoing N waves at the origin. The issue is illustrated by considering the image of RN under G 0 (26), with using the normalized vector spherical harmonics as described in Ref [197]. Near the origin, r → 0, the leading order radial dependencies of (28) and (A1) are ,m ,m From (26), the image of RN ,m under the Green function restricted to a spherical domain with radius R is where the final term is the δ-function contribution, and the RN co (r) and N co (r) terms are given by Exploiting the orthogonality of the vector spherical harmonics, the leading radial order for N co (r) is r 2 +1 . Therefore the N(r)N co (r) term has a leading radial order of r −1 , the same as the starting vector RN(r). At first sight, RN co (r) is more troubling. The dominate radial orders are r − for r 2 N(r ) and r −1 for RN(r ). Thus it would seem that the integrand has an r −1 dependence, which would result in a logarithmic divergence at the origin. A more careful consideration, however, shows that the leading order terms from A (2) ,m and A (3) ,m 033172-16 cancel as RN co (r) = i i ( + 1) 2 + 1 r −1 A (2) ,m · A (2) ,m The key to this cancellation is the ratio of the A (2) ,m and A (3) ,m terms, √ ( + 1)/ . So long as this ratio is maintained, the RN co factor does not generate logarithmic contributions, and in turn this causes the leading order ratio to remain intact under the further action of G 0 . By insuring that this does in fact occur, the Arnoldi process may continue to stably iterate until convergence is achieved. Consider any vector P = p r −1 A (2) ,m + + 1 Substituting back into (A8) then gives ,m ,m + C(R)RN(r).
Hence, as anticipated, all components retain a √ /( + 1) ratio. By induction, this argument extends to every step of the Arnoldi process, generating vectors well behaved at the origin. In implementation, care must be taken not to let numerical error push this component ratio away from √ /( + 1) at any step. (Otherwise, the logarithmic divergence will quickly destabilize new image vectors.) This precludes the use of spatial discretization based representations, since for finite grids discretization error is inevitable and leads to a rapidly growing instability. We have circumvented this issue by representing the radial dependence of the Green function and Arnoldi vectors by polynomials (Taylor series). For larger domain sizes, this approach demands a high level of numeric precision, and so, the PYTHON arbitrary precision floating-point arithmetic package MPMATH was used in all calculations [198]. When determining the image of a vector under the Green function, the tiny coefficient of r −1 due to numerical errors from the finite Taylor series and set floating-point precision were explicitly truncated (ignored). With sufficiently high precision and representation order the Arnoldi process can be performed stably and accurately up to convergence of each U matrix. Much of the difficulty, and inefficiency, associated with this method stems from working in spherical coordinate, which are inherently ill defined at the origin.

APPENDIX B: RADIATIVE (FOURIER) BASIS FOR FILMS
Working in Cartesian coordinates, consider a planar slab of thickness t with normal directionẑ. Let k p denote the in-plane wave vector with corresponding spatial vector x p . Following our notation of including an extra factor of k 2 as mentioned above, the Fourier representation of the Green function given in Ref. [103] below the light line (k p /k < 1) is given by and o = (k yx − k xŷ ). The four terms are, respectively, the upward traveling and downward traveling ordinary and extraordinary waves [199]. In (B1) and (B2), all wave vectors are normalized by the standard wave vector k = 2π/λ, and all spatial coordinates are multiplied by k; k z is defined to be the positive square root k z = 1 − k 2 p ; k (±) = k p ± k zẑ ; and when a wave appears on the right of a ⊗ symbol complex conjugation is implied. Antisymmetrizing, the skew-Hermitian component is The above basis waves are also valid for periodic films. The only change required is to replace the continuous index k p by its appropriate discrete counterpart.

APPENDIX C: INVERSE DESIGN PROCEDURE
All computational geometries were discovered following a standard topology optimization algorithm [44], where each pixel (susceptibility value) within the bounding sphere is considered as an independent design parameter, based on the method of moving asymptotes [200]. To begin, each pixel (index by x) is allowed to linearly explore a continuous spectrum of susceptibility values varying between the background χ x = 0 and the stated material χ x = χ . The resulting parameter values are then iteratively used as starting points for new optimizations subject to increasingly strict regularization filters, as described by Jensen and Sigmund [43], in order to enforce binarization (each pixel must correspond to either the material χ or background).
The primary challenge of applying this approach to maximize the scattered power from a planewave incident on a three-dimensional compact object, Eq. (5), lies in the computational cost of P sct flx , which is typically evaluated several thousand times during the course of any one optimization. Without simplification, each evaluation requires a new solution of Maxwell's equations in three-dimensional space. To avoid this otherwise significant computational challenge, the results presented in Sec. V exploit a fluctuating-volume current (FVC) formulation of electromagnetic scattering [69,201] and the efficiency of stochastic singular value decomposition. (Descriptions of the various subroutines entering into this program are given in Refs. [104,202,203].) This approach both limits the spatial solution domain to the bounding sphere, and allows the central matrix-vector multiplication, G 0 |J g , to be accelerated via the fast-Fourier transform.
More precisely, optimization is carried out using the first form of the scattered power quoted in (5), with |T = T ss |E i , as an objective. To compute the real number associated with this form for any particular structure, one inverse solve (sped up by applying a fast-Fourier transform) is needed to determine |T , which from (1) is proportional to the generated current |J g . The gradient of P sct flx with respect to the susceptibility degrees of freedom then follows from the defining relation T = (V −1 − G 0 ). (At this step, a nonzero but arbitrarily small scattering potential is assumed to exist at all points within the domain so that V −1 is well defined throughout the ball.) Parametrizing the susceptibility at each pixel by the single material linear equation χ x = t x χ , with t x ∈ (0, 1], Using this result, the total gradient (with well defined t x → 0 limits) is then Hence, only a single additional solve is required to obtain complete gradient information for the optimization objective. For the periodic films, the optimization is carried out using the rigorous coupled wave approach recently presented in Ref. [163].