Hierarchical Mean-Field $\mathbb{T}$ Operator Bounds on Electromagnetic Scattering: Upper Bounds on Near-Field Radiative Purcell Enhancement

We show how the central equality of scattering theory, the definition of the $\mathbb{T}$ operator, can be used to generate hierarchies of mean-field constraints that act as natural complements to the standard electromagnetic design problem of optimizing some objective with respect to structural degrees of freedom. Proof-of-concept application to the problem of maximizing radiative Purcell enhancement for a dipolar current source in the vicinity of a structured medium, an effect central to many sensing and quantum technologies, yields performance bounds that are frequently more than an order of magnitude tighter than all current frameworks, highlighting the irreality of these models in the presence of differing domain and field-localization length scales. Closely related to domain decomposition and multi-grid methods, similar constructions are possible in any branch of wave physics, paving the way for systematic evaluations of fundamental limits beyond electromagnetism.

Accelerating over the last decade, the adoption of inverse design techniques like "density" ("topology") or level-set optimization in photonics [1,2]-ideally matching structural degrees of freedom to the computational grid-has vastly simplified the challenge of discovering geometries with remarkable optical characteristics, leading to improved designs for wide band-gap photonic crystals [3][4][5], enhanced polarization control [6,7], ultra-thin optical elements [8][9][10], and topological materials [11][12][13]. However, because navigating the immense range of allowed structures in such formulations necessitates reliance on local information (approximations based on function evaluations, gradients, etc. [14]) and the relation between fields and structural variations set by Maxwell's equations is non-convex [15], it is rarely known how close these solutions are to true (global) optima, or to what extent they are determined by design choices (e.g. system length scales, material susceptibility, and properties of the algorithm) as opposed to physical principles.
Recently, a number of promising proposals for addressing these knowledge gaps have been put forward by combining Lagrange duality with physical consequences of scattering theory [16][17][18][19]-relaxing the true local constraints of wave physics to global conservation principles, Fig. 1. And in particular, exploiting an optical theorem requiring that real and reactive power be conserved on average [17,20], has been shown to produce performance limits for propagating waves that accurately anticipate the results of density optimization (usually within factors of unity) across a variety of examples [16]. Yet, when these same techniques are applied to situations where evanescent (near-field) wave effects dominate overall behavior, Fig. 2, calculated bounds and objective values obtained by density optimization often differ by several orders of magnitude, and exhibit markedly different trends.
In this Letter, we remedy this issue, and euclidate fundamental connections between dual bounds and struc-tural optimization, by proposing the adoption of constraint hierarchies originating from the central equality of scattering theory-the definition of the T operator. The approach functions, in essence, as a collection of successively refined mean-field theories [21]. At the base of the hierarchy, only the global (spatially integrated) real and reactive power constraints studied in Refs. [16,17] are imposed on the optimization objective, equating these prior results to a first-order mean-field approximation. In every subsequent refinement, the computational domain is decomposed (partitioned) into increasingly small, nested, subdomains, which, through projection, induces additional scattering constraints and results in a higherorder approximation, Fig. 1. Reminiscent of multi-grid and multi-scale methods [22,23], the order of the hierarchy thus acts as a resolution "knob" for systematically controlling the extent to which Maxwell's equations are respected locally, allowing multiple length scales (beyond the size of the domain) to be considered concurrently. Correspondingly, the method also presents a complementary top-down approach to inverse design: the solution of any optimization problem in the constraint hierarchy is always "more optimal" than what is conceptually possible if full wave physics are completely obeyed, whereas inverse design in a finite number of structural degrees of freedom is always suboptimal; in the limit of point subdomains (infinite mean-field order) and infinitesimal structural variations (vanishing "voxels") the two views agree, and, if strong duality holds [24], the bounds solution in fact determines a globally optimal structure. More concretely, proof-of-concept application of the method to the problem of enhancing radiative emission from a dipolar current source in the near field of a structured medium, crucial to optical sensing and quantum information technologies [25][26][27], yields bounds that come substantially closer to the values found by density optimization.
Notation-Throughout, I is used to denote the identity operator, and subscripts on operators (blackboard bold letters) are used as a booking device for the domains and codomains of definition. When only a single subscript is shown, the domain and codomain are identical. The  (1), implicitly contains relations that must be obeyed at every spatial point. Present approaches to electromagnetic limits, however, only impose that this equality be satisfied on average, allowing bound solutions to exhibit highly unphysical local features. (b) By enforcing that the definition of the T operator be respected on successively smaller subdomains or clusters, described by constraints α, the dual bound objective G given in (6) acquires a hierarchical structure that mirrors the primal problem of optimizing some objective L in (4), in terms of an increasing number of structural degrees of freedom contained in the scattering potential V. In the limit of point clusters and complete structural freedom, the two problem statements are equivalent; see discussion in main text. (c) Upper bounds on the plane-wave scattering cross-section σsca, relative to the geometric cross section of the domain σgeo, for any structure bounded by a sphere of radius R. Lighter lines result by enforcing that power be conserved globally, as in Refs. [16,17]. Similarly colored dark lines result when analogous equalities are asserted over eight, evenly spaced, radial subdomains. A profile of one of the density-optimized structures is shown as an inset. Shown below are logarithmic color maps of the corresponding violation in the real part of (2) for one, two and four (evenly spaced) shell clusters, in the plane perpendicular to both the incoming wave vector and direction of polarization.
subscripts b and s mark spatial locations as part of either the background (b) or scattering object (s) within in some predefined domain, Ω. When an operator appears without subscripts, its domain and codomain are Ω. G 0 refers to the background Green's function [16], and hence depends on Ω. V is used to denote the scattering potential, i.e. any polarizable medium not included in G 0 .
Constraints-In scattering theory [28,29], the role commonly played by wave equations (e.g. Maxwell's equations) is typified by the T operator, defined as which, together with knowledge of G 0 and V, abstractly determines all fields for a given source. This operational picture lies at the heart of the hierarchical construction. Any number of manifestly true relations can be generated by probing (1) with linearly independent combinations of fields from the right and linear functionals from the left, and, by Lagrange duality, any set of constraints can be used to produce calculable bounds on any given optimization objective (so long as the dual remains soluble) [24]. Therefore, because every relation derived from (1) is physical, each and every such collection of equalities generates some physical bound on any wave process.
As may be expected, certain choices are more naturally motivated than others, and the difficulty of the associated convex optimization problem that must be solved in each case depends closely on the particular constraints chosen. If all identities contained in (1) are imposed over some complete basis, then the associated primal problem is equivalent to completely free structural optimization, and in computing a solution to the dual (convex) system an actual T operator must be nearly constructed. (The collection of all input-output relationships determines any operator, and so, the only caveat that makes this statement inexact is that the solution of the dual problem may not satisfy every constraint.) Hence, making full use of (1) likely results in an optimization problem comparable to bottom-up structural inverse design [30,31]. On the other hand, if only a select subset of the information contained in (1) is kept, the simplicity of determining bounds can be greatly reduced, at the cost of allowing the discovered solution to inevitably violate some local scattering relation(s). But, unlike standard calculations where the implications of some "partially valid" model on actual solutions are seldom known a priori, solving a dual problem always determines a bound on performance, and therefore always yields useful information.
To begin, it is simplest to work with (1) using either the image or polarization field resulting from the action of T s on a single source field |S , T s |S → |T [16], or, with this image and the action of I s on |S , {T s |S → |T , I s |S → |R } [32]. Letting R = {Ω k } K denote the sets of chosen subdomains, the first choice leads to cluster constraints of the form where U = V −1 † − G 0 † so that Asym [U] is positive def-inite, and F|G = d 3 x F(x) * · G(x) is a complexconjugate inner product (spatial integration over the entire domain). The second choice allows for greater variety, and, (∀Ω k ∈ R), any combination of appear to be, at least presently, sensible.
Both (2) and (3) follow from (1) based on the properties of I s , T s , and the commutativity of spatial projection. (∀ U, V ⊂ R n ) U ∩ V = V ∩ U , and so, as I s and I Ω j both denote projections into spatial locations, I Ω j I s = I s∩Ω j = I s I Ω j , implying that (for any Ω j ⊂ Ω) I s I Ω j = I 2 s I Ω j = I s I Ω j I s , I Ω j T s = I Ω j I s T s = I s I Ω j T s , and I Ω j I s = I Ω j T † s UI s = T † s UI Ω j I s . If only a single cluster corresponding to the entire domain Ω is used, then (2) reduces to the constraints examined in Ref. [16], with the symmetric (Hermitian) and antisymmetric (skew-Hermitian) parts of (2) Every equality of the form given by (2) represents a similar requirement on how power may be transferred between an exciting field and the response (polarization) current it generates. As verified in Fig. 2, these additions are crucial for properly describing near-field interactions. When no constraints other than the global conservation of real and reactive power are imposed, the optimal |T discovered via the method presented in Refs. [16,17] conserves power only by canceling equally large positive and negative violations over Ω, Fig. 1: locally, the power drawn from the incident field is either far greater, or far smaller, than what can be accounted for by absorption and scattering. By successively subdividing the domain, the scale over which such cancellations of unphysical behavior can occur is continually reduced, and because there is always an implicit interaction length scale in U set by material loss and the background Green's function, these tighter requirements ultimately lead to increasingly physical fields. Intuitively, beyond a certain critical size, the effect produced by rapid spatial variations of violation is no different than an averaged "mean" field which respects the constraints at all spatial points.
Though notionally similar, the role of (1) in (3) (third line) requires a distinct interpretation. Rather than describing how power may be transferred between the source and the polarization field it generates within subregions of the domain, these relations state that within the scatterer U † must effectively "invert" |T , reproducing the projection of the source into the scattering object, |R . If (3) were true pointwise, instead of on average over some finite collection of subsets, then both |R and |T would be produced by a geometry with scattering material at all locations x where R(x) = S(x), equating this limit with structural inverse design. (2) and (3), any set of regions in Ω defines a collection of constraints on any observable property of the associated scattering theory; and in turn, these constraints define an optimization problem that bounds the observable, along with a dual solution |T d or {|T d , |R d }. Consider the collection C of all such sets of regions, R = {Ω k } K , with the properties that K Ω k = Ω and (∀k = j)

Hierarchy-Through
giving C a partial ordering. Directly, the map between collections of spatial sets and bounds on a given observable described above, restricted to C, is monotonic. If R ≥ R then the associated bound for R is necessarily tighter (bigger or smaller depending on the objective) than the bound for R. Every refinement R ≥ R in C results in a split of multipliers of the optimization Lagrangian, as each constraint is decomposed into a set of constraints over the matching subregions. Therefore, the codomain of the dual function corresponding to R contains the codomain of the dual function corresponding to R, and so the minima (resp. maxima) of the dual functions maintains the ordering of C, leading to a monotonically tighter bounds, Fig. 1. Hence, successive division of the spatial regions used in generating constraints via (1) indeed yields a well defined hierarchy of T operator bounds, approximating any optimization problem with increasingly better accuracy. Moreover, less refined solutions are always dual feasible points for more refined optimizations, and by evaluating the pointwise versions of (2) and (3) that would hold under complete compliance with the scattering theory for any given dual solution |T d (resp. {|T d , |R d }), as in the heat maps of Fig. 1 and Fig. 2, it is possible to assess where inconsistency is occurring, and use this knowledge to inform further regional decompositions.
The hierarchy construction amounts, conceptually, to a set of successively expanded mean-field theories [33][34][35]. Traditionally, one of the standard ways in which a mean-field theory is constructed is to consider the question of minimizing the Gibb's free energy of a system over the collection of all possible statistical distributions, parameterized by constraints on its moments (expectation values, two-point correlations, etc). To make such problems tractable, the form of the distribution is simplified in some way (e.g. partitioning the true system into effectively interacting spatial clusters) before carrying out a local optimization on the resulting (generally nonconvex) problem. The solution, stationary and selfconsistent with respect to the moments, is called a mean field since to lowest order it describes each component of the system as interacting with one other "averaged" body. By switching to a scattering description, the need for simplification is shifted from the objective to the constraints, but the interpretation of the solution is virtually unaltered. In either case, it is a field that is self-consistent when variations around average values are neglected.
Bounds-Generalizing the program given in Ref. [16], in any order of the hierarchy, the calculation of bounds for any net power transfer (scattering) objective can be equated with a domain monotonic optimization problem described by a sesquilinear Lagrangian. Grouping the source terms for the constraints and objective, |S and |Q , together as the super source |S L α where, supposing an objective of the form Im [ Q|T ] − T| O |T [16], the Z −,− operators are the linear couplings (depending on the multipliers) between the various fields k I Ω k , and α (1) k and α (2) k are the sets of Lagrange multipliers for all symmetric and anti-symmetric constrains respectively. For both scattering and radiative Purcell enhancement, O = Asym V −1 † ; see Ref. [16] for further details. Because optimization is always formulated in terms of real numbers, the total Z matrix is Hermitian.
The unique stationary point of (4) with respect to variations in |T occurs when |T = Z T T −1 Z T S |S = Z T T −1 (Z T S |S + Z T Q |Q ) , and so the dual of (4), G = max F L where the domain F is set by the criterion that max L is finite, is Note that a set α lies within F so long as Z T T −1 is positive definite [36].
Because the dual problem is convex and gradients of the dual reproduce the constraints set by (2), if a stationary point within the feasibility region is found, then strong duality holds and the maximum of the objective is Im [ Q|T ] − T| O |T , with |T determined as above using the multiplier values set by the simultaneous zero point of the gradients (constraints). If no such point exists in F, then the unique minimum value of G attained on the boundary, caused by Z T T becoming semi-definite, remains a bound on the optimization problem.
Example-As an initial exploration of T operator constraint hierarchies, under (2), Fig. 2 illustrates how the introduction of spherically symmetric (concentric) shell clusters alters bounds on Purcell enhancement (radiative emission from a dipolar current source in the vicinity of a structured medium normalized to emission in vacuum) compared to past global-conservation arguments [16,17]. Even with a choice of relatively simple clusters, which leads to considerable numerical simplifications, but also limit the extent to which violation can be localized, two promising trends are observed. First, the number of material and domain size combinations displaying resonant response characteristics, as qualified by a roughly inverse scaling between radiative enhancement and material loss Im [χ], is reduced. To the best of our knowledge, no compact single dielectric material device design exhibits such behavior for values of Im [χ] comparable to those considered here, leading to potentially gigantic gaps between calculated limits and the findings of density optimization, Fig. 2 (a). So long as strong duality is present, which holds for all of (a) and (b) below R/λ ≈ 0.3, the dependence of the cluster bounds on Im [χ] is greatly suppressed; and when a relation between these quantities is observed for Im [χ] 10 −2 , the corresponding optimal polarization field, |T , implies large local constraint violations, as shown by the spatial color maps of the real part of (2) included below Fig. 2 (c). Similarly, Fig. 2 (b) and (c), for separations as small as d/λ = 10 −2 and domains as large as R/λ = 1, cluster-bound values are substantially smaller (often an order of magnitude or more) for strong metallic and dielectric materials, |Re [χ] | 5. Nevertheless, the limits remain several orders of magni-tude larger than what is observed in devices discovered by density optimization when such values of χ are supposed, suggesting further room for improvement. This work was supported by the National Science Foundation under Grants No. DMR 1454836, DMR 1420541, DGE 1148900, EFMA 1640986, the Cornell Center for Materials Research MRSEC (award no. DMR1719875), and the Defense Advanced Research Projects Agency (DARPA) under Agreement No. HR00112090011. The views, opinions and/or findings expressed herein are those of the authors and should not be interpreted as representing the official views or policies of any institution. The authors thank Cristian L. Cortes for useful discussion.