\mathbb{T}-Operator Limits on Optical Communication: Metaoptics, Computation, and Input-Output Transformations
Sean Molesky, Pengning Chao, Jewel Mohajan, Wesley Reinhart, Heng Chi, Alejandro W. Rodriguez
TT -Operator Limits on Optical Communication:Metaoptics, Computation, and Input–Output Transformations S. Molesky, ∗ P. Chao, ∗ J. Mohajan, W. Reinhart, H. Chi, and A. W. Rodriguez Department of Electrical and Computer Engineering,Princeton University, Princeton, New Jersey 08544, USA Department of Materials Science and Engineering,Pennsylvania State University, University Park, Pennsylvania 16802, USA Siemens, Princeton, New Jersey 08540, USA
We present an optimization framework based on Lagrange duality and the scattering T operator ofelectromagnetism to construct limits on the possible features that may be imparted to a collection ofoutput fields from a collection of input fields, i.e., constraints on achievable optical transformationsand the characteristics of structured materials as communication channels. Implications of thesebounds on the performance of representative optical devices having multi-wavelength or multi-port functionalities are examined in the context of electromagnetic shielding, focusing, near-fieldresolution, and linear computing. As undoubtedly surmised since long before Shannon’spioneering work on communication [1], or Kirchhoff’s in-vestigation of the laws governing thermal radiation [2],physics dictates that there are meaningful limits onhow measurable quantities may be transferred betweensenders and receivers (collectively registers ) that applylargely independent of the precise details by which trans-mission is realized. The noisy-coding theorem, for in-stance, proves that probabilistically error free messagepassing is not possible at any rate larger than the “chan-nel capacity” [3]; while more recently, the possibility ofutilizing entanglement as a novel resource has motivatedthe development of a variety of limits on communicationin general quantum systems, e.g. Refs. [4–7]. Broadly,even at the coarsest levels of physical description, thereis generally some notion of invariants, and the existenceof such quantities implicitly precludes the possibility ofrealizing complete engineering control.Equating registers with electromagnetic fields andscattering objects with channels [8], using the languageof communication theory, it is thus reasonable to assumethat there are certain channel-based limits on the extentthat a collection of fields may be manipulated via mate-rial structuring. The wave nature of Maxwell’s equationssets a fundamental relation between wavelength and fre-quency (for propagation) that cannot be arbitrarily al-tered by any realizable combination of material and ge-ometry. Consequently, the scattered fields generated byany true object cannot be matched to any freely selectedmagnitude and phase profile, and so, certain transforma-tions cannot be achieved with perfect fidelity (e.g. knownlimits on light trapping [9–11], cloaking bandwidths [12–14], delay-bandwidth products [15–17], etc.)However, in establishing a means of evaluating thepotential of present and future electromagnetic devicesto address challenges requiring complex functionalities,such as artificial neural networks [18–21], spatial mul- ∗ Equal contribution
FIG. 1.
Investigation schematic and sample applica-tions.
The figure sketches the central question explored inthis article: given a specified design volume(s) and mate-rial(s), how effectively can some particular collection of fieldtransformations (described as input–output pairs) be real-ized? The panel on the left depicts an abstract optical de-vice, the imitable black-box, converting a known set of inputelectromagnetic field profiles into a desired set of output fieldprofiles. The partitions interior to the design domain suggesthow structural degrees of freedom or constraint clusters maybe refined to improve design performance and estimations oflimit performance (bounds), respectively. The smaller imagesto the right depict six (of the many) possible applicationsthat can be easily described within this framework. Workingfrom left-to-right, top-to-bottom, theses illustrations repre-sent a spatial multiplexer, a metaoptic lens, a cloaking shellfor an enclosed object, a light extractor (enhancing radiativeemission from a dipolar source), a waveguide bend, and adirectional antenna. tiplexing [22–24] and computing [25–27], the importantquestion is not whether such fundamental limitations ex-ist, but rather what quantitative implications can be de-duced in a given setting. Beyond well established con-siderations like the need to conserve power in passivesystems and the classical diffraction limit of vacuum [28],such specific channel features are typically unknown [29].Moreover, it is seldom clear what level of possible per-formance improvement could be sensibly supposed, evento within multiple orders of magnitude [30–32]. Archi- a r X i v : . [ phy s i c s . op ti c s ] F e b tectures possessing impressive field-transformation capa-bilities have already been demonstrated, from free-space(grating) couplers [33, 34] to beam steering [35, 36] andpolarization control [37, 38], suggesting that large-scaleoptimization methods may allow scattering attributes tobe tailored to a far greater degree than what has beenseen in past intuition-based designs. Simultaneously,ramifying from the core ideas of Lagrange duality and in-terpreting physical relations as optimization constraintsexpounded below [39–42], a string of recent articles onimproved bounds for scattering phenomena (including ra-diative heat transfer [32], absorbed power [43], scatteredpower [44], and Purcell enhancement [45]) have shownthat, in some cases, only modest improvements over stan-dard designs are even hypothetically attainable [44–53].In this article, we demonstrate that this rapidly devel-oping program for calculating bounds on scattering phe-nomena can be applied to determine limits on the accu-racy to which any particular set of field transformationsmay be implemented, providing an initial exploration ofchannel characteristics in the context of classical elec-tromagnetic (multi-wavelength and multi-port) devices.The article is broken into three sections. In Sec I, a con-densed overview of related work examining channel limitson electromagnetic devices prior to the articles referencedabove is provided. Section II then summarizes the guid-ing optimization outlook of the rest of the article alongwith current methods for constructing T -operator (scat-tering theoretic) bounds. During this brief overview, twoinnovations necessary for handling generic input–outputformulations are introduced: a further generalization onthe variety of operator constraints that can be imposedbased on the definition of the T operator, beyond thepossibility of local clusters examined in Refs. [45, 51, 52],and a generalization of the type of vector image con-straints that should be usually considered, as requiredto enforce that every transformation in an input–outputset references the same underlying structure. Finally,Sec. III provides a number of exemplary applications ofthe theory, including studies of model computational ker-nels, domain shielding and near-field focusing. For thechosen loss function (the quadratic distance or two-normbetween the set of desired target and realizable fields),the gap between calculated bounds and the performanceof device geometries discovered via topology (or “den-sity”) optimization is regularly found to be of unit order.The observed trade-offs between the size of the design do-main, the supposed material response of the device, thespecific transformations (channels) considered, the num-ber of constraint clusters, and the calculated limits revealseveral intuitive trends. Most notably, achievable perfor-mance is strongly dependent on the nature and num-ber of communication channels, with larger device sizesand indices of refraction resulting in greater achievablespatial resolution and transformations, in line with priorheuristics [REF]. For the wavelength-scale systems in-vestigated here, the variability in achievable performancespans multiple orders of magnitude and, as demonstrated in Sec. III, tied to the ability of the formulation to enforcethat all channel functionalities map to a unique (single)optimal geometry. I. CHANNEL LIMITS IN PRIOR ART
Excluding techniques based on Lagrange duality toconstrain electromagnetic design objectives, which arecovered in greater detail in Sec. II, three major threadsstudied in prior art have substantially informed the find-ings and discussion below.
Decomposition —Any finite dimensional linear opera-tor (or more generally any compact operator such as thescattering T operator of electromagnetism) has a corre-sponding singular value decomposition [54]. For almostany example of practical interest, particularly those withdiscrete representations, there are hence correspondingnotions of rank and pseudo rank [8, 55–58]. Rank —Potential field transformations between collec-tions of registers are inherently limited by any boundson maximal rank or pseudo rank. Any set of channelsconnecting registers that do not overlap in space is in-herently limited by the rank (pseudo rank) of its as-sociated free propagation operators (the Green’s func-tion) [30, 44, 59, 60].
Size —The rank (pseudo rank) characteristics of freepropagation in electromagnetics depend explicitly on ge-ometry. That is, there are cases where possible fieldtransformations are strongly limited by the spatial vol-ume occupied by the registers [43, 59].The main difference between past studies of channellimits utilizing these ideas and the T -operator approachgiven in Sec. II lies in the incorporation of additionalphysical constraints concerning the generation of polar-ization currents (within a scattering object) in order toeffect a desired transformation. More concretely, it isoften possible to abstractly describe a device in termsof (possibly intersecting) design and observation regions.Taking a simplistic description of a near-field microscopeas an example, the magnifying lens may be thought ofas a design volume, and the final connection between theoptical components and the electro–optic readout as anobservation region. (As seen in Fig. 2, the design andobservation regions are closely linked to the sender andreceiver registers in optical communication.) Under thisregional decomposition, the end goal of various applica-tions can be understood as minimizing the power differ-ence between the true total field created in the observa-tion region ( (cid:12)(cid:12) E to (cid:11) ) [42], resulting from a known incidentfield ( (cid:12)(cid:12) E i (cid:11) ), and a given target output ( (cid:12)(cid:12) E (cid:5) o (cid:11) ):min | J gd (cid:105) (cid:107) (cid:12)(cid:12) E to (cid:11) − (cid:12)(cid:12) E (cid:5) o (cid:11) (cid:107) . (1)(Here, subscripts mark domains of definitions and super-scripts label different types of fields.) Supposing that po-larizable media exists only within the design region (thedevice), following the notation of the upcoming section, FIG. 2.
Electromagnetic channel limits. (a) The properties of an electromagnetic channel(s) ( C ) depend on the environmentconnecting a particular (collection of) sender ( S ) and receiver ( R ) register(s), including characteristics like the volume eachregister occupies, their spatial separation, and the possibility of a transformation enabling polarizable device ( D ). Channellimits in prior art have been formulated predominately by analyzing the individual operators describing each of these aspects inisolation (e.g., the rank and pseudo rank of the free space Green’s function between the sender and receiver volumes G RS ). Here,in contrast, all environmental factors are treated simultaneously, accounting for the realities of interaction. (b) The difficultyof attaining device performance bounds using the formulation presented in Sec. II, for some predefined material, dependsmainly on three variables: the number of regional constraints (spatial clusters) imposed, the number of channels (simultaneoustransformations) considered, and the physical volume of the design region. The way in which these attributes determine howclosely the computed bounds approximate achievable device performance is generally unknown and problem dependent. Forradiative thermal emission and integrated absorption, increasing domain size leads to increasingly tight bounds when only asingle constraint, the conservation of real power, is employed [43]. For radiative Purcell enhancement, recent results suggestthat predictive bounds for larger design domains require a larger number of constraints [45]. Additional regional constraints,or concurrent analysis of additional channels, always produces tighter bounds. in terms of the generated polarization field (cid:12)(cid:12) J gd (cid:11) Eq. (1)becomes min | J gd (cid:105) (cid:107) (cid:12)(cid:12) E io (cid:11) − (cid:12)(cid:12) E (cid:5) o (cid:11) + i Zk o G od (cid:12)(cid:12) J gd (cid:11) (cid:107) . (2)The minimum of Eq. (2) with respect to (cid:12)(cid:12) J gd (cid:11) , i k o Z (cid:0) G od (cid:1) † (cid:0)(cid:12)(cid:12) E (cid:5) o (cid:11) − (cid:12)(cid:12) E io (cid:11)(cid:1) = (cid:0) G od (cid:1) † G od (cid:12)(cid:12) J gd (cid:11) , (3)illustrates the necessity of the three notions stated above.First, even if Eq. (3) can be satisfied exactly, the mini-mum of Eq. (1) may be nonzero if the span of the tar-get basis of G od does not contain (cid:12)(cid:12) E (cid:5) o (cid:11) − (cid:12)(cid:12) E io (cid:11) . Second,since in extending Eqs. (1)–(3) from single fields to col-lections nothing in the form of the minimum solution isaltered, the rank (pseudo rank) of the Hermitian operator (cid:0) G od (cid:1) † G od sets a limit on the number of channels thatcan connect (interact) with any register contained in theobservation domain. Third, as highlighted by the explicitdomain subscripts, the Green’s function ( G od ) connect-ing the observation and design domains, on which thetwo previous points are based, depends on the geometricfeatures of both regions.While of obvious importance, it is clear that these ob-servations do not fully encompass the difficulties associ-ated with the larger question of realizing a given set of field transformations. Directly, without additional con-straints, the challenge associated with attaining a par-ticular polarization current | J gd (cid:105) within the design regioncannot be inferred from Eq. (3). To include this infor-mation into the minimization of Eq. (1), additional phys-ical constraints must be taken into account. Achievingthis aim, without making the resulting problem state-ment computationally infeasible, is precisely the goal ofthe T -operator framework presented below. II. TECHNICAL DISCUSSIONA. Scattering theory, Lagrange duality, and limits
Following the sketch provided in Ref. [29] (support-ing details can be found in Refs. [61–65]), there are twosuperficial distinctions that distinguish a scattering the-ory view of electromagnetics from Maxwell’s equations.First, supposing a local, linear scattering potential, thefundamental wave relation from which all field propertiesare derived is I s = I s (cid:0) V − − G (cid:1) T s , (4)where the s subscript denotes projection into the scat-tering object, G is the background Green’s function ofthe design domain Ω s , throughout taken to be the vac-uum Green’s function (truncated to a problem-specificvolume) scaled by k o = (2 π/λ ) (so that all spatial di-mensions may be simply defined relative to the wave-length), V is the bulk (spatially local) scattering potentialof the material forming the design, I is the identity oper-ator, and T s is the scattering operator, formally definedby Eq.(4) [66] [67]. (Although more general descriptionsare possible, it is implicitly assumed throughout this ar-ticle that V = I χ , where χ is electric susceptibility ofsome isotropic medium.) Second, one of either the elec-tromagnetic field or polarization field is singled out asan initial [68], and all other fields are generated via thescattering operator T s . This leads to a division of thetheory into two classes: initial flux problems, where afar-field incident flux (cid:12)(cid:12) E i (cid:11) is assumed, generating a po-larization current | J g (cid:105) = − ( ik o /Z ) T s (cid:12)(cid:12) E i (cid:11) and scatteredfield (cid:12)(cid:12) E s (cid:11) = ( iZ/k o ) G (cid:12)(cid:12) J g (cid:11) , and initial source problems,where a polarization source (cid:12)(cid:12) J i (cid:11) is assumed, generat-ing an electromagnetic field | E g (cid:105) = ( iZ/k o ) G T s V − (cid:12)(cid:12) J i (cid:11) and total current (cid:12)(cid:12) J t (cid:11) = T s V − (cid:12)(cid:12) J i (cid:11) .Given that most quantities related to power trans-fer result from inner products between electromagneticand polarization fields [61], many common design objec-tives within scattering theory, encompassing applicationsvarying from enhancing the amount of radiation that canbe extracted from a quantum emitter [69–71] to increas-ing absorption into a photovoltaic cell [72–74], are givenby sesquilinear forms on the scattering operator T s . Asthe central example investigated in this article, continu-ing the discussion of Sec. I and covering all the specificapplications treated in Sec. III, suppose that a collectionof fluxes (cid:8)(cid:12)(cid:12) E ik (cid:11)(cid:9) , with k ranging over some finite indexingset, is incident on a scattering design region Ω. Assumethat the goal of the (input–output) device is to create aset (cid:110)(cid:12)(cid:12) E (cid:5) k (cid:11)(cid:111) of target fields within an observation regionΩ o , marked by a domain subscript o and spatial projec-tion operator I o , i.e. a collection of register mappings (cid:12)(cid:12) E ik (cid:11) (cid:55)→ (cid:12)(cid:12) E (cid:5) k (cid:11) . Taking the Euclidean norm as a metric ofdesign performance,dist (cid:16)(cid:110)(cid:12)(cid:12) E tk (cid:11)(cid:111) , (cid:110)(cid:12)(cid:12) E (cid:5) k (cid:11)(cid:111)(cid:17) = (cid:115)(cid:88) k (cid:107) I o (cid:0)(cid:12)(cid:12) E tk (cid:11) − (cid:12)(cid:12) E (cid:5) k (cid:11)(cid:1) (cid:107) , (5)the equivalent objective (in the sense (5) and (6) havethe same minimum)Obj (cid:16)(cid:110)(cid:12)(cid:12) E tk (cid:11)(cid:111) , (cid:110)(cid:12)(cid:12) E (cid:5) k (cid:11)(cid:111)(cid:17) = (cid:88) k (cid:107) I o (cid:0)(cid:12)(cid:12) E tk (cid:11) − (cid:12)(cid:12) E (cid:5) k (cid:11)(cid:1) (cid:107) (6) is effectively a sesquilinear form on T s .Obj = (cid:88) k (cid:10) E (cid:5) k (cid:12)(cid:12) I o (cid:12)(cid:12) E (cid:5) k (cid:11) − (cid:2)(cid:10) E (cid:5) k (cid:12)(cid:12) I o (cid:12)(cid:12) E tk (cid:11)(cid:3) + (cid:10) E tk (cid:12)(cid:12) I o (cid:12)(cid:12) E tk (cid:11) = (cid:88) k (cid:10) E (cid:5) k (cid:12)(cid:12) I o (cid:12)(cid:12) E (cid:5) k (cid:11) − (cid:2)(cid:10) E (cid:5) k (cid:12)(cid:12) I o (cid:12)(cid:12) E ik (cid:11)(cid:3) + (cid:10) E ik (cid:12)(cid:12) I o (cid:12)(cid:12) E ik (cid:11) + (cid:10) E ik (cid:12)(cid:12) T † s G † I o G T s + 2 Sym (cid:2) I o G T s (cid:3) (cid:12)(cid:12) E ik (cid:11) − (cid:2)(cid:10) E (cid:5) k (cid:12)(cid:12) I o G T s (cid:12)(cid:12) E ik (cid:11)(cid:3) , (7)and the objective of the design amounts to a minimiza-tion ofObj ( T s ) = (cid:88) k (cid:10) E ik (cid:12)(cid:12) T † s G † I o G T s + 2 Sym (cid:2) I o G T s (cid:3) (cid:12)(cid:12) E ik (cid:11) − (cid:2)(cid:10) E (cid:5) k (cid:12)(cid:12) I o G T s (cid:12)(cid:12) E ik (cid:11)(cid:3) . (8)Here, Sym [ O ] and Asym [ O ] denote symmetric (Her-mitian) and anti-symmetric (skew-hermitian) parts of O . (Other objective forms, e.g. originating from othernorms, could also be considered, but (8) offers many use-ful simplifications for both analysis and computation.)Acting with T † s from the left on (4), giving T † s = T † s U † T s with U † = V − − G so that Asym [ U ] is positive definite,shows that the physical constraints of scattering theorycan also be written as similar sesquilinear forms. As such,any electromagnetic design objective of this input–outputform (Ref. [42]) can be cast as a quadratically constrainedquadratic program [75, 76] (QCQP) by treating T s as avector, i.e.min T s Obj ( T s )such that ( ∀ l, m ) (cid:10) B l (cid:12)(cid:12) T s (cid:12)(cid:12) B m (cid:11) = (cid:10) B l (cid:12)(cid:12) T † s UT s (cid:12)(cid:12) B m (cid:11) , (9)over some basis [77] for the space of electromagneticfields over the design region Ω s . Although QCQP prob-lems are part of the non-deterministic polynomial-time(NP) hard complexity class, a number of relaxation tech-niques, as well as a range of well developed solvers,are known to often yield accurate solution approxima-tions [78] (that may in fact be exact). Principally, let-ting L ( T s , λ ) = Obj ( T s ) + λ C ( T s ) be the Lagrangianof Eq. (9), with C ( T ) denoting the collection of imposedscattering constraints and λ the associated set of La-grange multipliers, it always possible to relax Eq. (9),or any differentiable optimization problem for that mat-ter [79], to a convex problem by considering the uncon-strained dual optimization problemmax λ G ( λ ) , (10)with G ( λ ) = min T s L ( T s , λ ). Because any true solutionof Eq. (9) requires every constraint relation to be satis-fied, for any collection of multipliers λ , G ( λ ) is alwayssmaller (resp. larger if the optimization objective is tomaximize some quantity) than any solution of the primal(original) optimization. Thus, in addition to providing ameans of obtaining approximate solutions to Eq. (9) [78],solving Eq. (10) also determines a bound, or limit, onthe performance that could be possibly achieved by anyrealizable device geometry. By varying the constraintsincluded in (9) these bounds can be tailored to encom-pass select attributes of the channel(s) as one wishes,such as qualifications on material composition (via V ),the available design domain (via truncation of G ), andthe actual input fields (via (cid:110)(cid:12)(cid:12) E ik (cid:11)(cid:111) ). Further, althoughwe have not seriously examined the prospect, it is alsopossible to apply similar reasoning to objective forms dif-fering substantially from Eq. (6). B. Mean-field hierarchy (spatially localizedconstraints)
Building on the above perspective, as described ingreater detail in Refs. [45, 53], the T operator providesa natural means of obtaining arbitrarily accurate, con-text specific, mean-field solution approximations. Theoperator relation given by Eq. (4) is equally valid underevaluation with any linear functional or vector, or com-position with any other maps, and can hence be “coarse-grained” (averaged). Expressly, to construct an appro-priate mean-field approximation, take P Ω c to denote ageneralized spatial projection operator into the spatialcluster (subregion) Ω c , with the added freedom that, atany point within Ω c , P Ω c may transform the field in anyway. That is, taking the three dimensional nature ofthe electromagnetic and polarization vector fields intoaccount, P Ω c may be any operator of the form P Ω c ( x , y ) = (cid:40) (cid:54) = y or x (cid:54)∈ Ω c M x x = y and x ∈ Ω c , (11)where M x , at any point in Ω c , is some any linear operator(3 × x . Recalling that I s is defined as I s ( x , y ) = (cid:40) (cid:54) = y or x (cid:54)∈ Ω s = y and x ∈ Ω s , (12)with the identity operator (3 × I s , regardless of the actual ge-ometry of the scatterer, allowing Eq. (4) to be rewrittenas T † s P † Ω c I s = T † s P † Ω c I s (cid:0) V − − G (cid:1) T s ⇒ T † s I s P † Ω c = T † s P † Ω c I s (cid:0) V − − G (cid:1) T s ⇒ T † s P † Ω c = T † s P † Ω c (cid:0) V − − G (cid:1) T s ⇒ P Ω c T s = T † s UP Ω c T s . (13)Working with increasingly refined collections of clusters,sets of subdomains { Ω c } with c ranging over some index-ing set such that ∪ c Ω c = Ω, Eq. (13) leads to increasingly refined mean-field approximations that mirrors the “bot-tom up” optimization of a given objective with respectto structural degrees of freedom [80]. Letting | E (cid:105) de-note some predefined electromagnetic flux, and makingthe simplifying choice that M x = for all x ∈ Ω c actingon Eq. (13) with (cid:10) E (cid:12)(cid:12) . . . (cid:12)(cid:12) E (cid:11) gives (cid:90) Ω c d x E ∗ ( x ) · T ( x ) − (cid:90) Ω d y T ∗ ( y ) · U ( y , x ) · T ( x ) = (cid:10) E (cid:12)(cid:12) P Ω c (cid:12)(cid:12) T (cid:11) − (cid:10) T (cid:12)(cid:12) UP Ω c (cid:12)(cid:12) T (cid:11) = 0 , (14)where the generated polarization current (cid:12)(cid:12) T (cid:11) is the imageof (cid:12)(cid:12) E (cid:11) under the action of T s ( T s (cid:12)(cid:12) E (cid:11) (cid:55)→ (cid:12)(cid:12) T (cid:11) ). Returningto Eq. (14), the selection of any cluster Ω c therefore de-fines an averaging kernel for the scattering theory: overΩ c the otherwise free field must respect the true physi-cal equations of electromagnetism, however, at any point x ∈ Ω c unphysical fluctuations of the field (deviations ofthe integrand from zero) are permitted [81].As any physical polarization current field (cid:12)(cid:12) J g (cid:11) will au-tomatically satisfy any relation of the form of Eq. (14),any optimization that only imposes such relations over afinite set of clusters { P k } will automatically possess ob-jective values at least as optimal as what can be achievedif Eq. (4) is imposed in full. Hence, solving the mean-fieldoptimizationmin / max | T (cid:105) Obj (cid:0)(cid:12)(cid:12) T (cid:11)(cid:1) such that ( ∀ Ω k ) (cid:10) E (cid:12)(cid:12) P Ω k (cid:12)(cid:12) T (cid:11) − (cid:10) T (cid:12)(cid:12) UP Ω k (cid:12)(cid:12) T (cid:11) , (15)similar to the mean-field theories used in statisticalphysics [82, 83], game theory [84, 85] and machine-learning [86, 87], results in a bound (limit) on physicallyrealizable performance [88]. Due to the implicit wave-length scale contained in the Green’s function ( G in U ),there is, for most objectives, no meaningful differencebetween a mean-field fluctuating sufficiently rapidly anda physical field, satisfying Eq. (4) exactly. As a conse-quence, smaller clusters, effectively, lead to higher ordermean-field theories that more closely bound what is actu-ally possible. This tightening with decreasing cluster sizecauses (15) to act as an intriguing complement to stan-dard structural design. In “bottom-up” geometric opti-mization, the introduction of additional degrees of free-dom opens additional possibilities that improve achiev-able optima. In “top-down” mean-field optimization, theintroduction of additional clusters results in additionalconstraints that reduce the space of possibilities availableto the design field (cid:12)(cid:12) T (cid:11) . Our present high-level approachto solving Eq. (15) is described in Ref. [45]. C. Multiple transformations (channels)
When handling a collection of sources that span a rel-atively small subspace, an efficient approach to comput-ing T operator bounds through Eqs. 8 & 13 is to workwith individual inputs and outputs. Take {| S k (cid:105)} to bea given collection of sources, and let {| T k (cid:105)} be the col-lection of polarization fields resulting from the action of T s , T s | S k (cid:105) (cid:55)→ | T k (cid:105) . Following Ref. [45], for any pair ofindices (cid:104) k , k (cid:105) , and each P Ω c cluster operator, | S k (cid:105) and {| T k (cid:105) , | T k (cid:105)} must obey the relation (cid:104) S k | P Ω c | T k (cid:105) = (cid:104) T k | UP Ω c | T k (cid:105) , (16)where, again, U = V − † − G † , and (cid:104) F | G (cid:105) denotes thestandard complex-conjugate inner product over the com-plete domain ( (cid:104) F | G (cid:105) = (cid:82) Ω d x F ( x ) ∗ · G ( x ) in the spatialbasis). In Eq.(16), the extension to pairs of sources andpolarization fields, compared to the single source con-straints examined in Ref. [45], is necessary to accountfor the fact that a single scattering object (structuredmedia) simultaneously generates each | T k (cid:105) from each | S k (cid:105) . Over some set of clusters, using only the “diag-onal” constraint between each source field ( | S k (cid:105) ) andpolarization response ( | T k (cid:105) ) is equivalent to consideringeach source independently. But, by additionally includ-ing “off-diagonal” interactions between pairs, Eq. (16)introduces the requirement of a single consistently de-fined scattering object. Supposing a spatial basis in thelimit of “point” (vanishingly small) clusters and completefield mixing, Eq. (16) becomes ( ∀ x ∈ Ω & (cid:104) k , k (cid:105) ) S ∗ k ( x ) · M x · T k ( x )= (cid:90) Ω d y T ∗ k ( y ) · U ( y , x ) · M x · T k ( x ) . (17)Therefore, making use of the possible freedom in the def-inition of M x , whenever any T k n ( x ) is nonzero any ofthe three local vector coordinates we must have S ∗ k ( x ) = (cid:90) Ω s d y T ∗ k ( y ) · U ( y , x ) , (18)for all k and x ∈ Ω s , the collection of all spatial pointsin Ω where some T k ( x ) (cid:54) = . Taking the adjoint, thetotality of the point constraints given by Eq. (18) can becodified as | T k (cid:105) = (cid:16) U † Ω s (cid:17) − | S k (cid:105) . (19)Comparing with Eq. (4), Ω s thus determines the uniquegeometry of the scattering object (structured medium)for all sources. III. APPLICATIONS
In this section, the ideas presented in Sec. II are ap-plied to several example applications, expounding the useof input–output language to describe physical design asQCQP optimization problems, and demonstrating theutility and tightness of the associated bounds. Three main overarching features, which credibly apply beyondthe actual scope studied, are seen. First, distinct chan-nels can display widely varying (often unintuitive) char-acteristics. Second, both in terms of computed boundsand the findings of inverse design, attainable performancemay vary greatly as the number of channels increases—the design freedom offered by spatial structuring in awavelength scale device is by no means inexhaustible.Third, increasing the number of constraint clusters canimprove the accuracy of bounds calculations, but thisimprovement is problem dependent and varies, in par-ticular, with the material assumed. Throughout, as nofeature size can be imposed if one wishes to set boundson any possible geometry, quoted inverse design valuesrepresent the performance of “grayscale” structures, al-lowing each pixel to take on any susceptibility value ofthe form tχ with t ∈ [0 , A. Field screening
Objective —Minimize the spatially integrated field in-tensity over a subwavelength ball for two dipole fields,supposing that the distribution of scattering material (ofa predefined susceptibility) must be contained in a thinsurrounding shell.The observation (shielded) region is taken to have aradius of 0 . λ , while the design domain is confined withina shell of inner radius 0 . λ and outer radius 0 . λ . Assketched in the inset of Fig. 3(a), the “separation vectors”connecting the center of each dipole sources to the centerof the observation domain are supposed to be aligned.The loss objective is defined asLoss (cid:0)(cid:12)(cid:12) E ta (cid:11) , (cid:12)(cid:12) E tb (cid:11)(cid:1) = (cid:107) I obs (cid:12)(cid:12) E ta (cid:11) (cid:107) + (cid:107) I obs (cid:12)(cid:12) E tb (cid:11) (cid:107) (cid:107) I obs (cid:12)(cid:12) E ia (cid:11) (cid:107) + (cid:107) I obs (cid:12)(cid:12) E ib (cid:11) (cid:107) (20)with (cid:107) . . . (cid:107) denoting the Euclidean two-norm and I obs projection into the observation domain. Two clusters,evenly dividing the radius of the design region (0 . → . λ and 0 . → . λ ), are used in all bound compu-tations.In the limit of large separation, the dipole fields withinthe design and observation domains closely resemblecounter-propagating planewaves. Hence, the constantasymptote value seen as d approaches λ is, for all prac-tical purposes, a bound on shielding for two opposingplanewaves. The decrease of the bounds for small sep-aration, speculatively, is caused by a combination of lo-calization effects. As the dipoles are brought into thenear-field of the design, the bulk of the field magnitudeincreasingly shifts towards the edge of the observation re-gion. Consequently, inducing small shifts in the positionof the field results in a comparatively larger field expul-sion. The presence of relatively larger evanescent fieldslikely also simplifies the suppression of total emission viathe (inverse) Purcell effect.The imposition of cross-constraints (solid lines), whichin the limit of point clusters enforce that all considered FIG. 3.
Dipole screening and field transformations.
The figure depicts two representative three-dimensional ap-plications of input–output bounds: limits on the ability ofany structured geometry confined to a thin shell to (a) screendipole fields, and (b) mask the separation of dipole sourcepairs within the encompassing environment. In both sit-uations, dashed lines mark independent bounds, calculatedwithout the inclusion of cross-constraints which enforce aunique optimal structure across all field transformations,while the solid lines result when source interactions effectsare taken into account. A material loss value of Im [ χ ] = 0 . channels are realized within one common geometry (seeSec. II C), is seen to produce a relatively small correctionto the bounds attained for independently optimized chan-nels (dashed lines). This suggests that although perfectshielding is not possible at sufficiently large dipole sepa-ration, given the small shell and material constraints in-vestigated, many structures may potentially achieve sim-ilar performance, i.e. the considered channels have yetto “set” the scattering structure in any meaningful way. B. Dipole resolution and masking
Objective —Manipulate the field emanating from a pairof axially aligned dipoles, again limiting the volume avail-able for device design to a thin encompassing sphericalshell, so that exterior to the shell it appears that thedipoles are separated by either a larger or smaller dis-tance. If the field is altered so that the separation dis-tance appears larger than it actually is, it is easier toresolve that there are in fact two dipoles present, insteadof a single dipole of some effective strength.More exactly, the input field is generated by a pair of z polarized dipole sources, labeled 1 and 2, situated at+ d i ˆ z and − d i ˆ z for i = 1 ,
2, and the target field is taken toarise from equivalent dipolar configuration with a distinctseparation, ± d t ˆ z . The observation domain is set to be ashell of half wavelength thickness enveloping the designdomain, which extends from an inner radius of 0 . λ toan outer radius of 0 . λ . Bounds, in all cases, are com-puted using ten spherical clusters, evenly spanning thedesign domain with respect to radius. Additional techni-cal details concerning the example are given in Sec. V A.The bounds depicted in the three heat maps of Fig. 4,exploring variations in both input d = d and output d t pair separations, indicate that the challenge of ei-ther reducing or enlarging the perceived separation ofthe source varies sizably with supposed material prop-erties and gap distances. For larger target or sourceseparations (dipoles approaching the interior surface ofthe design domain) the increasingly dominant role ofevanescent fields makes almost any transformation chal-lenging, especially in weakly polarizable media. Focusingon Fig. 4(b), where the real part of the permittivity islimited to small positive numbers (weak dielectrics), thebounds prove that neither separation manipulation be-havior is achievable in any practical sense when the inputand target separations vary by more than ≈ . λ . Con-versely, even under the imposition of considerably largermaterial loss (Im [ χ ]), panel (a) suggests that much betterperformance is achievable with metals (Re [ χ ] < − d while keeping d = 0 . λ fixed, fora common target field corresponding to two co-locateddipoles at the center of the ball, d t = 0. The findingshighlight the need to enforce that performance bounds formultiple channels refer to a unique scattering geometry,achieved via the presence of cross-constraints in the opti-mization (see Sec. II C). When the two channels are ana-lyzed separately (dashed lines), the bound values match FIG. 4.
Dipole resolution and masking.
Considering the same three-dimensional system depicted in the inset of Fig. 3(b),but with a single dipole-pair source, the three heat maps depict bounds on the ability of an arbitrary geometry made froma specified material, restricted to a thin shell, to modify the perceived separation of a pair of axially aligned dipole sources,as implied by the electric field exterior to the design. Setup details are given in the main text. In moving from the trivialdiagonal (no transformation) to the top left, the target field is set to be the field of an axial dipole pair further separated thanthat of the source. Opposingly, in moving to the bottom right, the target field is set to be the field of an axial dipole pair lessseparated than that of the source. The observed asymmetry between these two directions provides evidence of the ostensiblygreater challenge of resolving closely separated fields, compared to masking separation. Field profiles for representative dipolepairs, along with the average field intensity within the observation domain (dashed lines), are given in panel (d). those displayed on the heat maps of Fig. 4, with thecontribution of source 1 to the loss function approach-ing zero as d →
0, since in this regime the source andtarget fields match. However, when cross-constraints areincluded, the lack of a need for a design in this regimefor the first source is at odds with what is needed to con-vert the field of the widely separated d source, leading todegraded performance. As d approaches d , the require-ment of a unique gemometry no longer introduces addi-tional requirements, and thus the “simultaneous” (cross-constrained) and “independent” results merge. C. Math kernels (integration and differentiation)
Objective —Within a bounding rectangle, design a two-dimensional scattering profile, for a known collection ofincident waves, such that the relation between each inci-dent field and total field, along specified input and out-put observation planes, reproduces the action of Volterra(integration) or differentiation operator. Given that theVolterra and differentiation kernels are two of the basicelements of differential equations, these examples are ofconsiderable interest in relation to recent proposals foroptical computing [25–27].The design region is chosen to have a transverse length (thickness) of 0 . λ and parallel length of 2 λ , insetsFig. 5(a) and Fig. 5(c). The input observation planeis set at distance of 0 . λ from the input edge of thedesign region, and the output observation plane is lo-cated at a distance of 0 . λ from the output edge ofthe design region. Incident fields are generated by dipole(pixel) sources placed on a source plane, running paral-lel to the long edge of the design region, situated at adistance 0 . λ from the input observation plane. Thesource plane has length of 1 . λ and is centered in rela-tion to the design domain, creating a 0 . λ offset fromthe top and bottom edges of the design region and in-put and output planes. For the Volterra operator, thetarget fields along the output plane, letting y denote thedirection running parallel to the long edge of the designregion, are calculated as E (cid:5) ( x (cid:5) , y ) = (cid:82) y d y (cid:48) E i ( x i , y (cid:48) ),with x (cid:5) and x i denoting the common coordinate of theoutput and input plane, respectively. For the differ-ential operator, keeping the above notation, the targetfields are calculated as E (cid:5) ( x (cid:5) , y ) = ∂ E i ( x i , y ) /∂y , im-plemented as a finite-difference approximation on a com-putational grid. The optimization objective is definedby Eq. (7) and the bounds formulated following the de-scription given in Sec. II. Further implementation detailsappear in Sec. V B. The loss function values appearing in FIG. 5.
Volterra and differentiation kernels.
Following the description in the main text, the figure depicts limits onthe degree to which a structured medium of a given electric susceptibility χ , confined to a rectangular design domain, mayimplement Volterra (top row) and differentiation (bottom row) kernels, as as function of the size of the subbasis (number ofchannels or transformations) over which the operator is defined. The cluster labels reference the number of cluster constraintsimposed on the design domain, following the discussion of Sec. II and with further details given in the main text. Notably,the calculated bounds generally agree to within an order of magnitude of the performance of inverse designed structures wheneither the rank of the specified subbasis, or number of imposed clusters, is large. In all four cases, a material loss value ofIm [ χ ] = 0 .
01 is assumed. Representative topology-optimized structures are included as insets in (b) and (d). The performanceof these binarized systems is roughly a factor of two worse than the associated “grayscale” geometries compared to in the figure.
Fig. 5 are determined byLoss (cid:16)(cid:12)(cid:12) E t (cid:11)(cid:17) = (cid:107) I out (cid:16)(cid:12)(cid:12) E t (cid:11) − (cid:12)(cid:12) E (cid:5) (cid:11)(cid:17) (cid:107) (cid:107) I out (cid:12)(cid:12) E i (cid:11) (cid:107) + (cid:107) I out (cid:12)(cid:12) E (cid:5) (cid:11) (cid:107) . (21)In Eq. (21), I out denotes projection onto the output ob-servation plane, and an overline on a vector indicatesvertical concatenation over all like named fields indexedby the incident waves, i.e. (cid:12)(cid:12) E t (cid:11) = (cid:2)(cid:12)(cid:12) E t (cid:11) , . . . , (cid:12)(cid:12) E t N (cid:11)(cid:3) for N sources. The cluster numbers appearing in Fig. 5 re-fer to specific, even and consistent, divisions of the designdomain rectangle along its long and short edges. Writingdivision of the short edge first, 8 = 2 ×
4, 16 = 2 × × ×
12. The most striking aspect of the findings displayed inFig. 5, especially when viewed in conjunction with Fig. 6,is the widely varying characteristics observed betweenthe four cases. For the Volterra and differentiation ker-nels, at least within the subbasis used here, the num-ber of channels considered has a remarkable influence onthe observed limit values. Given that analogous trendsare seen in performance of the associated structural opti-mizations, it can be safely concluded that this behavior isfundamental, hinting that there is a sort of Fourier (res-olution) limit at play (e.g., size-dependent constraints onspace-bandwidth products are known to limit the degreeto which a given optical feature may be resolved with afinite device [89, 90]). Moreover, presumably because dif-0ferentiation demands creation of more rapid profile varia-tions whereas Volterra operations demand only smooth-ing, loss function values (for both bounds and inversedesign) differ by multiple orders of magnitude for smallnumbers of channels ( ≈
3) before reaching “asymptotes”that differ roughly by a factor of 10, for χ = 10 . i . χ ]). For Re [ χ ] = 2, subdivid-ing the domain generally leads to substantial tightening,while for Re [ χ ] = 10 many more clusters are needed toproduce meaningful alterations owing to the decrease inthe effective wavelength in the medium. D. Near-field lens (metaoptics)
Objective —Working within a given bounding rectan-gular design domain, maximize the average field inten-sity within a specified focal region for a predefined set ofincident plane waves spanning a cone of incidence.The size of the rectangular design region is character-ized by the lengths L d = 0 . λ and L d = 2 λ , inset Fig. 6.The focal region is a square centered along L d , with sidelength L f = L f = 0 . λ . The middle of the focal regionis set to reside at a distance of 0 . λ from the outgoingedge of design region. The design objective isObj (cid:16)(cid:110)(cid:12)(cid:12) E tk (cid:11)(cid:111)(cid:17) = max N (cid:88) i =1 (cid:10) E ti (cid:12)(cid:12) I focal (cid:12)(cid:12) E ti (cid:11) (22)where the summation index runs over the N angles ofincidence, evenly distributed over a cone of ( N − × ◦ degrees centered on the midpoint of L d . As before, thecluster labels given in Fig. 6 reference even divisions ofthe rectangular design domain along its short and longedges: 16 = 2 ×
8, 32 = 4 ×
8, 36 = 3 ×
12, 100 = 5 × × FIG. 6.
Focusing metalens.
Highlighting the applicabil-ity of the proposed method to establish limits on metaopticalelements, the figure provides bounds on the ability of anystructure confined within a rectangular domain of subwave-length thickness to enhance the average field intensity overa predefined focal region in the near field of the structure,for a given collection of incident planewaves. Also shown arecomparative performances of inverse designed structures. Amaterial loss value of Im [ χ ] = 0 .
01 is supposed in all cases.Binarized structures for 3 , , . , . , .
56 and 1 . largely responsible for these findings. First, the focusingproblem does not impose any specific target profile forthe output fields, inset Fig. 6(b), leading to a greater va-riety of geometric possibilities with similar focusing per-formance. Second, given the assumed design domain, thespecified location and volume of the focal region makesthe problem more difficult than what might be naivelyexpected [91]; quite different results may be encounteredif the “spot size”, or separation between the focus anddesign domain, are reduced. Third, the extended natureof planewave sources naturally suggests the uniform clus-1ter constraints that have been employed. Contrastingly,for the Volterra and difference kernels, it is plausible thatnon-uniform cluster distributions, or simply smaller clus-ters, could lead to tighter bounds. A systematic study ofsuch questions, once a more resource-effective approachfor solving Eq. (9) is found [53], remains as an importantdirection for future work. IV. SUMMARY DISCUSSION
In summary, we have shown that recently developedmethods for calculating limits on sesquilinear electromag-netic objectives [42, 45, 51–53], can be applied to a widevariety of design problems, including spatial multiplex-ing [22–24], metaoptics [92–95], linear computation [25–27] and light extraction [96–99]. Employing the lan-guage of communication theory, this program stands asa significant extension of prior work on channel-basedelectromagnetic limits [8, 30, 55–60]. Namely, althoughhighly insightful, the characteristics of the backgroundGreen’s function connecting the volumes containing par-ticular sender and receiver registers are generally insuffi-cient for accurately assessing whether the extent to whichsome desired communication can occur, particularly inwavelength-scale devices. Rather, as may be confirmedby a survey of the results presented in Sec. III, the degreeto which communication between a predefined collectionof registers can occur may depend strongly on a range ofother environmental factors, such as the physical size andresponse parameters available for designing the channels,and the spatial profile of the register fields.The strong correspondence observed between boundscomputed with this approach and the findings of in-verse design exemplified in Sec. III continues many ofthe trends seen in the earlier works cited above. Regard-less of the particular objective and distribution of clusterconstraints considered, we have yet to encounter a sit-uation in which strong duality does not hold [100]. Assuch, it would seem that the outstanding difficulties ofconstructing a general toolbox for realizing (tight) limitsthat incorporate the full wave physics and limitations ofMaxwell’s equations, and in turn offer guidance for prac-tical designs, are computational [29, 45, 53]. All presentevidence indicates that optimization problems followingthe form of Eq. (15) can be solved exactly via the dualityrelaxation described in Sec. II. If this is indeed the case,then the matter of central importance is not whetherthere is some configuration of clusters that will captureall key effects that physically limit communication, but whether some reasonably good approximation can actu-ally be solved in an acceptable amount of time. Withthis in mind, we believe that the results presented aboveprovide substantial reason for optimism.Finally, the specificity of the program presented in themain text raises a prescient question pertaining to re-alizing abstract transformations. Expressly, there aremany instances in which the primary concern is how wellthe transformation can be implemented between a givennumber of undetermined inputs and outputs (channels).The distinction amounts, in essence, to the same shift inperspective advocated for in recent “end-to-end” inversedesign formulations [95, 101] and waveform optimizationinvestigations [102, 103]: the exact characteristics of thedomain and codomain sets for some particular discrim-ination or transformation are often mutable, and onlyimportant in so far as they implicitly alter achievableperformance. The duality program outlined above canbe extended to treat this problem in a variety of ways,and we intended to present this analysis in detail in anupcoming work.
ACKNOWLEDGMENTS
This work was supported by the National ScienceFoundation under the Emerging Frontiers in Researchand Innovation (EFRI) program, EFMA-1640986, theCornell Center for Materials Research (MRSEC) throughaward DMR-1719875, and the Defense Advanced Re-search Projects Agency (DARPA) under agreementsHR00112090011, HR00111820046 and HR0011047197.The views, opinions and findings expressed herein arethose of the authors and should not be interpreted asrepresenting official views or policies of any institution.
V. APPENDIXA. Off-origin dipole spherical wave expansions
Following results of Ref. [104], as explained in Refs. [64,105], off-origin, on-axis, outgoing spherical waves canbe expanded in terms of on-origin outgoing and regularspherical waves as W out plm ( r ) = (cid:80) p (cid:48) ,l (cid:48) U in ± p (cid:48) p,l (cid:48) lm ( d ) W reg p (cid:48) l (cid:48) m ( r ± d ˆ z ) r < d (cid:80) p (cid:48) ,l (cid:48) U out ± p (cid:48) p,l (cid:48) lm ( d ) W out p (cid:48) l (cid:48) m ( r ± d ˆ z ) r > d , (23)where2 U in ± p (cid:48) p,l (cid:48) lm ( d ) = (cid:88) ν (cid:20) l ( l + 1) + l (cid:48) ( l (cid:48) + 1) − ν ( ν + 1)2 δ pp (cid:48) ∓ i mdk o (1 − δ pp (cid:48) ) (cid:21) A h ± l (cid:48) lνm ( d ) , (24) U out ± p (cid:48) p,l (cid:48) lm ( d ) = (cid:88) ν (cid:20) l ( l + 1) + l (cid:48) ( l (cid:48) + 1) − ν ( ν + 1)2 δ pp (cid:48) ∓ i mdk o (1 − δ pp (cid:48) ) (cid:21) A j ± l (cid:48) lνm ( d ) , (25) A h ± l (cid:48) lνm ( d ) = ( − m i l − l (cid:48) ± ν (2 ν + 1) (cid:115) (2 l + 1)(2 l (cid:48) + 1) l ( l + 1) l (cid:48) ( l (cid:48) + 1) (cid:18) l l (cid:48) ν (cid:19) (cid:18) l l (cid:48) νm − m (cid:19) h ν ( dk o ) , (26) A j ± l (cid:48) lνm ( d ) = ( − m i l − l (cid:48) ± ν (2 ν + 1) (cid:115) (2 l + 1)(2 l (cid:48) + 1) l ( l + 1) l (cid:48) ( l (cid:48) + 1) (cid:18) l l (cid:48) ν (cid:19) (cid:18) l l (cid:48) νm − m (cid:19) j ν ( dk o ) . (27)In Eqs. (24)–(27), the large round brackets are Wigner-3j symbols, j ν and h ν denote the spherical Bessel andHankel functions of the first kind respectively, and theplus and minus signs indicate the sign of the necessarycoordinate transformation. Noting that the field of a ˆ z -polarized dipole radiation consists entirely of the N , outgoing wave, E z ( r, θ, φ ) = ik o √ π N , ( r, θ, φ ) , (28)the field of an off-axis dipole can thus be computed fromEqs. (24)–(27) by setting p = N, l = 1 , and m = 0implying, through the Wigner selection rules that theon-origin regular M waves, Rg M , waves need not be in-cluded. Using this knowledge, the expansion coefficientsfor an off-origin, on-axis, dipole can be reduce to E z − d ˆ z ( r ) = ik o √ π N , ( r + d ˆ z )= (cid:40) ik o √ π (cid:80) ∞ l =1 U in − NN,l, , Rg N l, ( r ) ( r < d ) , ik o √ π (cid:80) ∞ l =1 U out − NN,l, , N l, ( r ) ( r > d ) , (29) E z + d ˆ z ( r ) = ik o √ π N , ( r − d ˆ z )= (cid:40) ik o √ π (cid:80) ∞ l =1 U in+ NN,l, , Rg N l, ( r ) ( r < d ) , ik o √ π (cid:80) ∞ l =1 U out+ NN,l, , N l, ( r ) ( r > d ) , (30)where U in − NN,l, , = l +1 (cid:88) ν = l − L νl (cid:18) l ν (cid:19) h ν ( dk o ) , (31) U in+ NN,l, , = l +1 (cid:88) ν = l − ( − ν L νl (cid:18) l ν (cid:19) h ν ( dk o ) . (32) U out − NN,l, , and U out+ NN,l, , are equivalently defined but withthe Hankel functions of the first kind, h ν ( dk o ), replacedby Bessel functions of the first kind. In these expressions, summation limits are set by explicit evaluation of theWigner-3j selection rules and L νl = (2 ν + 1) 2 + l ( l + 1) − ν ( ν + 1)2 i − ν − l (cid:115) l + 32 l ( l + 1) . (33) B. Computation for two-dimensional examples
The Green’s function representations used for com-puting bounds in all two-dimensional examples were ob-tained by the open-source ceviche
FDFD package [106].The basis used for representing fields in all such cases arethe individual discretization pixels of the FDFD compu-tation. Correspondingly, the matrix representation of theGreen’s function connects every pixel with a numeric ap-proximation of the field generated by a dipole source atthe location of the pixel.For all studies with a free space background, the trans-lational invariance of the Maxwell is exploited to con-struct the matrix representation of the Green’s functionfrom a single dipole field solve. Concretely, consider a do-main of shape ( N x , N y ) in pixels. The required represen-tation of G is then a square matrix of dimension N x N y ,where the j -th column is the field created by an electricdipole at position j . This matrix can be obtained in asingle solve by placing a dipole at the center of a larger(2 N x − , N y −
1) domain and observing the calculatedfield in sliding a window of size ( N x , N y ).To account for the cross-source constraints in the pres-ence of multiple sources, compared to the approaches wehave used in prior works [45], the fields and matrices foreach of the individual sources are grouped into “super”vectors and matrices: (cid:12)(cid:12) T (cid:11) = (cid:12)(cid:12) T (cid:11) ... (cid:12)(cid:12) T N (cid:11) Z TT = Z T T . . . Z T TN ... . . . ... Z T N T . . . Z T N T N (34)As such, the dimension of the vectors and matrices in-volved in the calculation of a N -source bound scales lin-early with N . To facilitate computation of bounds for3large device footprints and values of N , an Arnoldi basisfor U is computed to more efficiently represent the | T (cid:105) image field. For any given problem, the only vectors thatinteract with the primal degrees of freedom, (cid:8)(cid:12)(cid:12) T n (cid:11)(cid:9) , arethe sources, (cid:8)(cid:12)(cid:12) E in (cid:11)(cid:9) , and the linear coefficients of (cid:8)(cid:12)(cid:12) T n (cid:11)(cid:9) in the primal objective (cid:8)(cid:12)(cid:12) O lin n (cid:11)(cid:9) . Hence, these vectors B = (cid:2)(cid:12)(cid:12) E i (cid:11) . . . (cid:12)(cid:12) E i N (cid:11) (cid:12)(cid:12) O lin1 (cid:11) . . . (cid:12)(cid:12) O lin1 (cid:11)(cid:3) , (35)posses favorable convergence characteristic for computingthe dual via Arnoldi iterations. For communications typeproblems (cid:12)(cid:12) O lin n (cid:11) = G † od (cid:0)(cid:12)(cid:12) E (cid:5) n (cid:11) − (cid:12)(cid:12) E in (cid:11)(cid:1) , while for the met-alens problem (cid:12)(cid:12) O lin n (cid:11) = G † od (cid:12)(cid:12) E in (cid:11) . To begin the Arnoldiprocedure, B is orthonormalized to give ¯ B . The i -thiteration is then computed by generating B i = U ¯ B i − ,and then orthonormalizing B i both internally and with respect to ¯ B , . . . , ¯ B i − , giving ¯ B i . The partial basis forthe generated block Krylov subspace at the end of the i -th iteration is the thus the aggregate of all the columnvectors of ¯ B , . . . , ¯ B i . Denoting U restricted to the i -thKrylov subspace as U i , the representation is effectivelycomplete when the column norms of U − i ¯ B converge towithin a certain tolerance, Ref. [29].Inverse design results for all two dimensional examplesare also computed using ceviche using a standard topol-ogy optimization [91]. For sake of comparison with theoptimization problem underlying the bound formulation,grayscale structures are treated as viable. That is, in anyparticular design, the permittivity value of a given pixelcan take any value in the convex set bounded by 0 andthe quoted material value. [1] Sergio Verdu. Fifty years of Shannon theory. IEEETransactions on information theory , 44(6):2057–2078,1998.[2] Pierre-Marie Robitaille. Kirchhoff’s law of thermalemission: 150 years.
Progress in Physics , 4:3–13, 2009.[3] Claude E Shannon. A mathematical theory of com-munication.
The Bell System Technical Journal , 27(3):379–423, 1948.[4] John Brian Pendry. Quantum limits to the flow of in-formation and entropy.
Journal of Physics A: Mathe-matical and General , 16(10):2161, 1983.[5] Carlton M Caves and Peter D Drummond. Quantumlimits on bosonic communication rates.
Reviews of Mod-ern Physics , 66(2):481, 1994.[6] Micha(cid:32)l Horodecki, Pawe(cid:32)l Horodecki, and RyszardHorodecki. Limits for entanglement measures.
Phys-ical Review Letters , 84(9):2014, 2000.[7] Stefano Pirandola, Riccardo Laurenza, Carlo Ottaviani,and Leonardo Banchi. Fundamental limits of repeater-less quantum communications.
Nature communications ,8(1):1–15, 2017.[8] David AB Miller. Waves, modes, communications, andoptics: a tutorial.
Advances in Optics and Photonics ,11(3):679–825, 2019.[9] Eli Yablonovitch. Statistical ray optics.
Journal of theOptical Society of America , 72(7):899–907, 1982.[10] Dennis M Callahan, Jeremy N Munday, and Harry AAtwater. Solar cell light trapping beyond the ray opticlimit.
Nano letters , 12(1):214–218, 2012.[11] Andrey E Miroshnichenko and Michael I Tribelsky. Ul-timate absorption in light scattering by a finite obstacle.
Physical Review Letters , 120(3):033902, 2018.[12] Hila Hashemi, Baile Zhang, John D Joannopoulos, andSteven G Johnson. Delay-bandwidth and delay-loss lim-itations for cloaking of large objects.
Physical ReviewLetters , 104(25):253903, 2010.[13] Francesco Monticone and Andrea Al`u. Invisibility ex-posed: physical bounds on passive cloaking.
Optica , 3(7):718–724, 2016.[14] Maxence Cassier and Graeme W Milton. Bounds onherglotz functions and fundamental limits of broadband passive quasistatic cloaking.
Journal of MathematicalPhysics , 58(7):071504, 2017.[15] H. John Caulfield and Tomas Hirschfeld. Optical com-munication at the source bandwidth limit.
Applied Op-tics , 16(5):1184–1186, 1977.[16] Kosmas L Tsakmakidis, Linfang Shen, Sebastian ASchulz, Xiaodong Zheng, John Upham, Xiaohua Deng,Hatice Altug, Alexander F Vakakis, and Robert WBoyd. Breaking lorentz reciprocity to overcome thetime-bandwidth limit in physics and engineering.
Sci-ence , 356(6344):1260–1264, 2017.[17] Sander A Mann, Dimitrios L Sounas, and Andrea Al`u.Nonreciprocal cavities and the time–bandwidth limit.
Optica , 6(1):104–110, 2019.[18] Yaser S Abu-Mostafa and Demetri Psaltis. Optical neu-ral computers.
Scientific American , 256(3):88–95, 1987.[19] Tyler W Hughes, Ian AD Williamson, Momchil Minkov,and Shanhui Fan. Wave physics as an analog recurrentneural network.
Science advances , 5(12):eaay6946, 2019.[20] Ying Zuo, Bohan Li, Yujun Zhao, Yue Jiang, You-Chiuan Chen, Peng Chen, Gyu-Boong Jo, Junwei Liu,and Shengwang Du. All-optical neural network withnonlinear activation functions.
Optica , 6(9):1132–1137,2019.[21] Wim Bogaerts, Daniel P´erez, Jos´e Capmany, David ABMiller, Joyce Poon, Dirk Englund, FrancescoMorichetti, and Andrea Melloni. Programmablephotonic circuits.
Nature , 586(7828):207–216, 2020.[22] David J Richardson, John M Fini, and Lynn E Nel-son. Space-division multiplexing in optical fibres.
Na-ture Photonics , 7(5):354–362, 2013.[23] Juhao Li, Fang Ren, Tao Hu, Zhengbin Li, YongqiHe, Zhangyuan Chen, Qi Mo, and Guifang Li. Re-cent progress in mode-division multiplexed passive op-tical networks with low modal crosstalk.
Optical FiberTechnology , 35:28–36, 2017.[24] Zhenshan Yang, Xiaoguang Zhang, Bin Zhang, XiaZhang, Zhentao Zhang, Xiangguo Meng, and ChenglinBai. Density-matrix formalism for modal coupling anddispersion in mode-division multiplexing communica-tions systems.
Optics Express , 28(13):18658–18680, Science , 363(6433):1333–1338, 2019.[26] Lianlin Li, Ya Shuang, Qian Ma, Haoyang Li, HantingZhao, Menglin Wei, Che Liu, Chenglong Hao, Cheng-Wei Qiu, and Tie Jun Cui. Intelligent metasurface im-ager and recognizer.
Light: Science & Applications , 8(1):1–9, 2019.[27] Hamid Rajabalipanah, Ali Abdolali, Shahid Iqbal, LeiZhang, and Tie Jun Cui. How do space-time digitalmetasurfaces serve to perform analog signal processing? arXiv:2002.06773 , 2020.[28] Vittorio Giovannetti, Seth Lloyd, Lorenzo Maccone,and Jeffrey H Shapiro. Sub-rayleigh-diffraction-boundquantum imaging.
Physical Review A , 79(1):013827,2009.[29] Sean Molesky, Pengning Chao, Weiliang Jin, and Ale-jandro W Rodriguez. Global T -operator bounds on elec-tromagnetic scattering: Upper bounds on far-field crosssections. Physical Review Research , 2(3):033172, 2020.[30] Owen D Miller, Steven G Johnson, and Alejandro WRodriguez. Shape-independent limits to near-field ra-diative heat transfer.
Physical Review Letters , 115(20):204302, 2015.[31] Weiliang Jin, Sean Molesky, Zin Lin, and Alejandro WRodriguez. Material scaling and frequency-selective en-hancement of near-field radiative heat transfer for lossymetals in two dimensions via inverse design.
PhysicalReview B , 99(4):041403(R), 2019.[32] Prashanth S Venkataram, Sean Molesky, Weiliang Jin,and Alejandro W Rodriguez. Fundamental limits to ra-diative heat transfer: the limited role of nanostructur-ing in the near-field.
Physical Review Letters , 124(1):013904, 2020.[33] Alexander Y Piggott, Jesse Lu, Thomas M Babinec,Konstantinos G Lagoudakis, Jan Petykiewicz, and Je-lena Vuˇckovi´c. Inverse design and implementation ofa wavelength demultiplexing grating coupler.
ScientificReports , 4(1):1–5, 2014.[34] Andrew Michaels and Eli Yablonovitch. Inverse designof near unity efficiency perfectly vertical grating cou-plers.
Optics express , 26(4):4766–4779, 2018.[35] Jianji Yang, David Sell, and Jonathan A Fan. Freeformmetagratings based on complex light scattering dynam-ics for extreme, high efficiency beam steering.
Annalender Physik , 530(1):1700302, 2018.[36] Prachi Thureja, Ghazaleh Kafaie Shirmanesh, Kather-ine T Fountaine, Ruzan Sokhoyan, Meir Grajower, andHarry A Atwater. Array-level inverse design of beamsteering active metasurfaces.
ACS Nano , 2020.[37] Oguzhan Akgol, Emin Unal, Olcay Altintas, MuharremKaraaslan, Faruk Karadag, and Cumali Sabah. Designof metasurface polarization converter from linearly po-larized signal to circularly polarized signal.
Optik , 161:12–19, 2018.[38] Haejun Chung and Owen D Miller. High-na achromaticmetalenses by inverse design.
Optics Express , 28(5):6945–6965, 2020.[39] Jesse Lu and Jelena Vuˇckovi´c. Inverse design ofnanophotonic structures using complementary convexoptimization.
Optics Express , 18(4):3793–3804, 2010.[40] Jesse Lu and Jelena Vuˇckovi´c. Objective-first designof high-efficiency, small-footprint couplers between ar- bitrary nanophotonic waveguide modes.
Optics express ,20(7):7221–7236, 2012.[41] Owen D Miller, Athanasios G Polimeridis, M. T. HomerReid, Chia Wei Hsu, Brendan G DeLacy, John DJoannopoulos, Marin Soljaˇci´c, and Steven G Johnson.Fundamental limits to optical response in absorptivesystems.
Optics Express , 24(4):3329–3364, 2016.[42] Guillermo Angeris, Jelena Vuˇckovi´c, and Stephen PBoyd. Computational bounds for photonic de-sign.
ACS Photonics , 6(5):1232, 2019. doi:10.1021/acsphotonics.9b00154.[43] Sean Molesky, Weiliang Jin, Prashanth S Venkataram,and Alejandro W Rodriguez. T -operator bounds onangle-integrated absorption and thermal radiation forarbitrary objects. Physical Review Letters , 123:257401,2019.[44] Sean Molesky, Prashanth S Venkataram, Weiliang Jin,and Alejandro W Rodriguez. Fundamental limits toradiative heat transfer: theory.
Physical Review B , 101(3):035408, 2020.[45] Sean Molesky, Pengning Chao, and Alejandro W Ro-driguez. Hierarchical mean-field T -operator bounds onelectromagnetic scattering: Upper bounds on near-fieldradiative purcell enhancement. Physical Review Re-search , 2(4):043398, 2020.[46] Mats Gustafsson, Kurt Schab, Lukas Jelinek, andMiloslav Capek. Upper bounds on absorption and scat-tering.
New Journal of Physics , 2020.[47] Prashanth S Venkataram, Sean Molesky, PengningChao, and Alejandro W Rodriguez. Fundamental limitsto attractive and repulsive casimir-polder forces.
Phys-ical Review A , 101(5):052115, 2020.[48] Zeyu Kuang, Lang Zhang, and Owen D Miller.Maximal single-frequency electromagnetic response. arXiv:2002.00521 , 2020.[49] Kurt Schab, Austin Rothschild, Kristi Nguyen, MiloslavCapek, Lukas Jelinek, and Mats Gustafsson. Trade-offs in absorption and scattering by nanophotonic struc-tures.
Optics Express , 28(24):36584–36599, 2020.[50] Rahul Trivedi, Guillermo Angeris, Logan Su, StephenBoyd, Shanhui Fan, and Jelena Vuckovic. Fundamentalbounds for scattering from absorptionless electromag-netic structures. arXiv:2003.00374 , 2020.[51] Zeyu Kuang and Owen D Miller. Computational boundsto light-matter interactions via local conservation laws. arXiv preprint arXiv:2008.13325 , 2020.[52] L Jelinek, M Gustafsson, M Capek, and K Schab. Sub-structure limits to optical phenomena. In , pages 400–402. IEEE.[53] Guillermo Angeris, Jelena Vuˇckovi´c, and Stephen Boyd.Heuristic methods and performance bounds for photonicdesign.
Optics Express , 29(2):2827–2854, 2021.[54] Walter Rudin.
Real and complex analysis . TataMcGraw-hill education, 2006.[55] David AB Miller. All linear optical devices are modeconverters.
Optics express , 20(21):23985–23993, 2012.[56] David AB Miller. Self-configuring universal linear opti-cal component.
Photonics Research , 1(1):1–15, 2013.[57] David AB Miller, Linxiao Zhu, and Shanhui Fan. Uni-versal modal radiation laws for all thermal emitters.
Proceedings of the National Academy of Sciences , 114(17):4336–4341, 2017. [58] Sunil Pai, Ben Bartlett, Olav Solgaard, and David ABMiller. Matrix optimization on universal unitary pho-tonic devices. Physical Review Applied , 11(6):064044,2019.[59] David A B Miller. Fundamental limit for optical com-ponents.
Journal of the Optical Society of America B ,24(10):A1–A18, 2007.[60] David AB Miller. How complicated must an opticalcomponent be?
JOSA A , 30(2):238–251, 2013.[61] Leung Tsang, Jin Au Kong, and Kung-Hau Ding.
Scat-tering of electromagnetic waves: Theories and applica-tions , volume 27. John Wiley & Sons, 2004.[62] Andreas Kirsch and Armin Lechleiter. The operatorequations of Lippmann–Schwinger type for acoustic andelectromagnetic scattering problems in L2.
ApplicableAnalysis , 88(6):807–830, 2009.[63] Martin Costabel, Eric Darrigrand, and Hamdi Sakly.The essential spectrum of the volume integral operatorin electromagnetic scattering by a homogeneous body.2012.[64] Matthias Kr¨uger, Giuseppe Bimonte, Thorsten Emig,and Mehran Kardar. Trace formulas for nonequilibriumcasimir interactions, heat radiation, and heat transferfor arbitrary objects.
Physical Review B , 86(11):115423,2012.[65] David Colton and Rainer Kress.
Inverse acoustic andelectromagnetic scattering theory , volume 93. SpringerNature, 2019.[66] Note1. All operators and fields are supposed to be fre-quency dependent, and every stated relation holds atany given frequency.[67] Note2. The potential of the scatterer is constructedvia the explicit spatial projection operator I s , which isdefined to be identity at spatial locations within thescatterer and zero at all other points in the domain Ω,and the implicit spatial projection into the scatteringobject encoded in the definition the scattering operator, T s .[68] Note3. The specification of an initial field is essentiallyequivalent to the introduction of a source, or flux bound-ary condition, in Maxwell’s equations.[69] Jingfeng Liu, Ming Zhou, Lei Ying, Xuewen Chen, andZongfu Yu. Enhancing the optical cross section of quan-tum antenna. Physical Review A , 95(1):013814, 2017.[70] A Femius Koenderink. Single-photon nanoantennas.
ACS Photonics , 4(4):710–722, 2017.[71] Kevin C Cox, David H Meyer, Fredrik K Fatemi, andPaul D Kunz. Quantum-limited atomic receiver in theelectrically small regime.
Physical Review Letters , 121(11):110502, 2018.[72] Sudha Mokkapati and KR Catchpole. Nanophotoniclight trapping in solar cells.
Journal of Applied Physics ,112(10):101101, 2012.[73] Xing Sheng, Juejun Hu, Jurgen Michel, and Lionel CKimerling. Light trapping limits in plasmonic solarcells: an analytical investigation.
Optics Express , 20(104):A496–A501, 2012.[74] Vidya Ganapati, Owen D Miller, and Eli Yablonovitch.Light trapping textures designed by electromagnetic op-timization for subwavelength thick solar cells.
IEEEJournal of Photovoltaics , 4(1):175–182, 2013.[75] E Phan-huy Hao. Quadratically constrained quadraticprogramming: Some applications and a method for so-lution.
Zeitschrift f¨ur Operations Research , 26(1):105– 119, 1982.[76] Subhonmesh Bose, Dennice F Gayme, K Mani Chandy,and Steven H Low. Quadratically constrained quadraticprograms on acyclic graphs with application to powerflow.
IEEE Transactions on Control of Network Sys-tems , 2(3):278–287, 2015.[77] Note4. Technically, B should be complete, and henceinfinite dimensional. In practice, however, B must trun-cated to some finite set in order to apply numericalmethods. A similar approximation is needed in any nu-meric solver for electromagnetics.[78] Jaehyun Park and Stephen Boyd. General heuristicsfor nonconvex quadratically constrained quadratic pro-gramming. arXiv:1703.07870 , 2017.[79] Stephen Boyd and Lieven Vandenberghe. Convex opti-mization . Cambridge university press, 2004.[80] Note5. Cluster refinements, as shown in Ref. [45], canbe given a mathematical ordering that is preserved inthe limit values obtained in the associated optimization,whether or not the dual is used to relax the optimiza-tion.[81] Note6. In particular, if Ω c = Ω, so that there isonly a single cluster filling the entire domain, then thesymmetric (Sym [ A ] = (cid:0) A + A † (cid:1) /
2) and antisymmetric(Sym [ A ] = (cid:0) A + A † (cid:1) /
2) parts of the scattering con-straints impose that both real and reactive power mustbe conserved, as was studied in Refs. [43, 46].[82] Timm Plefka. Convergence condition of the tap equa-tion for the infinite-ranged ising spin glass model.
Jour-nal of Physics A , 15(6):1971, 1982.[83] Jiasen Jin, Alberto Biella, Oscar Viyuela, LeonardoMazza, Jonathan Keeling, Rosario Fazio, and DavideRossini. Cluster mean-field approach to the steady-statephase diagram of dissipative spin systems.
Physical Re-view X , 6(3):031011, 2016.[84] Mojtaba Nourian and Peter E Caines. (cid:15) -nash mean fieldgame theory for nonlinear stochastic dynamical systemswith major and minor agents.
SIAM Journal on Controland Optimization , 51(4):3302–3331, 2013.[85] Nevroz S¸en and Peter E Caines. Mean field game theorywith a partially observed major agent.
SIAM Journalon Control and Optimization , 54(6):3174–3224, 2016.[86] Hilbert J Kappen and Francisco de Borja Rodr´ıguez Or-tiz. Boltzmann machine learning using mean field theoryand linear response correction. In
Advances in NeuralInformation Processing Systems , pages 280–286, 1998.[87] Minmin Chen, Jeffrey Pennington, and Samuel SSchoenholz. Dynamical isometry and a mean field the-ory of rnns: Gating enables signal propagation in recur-rent neural networks. arXiv:1806.05394 , 2018.[88] Note7. Duality can also be used to solve a relaxationof the mean-field optimization that in many cases is infact a true solution [79].[89] David Mendlovic and Adolf W Lohmann. Space–bandwidth product adaptation and its application to su-perresolution: fundamentals.
JOSA A , 14(3):558–562,1997.[90] Mark A Neifeld. Information, resolution, and space–bandwidth product.
Optics Letters , 23(18):1477–1479,1998.[91] Rasmus E Christiansen and Ole Sigmund. Inverse de-sign in photonics by topology optimization: tutorial.
JOSA B , 38(2):496–509, 2021. [92] Sergey Kruk and Yuri Kivshar. Functional meta-opticsand nanophotonics governed by mie resonances. ACSPhotonics , 4(11):2638–2649, 2017.[93] Tomer Lewi, Nikita A Butakov, Hayden A Evans,Mark W Knight, Prasad P Iyer, David Higgs, HamidChorsi, Juan Trastoy, Javier Del Valle Granda, IlyaValmianski, et al. Thermally reconfigurable meta-optics.
IEEE Photonics Journal , 11(2):1–16, 2019.[94] Isabelle Staude, Thomas Pertsch, and Yuri S Kivshar.All-dielectric resonant meta-optics lightens up.
ACSPhotonics , 6(4):802–814, 2019.[95] Zin Lin, Charles Roques-Carmes, Rapha¨el Pestourie,Marin Soljaˇci´c, Arka Majumdar, and Steven G John-son. End-to-end inverse design for inverse scatteringvia freeform metastructures. arXiv:2006.09145 , 2020.[96] Philipp-Immanuel Schneider, Nicole Srocka, Sven Rodt,Lin Zschiedrich, Stephan Reitzenstein, and SvenBurger. Numerical optimization of the extraction ef-ficiency of a quantum-dot based single-photon emitterinto a single-mode fiber.
Optics Express , 26(7):8479–8492, 2018.[97] Lian Shen, Xiao Lin, Mikhail Y Shalaginov, Tony Low,Xianmin Zhang, Baile Zhang, and Hongsheng Chen.Broadband enhancement of on-chip single-photon ex-traction via tilted hyperbolic metamaterials.
AppliedPhysics Reviews , 7(2):021403, 2020.[98] Raymond A Wambold, Zhaoning Yu, Yuzhe Xiao, Ben-jamin Bachman, Gabriel Jaffe, Shimon Kolkowitz, Jen-nifer T Choy, Mark A Eriksson, Robert J Hamers,and Mikhail A Kats. Adjoint-optimized nanoscalelight extractor for nitrogen-vacancy centers in diamond.
Nanophotonics , 10(1):393–401, 2021. [99] Srivatsa Chakravarthi, Pengning Chao, Christian Ped-erson, Sean Molesky, Karine Hestroffer, Fariba Hatami,Alejandro W Rodriguez, and Kai-Mei C Fu. Inverse-designed photon extractors for optically addressable de-fect qubits. arXiv:2007.12344 , 2020.[100] Note8. It should be mentioned that we have encoun-tered situations in which very high numeric precisionwas required to verify the presence of strong duality.[101] Jack Valmadre, Luca Bertinetto, Joao Henriques, An-drea Vedaldi, and Philip HS Torr. End-to-end represen-tation learning for correlation filter based tracking. In
Proceedings of the IEEE Conference on Computer Vi-sion and Pattern Recognition , pages 2805–2813, 2017.[102] Lisa A Poyneer and Jean-Pierre V´eran. Optimal modalfourier-transform wavefront control.
JOSA A , 22(8):1515–1526, 2005.[103] Vincent Sitzmann, Steven Diamond, Yifan Peng, XiongDun, Stephen Boyd, Wolfgang Heidrich, Felix Heide,and Gordon Wetzstein. End-to-end optimization ofoptics and image processing for achromatic extendeddepth of field and super-resolution imaging.
ACMTransactions on Graphics (TOG) , 37(4):1–13, 2018.[104] Ronald C Wittmann. Spherical wave operators and thetranslation formulas.
IEEE Transactions on Antennasand Propagation , 36(8):1078–1087, 1988.[105] Yu-lin Xu. Calculation of the addition coefficients inelectromagnetic multisphere-scattering theory.
Journalof Computational Physics , 127(2):285–298, 1996.[106] Tyler W Hughes, Ian AD Williamson, Momchil Minkov,and Shanhui Fan. Forward-mode differentiation ofmaxwell’s equations.