Maximum Shannon Capacity of Photonic Structures
Alessio Amaolo Pengning Chao Benjamin Strekha Stefan Clarke Jewel Mohajan Sean Molesky Alejandro W. Rodriguez


APS/123-QED




alessioamaolo@princeton.edu
Department of Chemistry, Princeton University, Princeton, New Jersey 08544, USA
Equal contribution

pchao827@mit.edu
Department of Mathematics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Equal contribution

Department of Electrical and Computer Engineering, Princeton University, Princeton, New Jersey 08544, USA

Department of Operations Research and Financial Engineering, Princeton University, Princeton, New Jersey 08544, USA

Department of Electrical and Computer Engineering, Princeton University, Princeton, New Jersey 08544, USA

Department of Engineering Physics, Polytechnique Montréal, Montréal, Québec H3T 1J4, Canada

Department of Electrical and Computer Engineering, Princeton University, Princeton, New Jersey 08544, USA


§ ABSTRACT 
Information transfer through electromagnetic waves is an important problem that touches a variety of technologically relevant application areas, including computing, telecommunications, and power management. 
Prior attempts to establish limits on optical information transfer have focused exclusively on waves propagating through vacuum or a known photonic structure (prescribed wave sources and receivers).
In this article, we describe a mathematical theory that addresses fundamental questions relating to optimal information transfer in photonic devices. 
Combining the mathematics of information theory, wave scattering, and optimization theory, we derive a method to compute upper bounds on the maximum Shannon capacity that may be achieved by structuring senders, receivers, and their surrounding environment.
This approach provides a means to understand how material selection, device size, and general geometrical features impact power allocation, communication channels, and bit-rate in photonics.  
Arbitrary structuring provides a means to optimize the power requirements and channel capacities for communication and leads to a non-convex problem that is significantly more difficult to bound than its fixed structure counterpart, which is convex and thus satisfies a known “water-filling” solution. 
We derive a geometry-agnostic convex relaxation of the problem that elucidates fundamental physics and scaling behavior of Shannon capacity with respect to device parameters, as well as the importance of structuring for enhancing capacity.
We also show that in regimes where communication is dominated by power insertion requirements, maximum Shannon capacity maps to a biconvex optimization problem in the basis of singular vectors of the Green's function. 
This problem admits analytical solutions that give physically intuitive interpretations of channel and power allocation and reveals how Shannon capacity varies with signal-to-noise ratio.
Proof of concept numerical examples show that bounds are within an order of magnitude of achievable device performance in a variety of settings, with even better agreement in high signal-to-noise regimes, and successfully predict the scaling of performance with channel noise.
The presented methodologies have implications for the optimization of antennas, integrated photonic devices, metasurface kernels, MIMO space-division multiplexers, and waveguides to maximize communication efficiency and bit-rates.


Maximum Shannon Capacity of Photonic Structures
    Alejandro W. Rodriguez
    September 9, 2024
===============================================





§ INTRODUCTION


From the 19th century electrical telegraph to modern 5G telecommunications, electromagnetics has played an increasingly important role in human communication. 
Paralleling this rapid improvement in physical technologies, our understanding of the fundamental meaning of information has also steadily improved.
Building upon earlier work by Nyquist <cit.> and Hartley <cit.>, Shannon demonstrated in two landmark papers <cit.> that not only is it possible to transmit information at a finite rate with arbitrarily small error through noisy channels, but that the maximum transmission rate—now known as the Shannon capacity—can also be explicitly calculated as a function of the channel bandwidth and signal-to-noise ratio.  
Shannon's original derivation focused on scalar, time-varying signals subject to additive white Gaussian noise, but the principles he elucidated are universal and can be applied to more general physical settings, including electromagnetic wave phenomena <cit.>. 
Information capacity results have been obtained for a variety of systems including fiber optics <cit.> and wireless communications networks <cit.>; in particular, a major area of research are multiple-input multiple-output (MIMO) systems that make use of spatial multiplexing, wherein the signal fields have spatial as well as temporal degrees of freedom <cit.>. 
Communication in the presence of a scatterer has also been analyzed by choosing an appropriate optical characteristic proxy, e.g., the number of waves escaping a region <cit.> or the sum of field amplitudes at the receiver <cit.>.
These figures of merit are intimately related to other quantities of power transfer such as thermal radiation <cit.> and near-field radiative heat transfer <cit.> and establish intuitive connections between fields at the receiver and communication, but do not directly maximize Shannon capacity. 

These prior results generally work well if the wave propagation medium is fixed, be it antenna design in free-space <cit.> or guided modes in a fiber <cit.>. 
Yet, they do not capture the range of physics and possibilities of control and optical response achievable via (nano)photonic structuring. 
Antennas at either the sender or receiver may be designed to enhance gain and directivity (therefore bit-rate), and similarly their environment may be structured to achieve greater communication and field enhancements (e.g., parabolic reflectors, metasurfaces, waveguides, fibers, etc.).
In fact, there is emerging consensus on the increasing importance of complex structured environments for electromagnetic communication in a wide range of contexts, including silicon photonics <cit.>, 
on-chip optical interconnects <cit.>, and integrated image processing <cit.>.
Naturally, for each of these systems, the question of how much environmental structuring can improve their information transfer capabilities is key to quantifying future potential.
Recently, significant progress has been made in computational photonic design <cit.> by utilizing techniques from mathematical optimization to maximize field objectives over many degrees of freedom. 
In particular, these designs' ability to concentrate power into selective wavelengths (e.g., resonant modes) and manipulate fields in counterintuitive ways has a wide array of applications and implications for photonic templates of communication: metasurface kernels <cit.>, integrated resonators <cit.>, and channel engineering <cit.>. 

In this work, we leverage recent developments in photonic optimization to present a unified framework for investigating the extent to which structured photonic devices may enhance information transfer. First, we show how the structure-dependent, non-convex problem of maximizing information transfer can be relaxed into a shape-independent, convex problem to arrive at a bound on Shannon capacity subject to drive-power and physical wave constraints derived from Maxwell's equations. 
Second, we show that in regimes where power requirements are dominated by insertion impedance (i.e., internal resistance or contact impedance in the driving currents), the maximization of the Shannon capacity can be significantly simplified in the basis of singular vectors of the Green's function, written in terms of a single signal-to-noise free parameter, and further relaxed to a biconvex problem that admits analytical solutions.  
We present several illustrative examples where we compare these bounds to inverse designs and show that not only do the bounds predict trends but are within an order of magnitude of achievable performance for a wide range of signal-to-noise ratios. 

Bounds on the Shannon capacity provide performance targets for communication and elucidate how information capacity scales with key device parameters (e.g., material choices, device sizes, or bandwidth). 
Even infinitely large spaces are guaranteed to have decaying singular values by the compactness of the Green's function <cit.>, meaning that only a finite number of channels will be useful in transferring information.
Additionally, the fundamental limitations encoded in Maxwell's equations (spectral sum rules, finite speed of light, space-bandwdith products, etc.) dictate a finite limit on the ability to concentrate power over selective wavelengths (e.g., via resonances) or propagate fields over a given distance <cit.>. The proposed method paves the way for understanding how such trade-offs affect electromagnetic information transfer.



§ THEORY
 



In this section, we begin with a brief review of several defining relations of information capacity that pertain to wave communication. We show how the standard formula for MIMO communication can be derived from the fundamental definition of the Shannon capacity, subject to a constraint on the magnitude of the input signal. 
This constraint manifests electromagnetically as a constraint on input current amplitudes; we show how we can also derive a constraint with respect to drive power, connecting to prior work on vacuum information transfer in the context of antenna design <cit.>. 
By incorporating the freedom to structure senders, receivers, and their surrounding environments, we formulate a non-convex structural optimization problem that bounds the maximum Shannon capacity over any photonic structure. 
We show how the problem can be relaxed to a convex optimization problem or, alternatively, to a biconvex problem on the channel capacities as given by the electromagnetic Green's function.

We consider a MIMO system with input signal x∈ℂ^n and output signal y∈ℂ^m related by 

    =  +

where the ∈ℂ^m × n is the channel matrix and the channel has additive zero-mean complex Gaussian noise with zero mean and covariance (^†) = N.
The channel matrix and noise prescription relate “input" signals x to “output" signals y, which are viewed as random variables given the randomness of noise and the probabilistic argument used by Shannon to derive the Shannon capacity <cit.>. Communication is possible when x and y are dependent: an observation of y changes the a posteriori probability for any given x. Denote the joint probability distribution of x and y as γ( x,  y); the marginal and conditional probability distributions are given by 

    
    γ_ x( x)    = ∫γ( x,  y)  y 
    γ( x| y)    = γ( x,  y) / γ_ y( y)


and symmetrically in .
The differential and conditional differential entropy of  can now be written as

    ()    = ∫γ_ x() log_2(1/γ_x()) , and
    (|)    = ∫γ_ y() ∫γ(|) log_2(1/γ(|)) .

() is a continuous analog of the discrete information entropy which encodes the level of uncertainty in x 

and is a measure of the information gained upon observation of x; (|) in turn is the expected entropy of  conditioned on observation of  <cit.>. 
Their difference is the mutual information between  and , denoted ( ; ): 

    (;) = () - (|) = () - (|).

Intuitively, (;) is the expected reduction in uncertainty (and hence gain in information) of x upon observation of . The second equality follows from Eq. (<ref>) and indicates that (;) is symmetric between  and ; note that for our MIMO channel (<ref>) the second term is simply ( | ) = ().
Shannon's noisy coding theorem states that a tight upper limit on the error-free information transfer rate across the channel achievable via encoding is the Shannon capacity C = sup_γ_x (;) <cit.>. 
Over a channel with continuous signals and finite additive noise as in (<ref>), C diverges without further constraints on the magnitude of : 
intuitively, a larger signal magnitude will make the corruption of additive noise n proportionally smaller. Without further details of the underlying physics of the channel, this is addressed by enforcing a cap on the average magnitude of , which is equivalent to a constraint on the covariance Q of γ_x:

    _γ_x[^†] = (_γ_x[^†]) ≡ Q ≤ P.


The Shannon capacity of our MIMO system is thus

    C(, P) = sup_γ_x, ( Q)≤ P  (;)

where the supremum is taken over all possible probability distributions γ_x() for .
Given (<ref>), one can show that among all γ_x with a particular covariance Q, (;) is maximized by the corresponding Gaussian distribution with covariance Q <cit.>. 
Thus when solving (<ref>) we can only consider Gaussian distributions, leading to the well-known expression for the Shannon capacity per unit bandwidth C for a MIMO  system <cit.>
 
    C(,P)    = max_Q≻0, (Q)≤ P  () - () 
       = max_Q≻0, (Q)≤ P  log_2 (Q^† + N) 
       - log_2 (N) 
       = max_Q≻0, (Q)≤ P  log_2 ( + 1/NQ^†).
 
This is a convex optimization problem and has a analytical solution that is often referred to as “water-filling" in the literature. For further mathematical details the reader is referred to <cit.>.

We will now consider Shannon capacity in the physical context of information transfer via electromagnetic fields and describe how to incorporate the possibility of photonic structuring into bounds on the information transfer rate. 
The schematic setup is illustrated in Fig. <ref>, where the input  is a free current distribution _i within a sender region S, and the output  is the electric field _R generated in a receiver region R. 
Working in the frequency domain with dimensionless units of μ_0=ϵ_0=1 and vacuum impedance Z = √(μ_0 / ϵ_0) = 1, the electric field e in a region obeys Maxwell's wave equation ( ∇×∇× - ϵ( r) ω^2 ) ( r; ω) = iω_i( r), where ω is the angular frequency of the source current _i, k= ω, and ϵ( r) is the permittivity distribution in space. 
The total Green's function _t(r,r', ϵ(·); ω) maps input currents j_i to electric fields  via ( r) = i/ω∫_t( r, r') _i( r')  r' and satisfies ( ∇×∇× - ϵ( r) ω^2 )_t( r,  r',ϵ(·);ω) = ω^2 Iδ(r - r'), where I is the unit dyad. 
We also denote χ_D a constant susceptibility factor associated with a photonic device and related to the permittivity profile via ϵ( r) -  = P_ϵ( r) χ_D where P_ϵ( r) is a projection operator containing the spatial dependence of ϵ( r).
In vacuum, there is no photonic structuring so the vacuum Green's function is given by _0 ≡_t( r,  r', ; ω). We denote _t,RS as the sub-block of the Green's function that maps currents in the source region to fields in the receiver region, so the channel matrix H = i/ω_t,RS <cit.>. 
As shown in Fig. <ref>, we can modify _t and H through the engineering of a photonic structure ϵ( r) within the design region D, which may also incorporate either the S or R regions; the central question we seek to address is the extent to which photonic engineering may boost the Shannon capacity C. 

Some constraints on the input currents _i are needed for C to be finite; instead of the magnitude constraint prescribed in the more abstract formulation of (<ref>), we instead consider a budget P on the total drive power needed to maintain _i: 

    [-1/2i/ω_i^†_t,SS_i] + [ α |_i|^2 ] ≤ P.
 
The total drive power in (<ref>) is comprised of two parts: the first term is the power extracted from _i by the electric field, which includes radiated power and material absorption: optical loss processes that are modeled explicitly via Maxwell's equations. 
The second term is an “insertion cost” for the (abstract) drive producing the input currents, which we assume to be proportional to the squared norm of _i with proportionality factor α.
Depending on the specific drive mechanism, this term may represent, for example, the internal resistance of a power source or contact impedance leading into an antenna <cit.>. 


Adapting the Shannon capacity expression (<ref>) to the specifics of our electromagnetic channel along with the drive power budget (<ref>), keeping in mind that we are also designing the propagation environment (including sender and receiver) via ϵ( r), we arrive at the following structural optimization problem


    max_ϵ(r)     C(ϵ, P) = max_γ__i (_i;_R;ϵ) 
    s.t.     _γ__i[-1/2i/ω_i^†_t,SS(ϵ) _i  + α|_i|^2] ≤ P.



Given that (<ref>) can be written as a linear trace constraint on the covariance of γ__i, we only need to consider Gaussian distributions for the inner optimization, allowing us to rewrite (<ref>) more explicitly as 
 
    
    max_ϵ(r), _i     log_2 ( + PNω^2_t, RS(ϵ) _i _t, RS^† (ϵ)) 
    s.t.     (( α_S + 12ω_t, SS(ϵ)) _i) ≤ 1 
       _i ≽0
        ( - ϵ(r)ω^2) _t(r,r';ϵ) = ω^2 δ(r-r')


where _i is the covariance matrix [_i _i^†]/P and A≡1/2i (A - A^†).
This scaling by P allows the expression of the objective in terms of a familiar signal-to-noise ratio P/N. The two free parameters in Eq. (<ref>) are P/N and α. 
We note that α formally has units of resistivity: in this article, we take these units to be Z λ. Similarly, P/N has units of ϵ_0 c λ^2 in 3D and ϵ_0 c λ in 2D.
Note that the optimization degrees of freedom ϵ( r) may be restricted to any sub-region of the domain (e.g., one may consider a prescribed set of antennas and instead seek only to optimize the region between them).
This problem is non-convex owing to the nonlinear dependence of  on ϵ( r). 
In particular, the water-filling solution does not apply: a bound on Shannon capacity is given by the co-optimization of current covariance and the photonic structure. 


In the following sections, we will derive bounds to (<ref>) via two different methods and discuss our results in the context of prior work. In the first method of Sec. <ref>, we will reinterpret (<ref>) from a maximization over ϵ( r) and γ__i to a maximization over joint probability distributions of _i and the induced polarization current _s = ( ϵ( r) - ) _t(ϵ( r)) _i. 
By relaxing physical consistency requirements on this joint probability distribution, we obtain a convex problem that can be solved with methods in semidefinite programming (SDP) or via a generalized waterfilling approach described in Appendix <ref>. 
Second, in Sec. <ref>, we consider a regime in which the drive power (<ref>) is dominated by insertion losses in the sender, and utilize a convenient basis transformation to express the Shannon capacity in terms of the singular values of _t,RS, yielding an intuitive decomposition into orthogonal channels. 
By relaxing Maxwell's equations to encode wave behavior directly into structurally agnostic bounds on these singular values, we derive a biconvex problem that for a limited number of constraints admits analytical solutions and physically intuitive interpretations of channel and power allocation. 
Both relaxations can in principle be applied in the appropriate regimes of validity, but theoretically understanding their differences is challenging: while both encode the physics of Maxwell's equations, the first does so via power-conservation constraints on polarization currents while the second through the dependence of the channel matrix on its singular values. 
Surprisingly, however, we find that the predicted bounds are close in regimes where they can be compared. 



 §.§ Convex Relaxation via Expected Power Conservation over Joint Distributions
 


For notational convenience we define the combined source-induced current = [ j_i; j_s ] and a function that maps input currents to induced polarization currents S(_i; ϵ) = (ϵ-) _t(ϵ) _i. 
The set of physically feasible joint distributions is then given by

    𝒫_physical = {γ_ | ∃ϵ ∀_i   γ_(_s|_i) = δ(_s - S(_i; ϵ)) }.

Sampling any joint distribution γ_∈𝒫_physical gives  such that _s is the induced polarization current in some ϵ(r) given _i; furthermore for a fixed γ_ this underlying ϵ(r) is consistent across all _i.  
Given _s, the output field _R and the extracted drive power P can be expressed solely using the vacuum Green's function

    _R = i/ω_RS(ϵ) _i = i/ω(_0,RS_i + _0,RD_s), 
    
    P = -1/2i/ω_i^† (_0,SS_i + _0,SD_s) + α |_i|^2.
 
Observe that the joint distribution optimization
      
    
    max_γ_∈𝒫_physical (;_R)  
       s.t. _γ_[-1/2i/ω_i^† (_0,SS_i + _0,SD_s)  + α |_i|^2 ] ≤ P


is exactly equivalent to (<ref>): here, the (ϵ( r), γ__i) pairs have a one-to-one correspondence with γ_∈𝒫_physical, both producing the same output distribution γ__R and hence, the same ℋ(_R) and ℐ. 

So far, (<ref>) is just a rewriting of (<ref>) with the complexity of the problem buried in the description of the set 𝒫_physical; we now relax (<ref>) by considering a superset of 𝒫_physical that has a simpler mathematical structure. 
Instead of insisting that there exists some specific ϵ(r) for a member joint distribution as in (<ref>), we simply require that on average, _i and _s satisfy structure agnostic local energy conservation constraints (see Appendix <ref>):

    _γ_[_i^†_0,DS^†_k _s - _s^†_k _s] = 0

where ≡( χ_D^-† - _0,DD^†), _k is the indicator function over some design sub-region D_k ⊂ D and χ_D is the constant susceptibility of the design region that does not depend on geometry. 
This relaxed requirement gives us the relaxed distribution set

    𝒫_relaxed = {γ_ | _γ_[_i^†_0,DS^†_k _s - _s^†_k _s] = 0}

with an associated Shannon capacity problem

    max_γ_∈𝒫_relaxed (;_R) 
       s.t. _γ_[-1/2i/ω_i^† (_0,SS_i + _0,SD_s) + α |_i|^2] ≤ P.

Given that 𝒫_physical⊂𝒫_relaxed, the optimum of (<ref>) bounds (<ref>) and hence (<ref>) from above. Furthermore, the restrictions of (<ref>) are just linear constraints on the covariance of γ_; we can again restrict attention to just Gaussian distributions to solve (<ref>), yielding

    
    max_     log_2 ( + P/Nω^2[_0,RS   _0,RD] [_0,RS^†
    _0,RD^†] ) 
    s.t.     [α_S + 1/2ω(_0,SS)    1/2ω1/2i_0,SD
    1/2ω-1/2i_0,SD^†   0] ≤ 1 
       P [0    - _0,DS^†_k 
     
            - _k _0,DS   _k]  = 0,  ∀_k  
       ≽0

    

This convex optimization problem over the joint covariance ≡[^†]/P is guaranteed to give a finite bound on the Shannon capacity given finite α, χ_D > 0 (Appendix <ref>) and can be solved using methods in semidefinite programming. 
Although (<ref>) is not strictly an SDP due to the log det objective, -log((·)) is the standard logarithmic barrier function for SDP interior point methods and certain solvers support this objective natively (see for example SDPT3 <cit.>). In our numerical testing, we have found that using existing general convex optimization solvers for (<ref>) is computationally expensive for even wavelength-scale domains, and instead we solve (<ref>) using a “dual waterfilling” approach that is a generalization of the method used in <cit.> (see Appendix <ref>).  



 §.§ Insertion-Dominated Power Requirements and Singular-Value Decompositions
 

Although the formulation above rigorously includes drive-power constraints on the Shannon capacity to all orders in scattering by the photonic structure, it is both practically and pedagogically useful to consider a relaxation of Eq. (<ref>) that assumes drive-power considerations are dominated by insertion losses. 
This yields a solvable biconvex problem that elucidates fundamental principles of power transfer and channel allocation and may prove accurate in cases where structuring (local density of states enhancement) does not significantly influence extracted power, e.g., a laser incident on a scatterer or communication in an electrical wire. 
Mathematically, the contributions on the power constraint given by _t,SS are negligible compared to the insertion losses given by ∼α_i^2, which results in a current amplitude constraint of the form _i^2 ≤ P / α.

We can simplify Eq. (<ref>) by decomposing the Green's function via the singular value decomposition _t,RS = ∑_j σ_j μ_j ν_j^†.
We then define each μ_j, ν_j as the available “channels" for communication and σ_j as its corresponding channel strength <cit.>. 
The compactness of  ensures that inclusion of a finite number of channels 
n_c approximates the operator to high precision (i.e., the singular values decay) <cit.>. 
This operator maps unit amplitude currents ν_j in S to fields 𝐞_R,j = σ_j μ_j in R.
Diagonalizing _t,RS^†_t,RS = U^†Σ^2 U, where Σ^2 is a diagonal matrix of squared singular values, the Shannon capacity objective of Eq. (<ref>) may then be written as <cit.> 

    log_2 (  + PNω^2Σ U  U^†Σ).

By Hadamard's inequality <cit.>, this is maximized when U  U^† is diagonal, so we bound Eq. (<ref>) by taking the degrees of freedom J_i to be the diagonal entries of U U^†.

The structural optimization problem (<ref>) then takes the simpler form

    max_ J, ϵ( r)     ∑_i=1^n_clog_2 ( 1 + γ J_i σ_i^2 ) 
    s.t.      J ≥ 0 
        1^T  J ≤ 1  
       ( ∇×∇× - ϵ( r) ω^2 )_t( r,  r',ϵ(·);ω) = ω^2 Iδ(r - r') 
       _t,RS(ω, ϵ( r)) = ∑_i^n_cσ_j μ_j ν_j^†,

where neglecting the work done by _i allows writing the problem in terms of a single free parameter γ = P/αω^2 N, and J is the vector of current amplitudes J_i. 
The utility of this formulation is now apparent: for problems where the impact of “back-action" from the photonic environment on power requirements is negligible compared to losses in the drive currents, this relaxation provides a means to investigate Shannon capacity in terms of a “scale-invariant" SNR parameter γ.

For a known structure and therefore Green's function, (<ref>) can be solved via the water-filling solution <cit.>.
The highly nonlinear dependence of the singular values σ_i on ϵ( r), however, makes this problem non-convex. Numerical solutions of (<ref>) via gradient-based topology optimization, which relaxes the discrete optimization problem over ϵ( r) to a continuous counterpart <cit.>, can lead to useful designs, but do not provide bounds on the Shannon capacity over all possible structures. To obtain a bound on the Shannon capacity that is independent of geometry, one may relax (<ref>) by ignoring the precise relation between singular values and the Green's function and instead encode structural information by imposing the same energy conservation constraints of Eq. (<ref>) on the calculation of singular values. Thus, one may write

    𝐉, σ^2max       ∑_i=1^n_clog_2 ( 1 + γ J_i σ_i^2 ) 
    s.t.        J_i ≥ 0, σ_i^2 ≥ 0,   ∀ i ≤ n_c 
       ∑_i J_i ≤ 1 
       ∑_i σ_i^2 ≤ M_𝙵
       σ_max^2 ≤ M_1

where M_𝙵, M_1 are bounds on the square of the Frobenius norm _t,RS_𝙵^2 ≡∑_i σ_i^2 and spectral norm _t,RS_2^2 ≡σ_max^2 of the channel matrix, respectively.
Each of these bounds can be calculated in a structurally independent way while maintaining crucial information on the underlying physics of Maxwell's equations (See Appendix <ref>). 
As long as constraints on σ_i^2 are bounds over all possible structures allowed in Eq. (<ref>), solutions to Eq. (<ref>) represent structure-agnostic limits on Shannon capacity between the S and R regions. 

This objective is biconvex and therefore difficult to solve <cit.>.
While numerical methods to find global optima of biconvex functions exist (see for example <cit.>),  existing implementations do not easily allow for this objective function (see for example the cGOP package). Development and application of techniques to solve for the global optimum of (<ref>) is a promising avenue of future work. 
Surprisingly, however, incorporating the Frobenius norm bound M_𝙵 but excluding the maximum singular value bound M_1 makes  (<ref>) amenable to closed-form solutions. As shown in Appendix <ref>, the optimal solution to this problem may be written as follows:  

    C(M_F) = f_n=n_opt = n_optlog_2 (1 + γ'/ n_opt^2 ) 
     n_opt = 
    r√(γ') or r√(γ'),      r√(γ')≤ n_c 
     
         n_c,      r√(γ')≥ n_c

where n_opt≥ 1, r = max_x∈ [0, 1][ x log_2 (1+1/x^2) ] ≈ 0.505, and we have introduced the total signal-to-noise ratio γ' = γ M_𝙵; note that the floor and ceiling functions are a consequence of the number of populated channels being an integer allocation.
This solution corresponds to evenly allocating capacity and power among a subset of all channels.
Notably, r√(γ')≤ 1  n_opt = 1 and r√(γ')≥ n_c  n_opt = n_c. 
In the high-noise regime (i.e., SNR dominated by increasing channel capacities), the Shannon capacity is maximized by allocating all available channel capacity and current to a single channel. 
As noise is lowered, the number of channels increases by integer amounts until the low-noise regime (i.e., SNR dominated by low noise, regardless of channel capacities) is reached, in which case every channel has capacity M_𝙵/n_c and power P/n_c.
In both cases, current allocations naturally satisfy the water-filling solution given their corresponding channel capacity allocation. 
As shown in Appendix <ref>, asymptotic solutions in the low- and high SNR regimes take the form

    C(M_𝙵)    γ' / log 2, 
     
        C(M_𝙵)     n_c log_2 (γ' / n_c^2),
 
in agreement with the analysis above. 
Intuitively, this solution is unrealistic in the low-noise regime: 
the compactness of  implies that singular values must decay, making even allocation of singular values among all channels impractical.
Regardless, this method provides a simple way to calculate limits on Shannon capacity given only a sum-rule on the Frobenius norm M_𝙵 of the channel matrix.
In the high-noise and low-noise regimes, it is straightforward to add the maximum singular value bound M_1 to this analysis; this solution is denoted C(M_𝙵, M_1).
If r√(γ')≤ 1 (the regime where n_opt = 1) a simple substitution of M_𝙵→ M_1 yields a bound. 
In general, if M_𝙵/n_opt≤ M_1, the maximum singular value constraint is not active in the solution C(M_𝙵, M_1) and so it coincides with C(M_𝙵). The amenability of this rigorous formulation to analysis via the familiar language of channel capacity and power allocation makes this conceptual ground for further studies of the impact of structuring and SNR on optimal channel design.



 §.§ Relation to Prior Work


The challenge of (<ref>) from an optimization perspective is the structural optimization over ϵ( r) which is both high-dimensional and non-convex. The main novelty of our results is the ability to account for the deliberate engineering of the propagation environment through such structural optimization, which is increasingly relevant given the rise of integrated photonics platforms. Mathematically, prior work on bounding structural optimization in photonics focused predominantly on objectives that are quadratic functions of the fields <cit.>; here we bound the Shannon capacity log objective using a relaxation with a probability distribution based interpretation that ties directly into the probabilistic proof of Shannon's noisy channel coding theorem <cit.>. 

If the structure ϵ( r) and therefore the channel matrix _t,RS is fixed, the Shannon capacity is given by a convex optimization problem that can be solved via the water-filling solution. Prior work on Shannon capacity in electromagnetism have generally made such a simplifying assumption, either taking the propagation environment to be fully vacuum <cit.> or considering a fixed antenna structure in free space <cit.>. 
Specifically, Refs. <cit.> investigated a planar antenna sheet surrounded by an idealized spherical shell receiver that captures all waves propagating through free space. 
Ref. <cit.> calculated the Shannon capacity in vacuum between concentric sphere-shell sender and receiver regions: due to the domain monotonicity of the singular values of _0 <cit.>, this bounds the capacity between smaller sender-receiver regions completely enclosed within the larger sphere-shell. 
Related works <cit.> have established sum-rule bounds on these vacuum singular values, establishing a qualitative connection between information transfer and the notion of maximizing field intensity at the receiver. 

There are also differences in the exact constraint used to restrict the input current magnitude: Ref. <cit.> implicitly used a spectral norm constraint on the input currents, which is similar to the insertion-loss dominated case of <ref>. We explicitly impose a constraint on the total power required to maintain the input currents, including not just insertion loss but also the power going into propagating fields, similar to Ref. <cit.>. Refs. <cit.> also considered further restrictions on the antenna efficiency, defined as the ratio of radiated power to total power; Ref. <cit.> included constraints on the ratio of stored vs. radiated power in the antenna as a proxy for the operating bandwidth. 
Our bound formulation (<ref>) given in the previous section is flexible enough to accommodate such additional constraints; given the conceptual richness that structural optimization adds to the problem, we have opted to only include the drive power constraint (<ref>) in this manuscript for simplicity. 



§ ILLUSTRATIVE EXAMPLES
 





In the prior section, we derived two methods to compute upper bounds on the Shannon capacity between sender and receiver regions containing any possible photonic structure in either region or their surroundings. 
As described in Section <ref>, the total drive power constraint Eq. <ref> can be considered in two regimes. 
In one regime, handled by the direct convex relaxation, the extracted power (determined by Maxwell's equations) and the insertion power (as abstracted by a resistivity associated with driving the initial currents) are comparable: 
1/2ω_t,SS_2^2 ∼α.
In the other, treated by both the direct convex relaxation and the simpler biconvex problem, the cost of driving initial currents dominates power consumption (i.e., insertion-dominated): 1/2ω_t,SS_2^2 ≪α.
In this section, we evaluate bounds in both of these contexts for several illustrative 2D configurations, showcasing the possibility of enhancing information transfer via photonic design and also revealing scaling characteristics with respect to material and design parameters. 
We compute bounds both for out-of-plane (scalar electric field) polarization (TM) and in-plane (vector electric field) polarization (TE). 

In the particular case of insertion-dominated power, we also compare the bounds with structures obtained via topology optimization, using an in-house automatic differentiable Maxwell solver written in JAX <cit.> to compute gradients of the Shannon capacity with respect to structural degrees of freedom. 
    
Direct Convex Relaxation: We consider a scenario, shown schematically in Fig. <ref>A, involving sender (S) and receiver (R) square regions of size λ×λ with centers aligned along the x-axis and separated by a surface–surface distance of 1.2λ. A rectangular mediator region M of length L and width λ (along the y dimension) is placed between S and R. Bounds are computed by numerically solving (<ref>) with energy conservation constraints imposed globally for a specified design region, i.e., _k = _D. 

To investigate the impact of structuring different regions on the Shannon capacity, Fig. <ref>B shows bounds for TM polarized sources, assuming a fixed impedance α=0.01 Zλ and SNR P/N=10^7 ϵ_0 c λ, as a function of L and for different combinations of designable S, M, and R regions. As L increases, the mediator eventually overlaps with and completely encloses the S and R regions. 
The results reveal a striking feature: in contrast to designing the sender, structuring the receiver has by far the greatest impact on the Shannon capacity. 
Intuitively, designing R allows large field concentrations _R. Such density of states enhancements are further corroborated by the dependence of the bounds on material susceptibility χ seen in Fig. <ref>C. In particular, bounds grow rapidly with χ for weak materials before saturating, and exhibit highly subdued scaling with respect to material loss (χ) consistent with resonant enhancements also seen in prior works on photonic limits <cit.>: decreasing χ by two orders of magnitude yields an increase in C of less than a factor of two. The exact nature of this scaling behavior is the subject of future work.

The apparent ineffectiveness of structuring the sender to increase C stems from the implicit freedom associated with designing the input drive, which allows for arbitrary input currents j_i within the S region, capped in magnitude only by the finite drive power constraint. In particular, the field e_R produced by a free current source j_i in a structured S region containing a polarization current j_s is in fact equivalent to one produced by a corresponding “dressed" input current j_t = j_i + j_s in vacuum. 
The imposition of a finite drive power requirement with an insertion power contribution does, however, create a distinction between these two scenarios. In the α→ 0 limit where there is no power cost to modifying the input current j_i, there is no advantage in structuring since non-zero j_s incurs further material losses compared to inserting j_t directly in vacuum. Conversely, in the insertion dominated regime α≫_t,SS_2/2ω increasing j_i has a steep cost, making it more advantageous to achieve larger field enhancements by increasing j_t via structuring. Clearly, the relative importance of insertion- versus extraction-dominated power as determined by α plays a central role in determining how C should be optimized.

Figure <ref>D shows bounds on C as a function of 1/α for both TM and TE sources, allowing for structuring only in R but fixing the distance between the S and R regions to be 1.2λ (see inset). In the extraction dominated regime (α→ 0), the Shannon capacity exhibits a logarithmic dependence on α attributable to “dark currents” which do not radiate into the far field as input signals. 
Mathematically, these dark currents lie within the nullspace of _SS, so the extracted power contribution -1/2ωj_i^†j_i to the drive power constraint (<ref>) is 0.  While dark currents do not radiate nor consume power, they nonetheless produce evanescent fields that are picked up by a receiver at a finite distance, allowing for information transfer. 
The dark current amplitudes are limited by α, so for fixed P/N their contribution to the covariance  is inversely proportional to α, resulting in a logarithmic dependence of C on 1/α via the log_2 objective.  This leads to a logarithmic divergence in Shannon capacity as α→ 0, even for finite power. As α increases, C decreases due to power restrictions on input currents, eventually reaching the insertion dominated regime (α→∞) where the extracted power is dwarfed by the insertion impedance term αj_i^2. 
Neglecting the extracted power, the resulting power constraint on the input currents ∼_i^2 ≤ P/α leads to a problem that can be bounded via both the direct convex relaxation and the biconvex relaxation of (<ref>), as illustrated below. 

Insertion-Dominated Power: 
In this section we study the insertion dominated regime (α≫1/2ω_t,SS), where one can also apply the biconvex relaxation given in Eq. (<ref>). 
As described in Section <ref>, this simplification permits a definition of channels as the singular vectors of _t,RS, with the relevant physics encoded via sum-riules on the associated singular values. This leads to a simplified, biconvex optimization problem (<ref>) with two analytical solutions: C(M_𝙵) maximizes Shannon capacity subject to a bound M_F on the maximum sum of all channel capacities (Frobenius norm), and C(M_𝙵, M_1) additionally incorporates a bound on the spectral norm M_1 or the maximum channel capacity of the Green's function. The resulting bounds depend only on the single universal total signal-to-noise ratio parameter γ' = P/αω^2 N× M_F.

We consider a simplified scenario shown schematically in Fig. <ref>, consisting of discrete S and R regions comprising 4× 4 vacuum pixel arrays separated by a designable contiguous M mediator of material χ.  Hence, there are 16 orthogonal singular vector channels available for communication. 
The figure shows bounds C(M_𝙵) and C(M_𝙵, M_1) alongside comparable inverse designs. 
Notably, the tightest bounds are found to be within or near an order of magnitude of achievable performance for all studied noise regimes, with particularly realistic performance for high SNR.
The additional tightness in C(M_𝙵, M_1) compared to C(M_𝙵) for high-noise regimes motivates extending the domain of this solution via further analysis. 
The corresponding channel occupation n_opt and power allocation P_i = 1/n_opt for C(M_𝙵) (shown above the main plot) demonstrate the nonlinear dependence of these two values on noise. 
Even allocation f_n∈[1,16] are also shown, demonstrating how the solution C(M_𝙵) is a pointwise maximum of these uniform allocations. 
As discussed earlier, higher SNR pushes allocation of J_i and channel capacity to a single channel leading to  logarithmic dependence of Shannon capacity on γ', while low SNR pushes even allocations of J_i and channel capacity across all channels, leading instead to linear scaling.  This is consistent with the results presented in Fig. <ref>, where large and small 1/α result in logarithmic and linear scaling of Shannon capacity, respectively. The biconvex formulation, however, suggests that while it may not be physically possible to optimize for any desired channel matrix, the SNR determines a target number of high-capacity channels. 
Accordingly, high SNR structures display higher symmetry than the seemingly random structure at low SNR owing to maximizing information transfer among one or many channels, respectively.

Finally, the figure also shows bounds C_convex pertaining to the convex relaxation of Eq. (<ref>) and obtained by enforcing large α and P so as to reach the insertion-dominated power regime. The latter is found to be looser than C(M_𝙵, M_1) in the low SNR regime and tighter in the high SNR regime, showcasing the versatility of the convex formulation and, despite its complexity, its general agreement with the conceptually simpler biconvex relaxation.



§ CONCLUDING REMARKS


In this article, we formalized an appropriate figure of merit for maximizing electromagnetic communication subject to arbitrary photonic structuring, and connected it to existing work restricted to vacuum propagation and/or fixed antenna geometries. 
The proposed framework rigorously connects the mathematics of electromagnetic wave propagation in photonic media to key information-theoretic quantities such as the mutual information and the Shannon capacity. 
We derived two possible relaxations of this problem (one valid in general settings and the other in regimes where insertion losses dominate) that encode wave physics from Maxwell's equations along with drive-power considerations, to arrive at structure-agnostic upper bounds on Shannon capacity. 

These bounds not only yield quantitative performance targets, but also critical scaling information with respect to device size, SNR, and material strength. In particular, our numerical evaluation of illustrative examples has highlighted the importance of optimizing the receiver, with more than an order of magnitude potential improvement of the Shannon capacity. The bounds have weak scaling with material loss, a fact with both theoretical and practical significance and likely related to similar sub-linear (logarithmic) scaling with loss seen in prior works involving distributed sources such as in near-field radiative heat transfer <cit.>. 
When radiated power costs are dominant, the Shannon capacity displays logarithmic growth as α→ 0 driven by nonradiative dark currents. 
As α is increased and insertion impedance restricts input currents, we observe capacity scaling with SNR transitioning from linear for low SNR to logarithmic for high SNR, alongside a gradual increase in the number optimal sub-channels used. Furthermore, topology optimized structures are shown to generally achieve performance within an order of magnitude of our bounds across a wide range of SNR. 

Looking ahead, we believe there is significant promise for further research into photonics design for information transfer. 
The results in this manuscript were computed in 2D and serve mainly as a proof of concept; we anticipate that our formulation may be adapted to study the information capacity of practical devices of increasing interest such as on-chip optical communication and image processing using integrated photonics. 
Doing so would also require improving computational techniques for evaluating the bounds at larger scale and in 3D, along with the development of new inverse design schemes specifically for optimizing information transfer. 
Our current results are also derived at a single frequency; further work remains to be done on generalizing the bounds for finite bandwidths, perhaps in the spirit of spectral sum rules <cit.> or delay-bandwidth products <cit.>. Another promising line of study is to incorporate macroscopic quantum electrodynamics <cit.> into the formulation for investigating information transfer in quantum communication systems.  

Funding:  we acknowledge the support by a Princeton SEAS Innovation Grant and by the Cornell Center for Materials Research (MRSEC). 
S.M. acknowledges financial support from the Canada First Research Excellence Fund via the Institut de Valorisation des Données (IVADO) collaboration.
The simulations presented in this article were performed on computational resources managed and supported by Princeton Research Computing, a consortium of groups including the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology's High Performance Computing Center and Visualization Laboratory at Princeton University.





§ PHYSICAL BOUNDS ON GREEN'S FUNCTION CHANNEL CAPACITIES
 

We will now show how the Lagrange dual relaxation as described in Refs. <cit.> can be utilized to obtain limits on the Frobenius norm and largest singular value of the total Green's function under arbitrary structuring. 
Ultimately, we seek to place sum-rule bounds on the squares of the singular values of _t,RS. 
These limits are combined with solutions to Eq. (<ref>) to find limits on the Shannon capacity. 

Unlike the prior algebraic relaxations exploiting passivity (e.g., <cit.>), these incorporate significantly richer physics through the enforcement of local energy conservation constraints and thus yield tighter limits. 
Note that the trace of _t,RS^†_t,RS can be evaluated by ∑_j v_j^†_t,RS^†_t,RSv_j in the basis of the singular vectors of _0,RS = ∑_j s_j u_j v_j^†.
Defining the vacuum “source" fields in the receiver region S_j ≡_0  v_j = s_j  u_j with amplitude s_j^2 and the polarization field in the photonic structure p_j ≡( χ( r)^-1 - _0 )^-1 S_j ≡ S_j <cit.>, we find 

    v_j^†_t,RS^†_t,RS v = v_j^†_0,RS^†_0,RSv_j + 
         p_j^†_0,RD^†_0,RDp_j 
     + 
         2v_j^†_0,RS^†_0,RDp_j

where _0,RD acts on polarization fields in the design region D (which may contain any region) and gives fields in the receiver region. 
The first term describes the square of the j-th singular value of the vacuum Green's function s_j^2, while the rest describe the contribution from the photonic structure.
In this picture, p_j,  S_j are vectors in the design region, while v_j are vectors in the source region.
Generalized energy conservation scalar identities <cit.>, can be written 
S_j^†_v  p_k -  p_j^†(χ^-† - _0,DD^†) _v  p_k = 0, ∀ j,k (see Appenidx <ref>), where _v is projection into any spatial subdomain. 
For the purposes of clarity, we take _v →, noting that constraints will be enforced in every region for the numerical calculation of bounds.
From these components, quadratically constrained quadratic programs (QCQPs) can be written to maximize different objectives related to the singular values of _t,RS, from which bounds on their optimal values can be computed with duality relaxations as described below. 

 

Frobenius norm ‖𝔾_t,RS‖_𝙵^2 = (_t,RS^†_t,RS) = ∑_i σ_i^2 ≤ M_𝙵:

    max_ p_j     ∑_j  v_j^†_t,RS^†_t,RSv_j  
    s.t.      S_j^† p_k -  p_j^†(χ^-† - _0,DD^†)  p_k = 0, 
       ∀ j,k

L_2 norm ‖𝔾_t,RS‖_2^2 = σ_max≤ M_1:

    max_v,p_j     v^†_t,RS^†_t,RSv
    s.t.     v^†_0,DD^†p - p^† (χ^-1† - _0,DD^†) p = 0,   
       v^†v = 1,  
         ∀ j ≤ k


The first problem is simply maximizing the Frobenius norm of _t,RS subject to energy conservation constraints.  In the second problem, we maximize the largest singular value by co-optimizing the source currents v and the polarization currents p and enforcing the singular vector is normalized. 
We note that in all calculations, “cross"-constraints between sources j≠ k are not enforced for computational efficiency. 
This relaxes the problem: full incorporation of these constraints is expected to tighten limits. 

To compute limits on these problems, we note that the Lagrangian ℒ of a QCQP with primal degrees of freedom ψ is given by

    (ψ, λ) =  -ψ^†𝔸(λ) ψ   + 2 ( ψ^†𝔹(λ)  S ) 
       +  S^†ℂ(λ)  S

where 𝔸, 𝔹, and ℂ represent the quadratic, linear, and constant components of the Lagrangian and contain both objective and constraint terms.  
The Lagrange dual function 𝒢 and its derivatives can be written, in terms of ψ_opt≡𝔸^-1𝔹 S and for positive-definite 𝔸,


    𝒢(λ)    =  S^†( 𝔹^†𝔸^-1𝔹 + ℂ)  S, 
    ∂𝒢∂λ   = 2 ( ψ_opt^†∂𝔹∂λ S ) - ψ_opt^†∂𝔸∂λψ_opt
         +  S^†∂ℂ∂λ S.

The positive-definiteness of 𝔸 ensures that the dual problem has feasible points with finite objective values. 
We must therefore show that for each primal problem, there exist dual feasible Lagrange multipliers λ such that 𝔸 is positive-definite. 
In the case of (1), we simply note that (_0,DD) is positive-semidefinite <cit.>. 
Therefore, for a lossy design material χ > 0, we can increase the Lagrange multiplier for the imaginary part of this constraint to ensure 𝔸 is positive-definite at some λ. All other multipliers can be initialized to zero. 
In the case of (2), the same strategy can be employed to make 𝔸 positive-definite in the p,  p sub-block. 
To ensure positive-definiteness in the v,  v sub-block, the semidefinite quadratic constraints v_i^† v_i = 1 can be employed. 
Overall, the existence of these points proves the existence of a bound on their respective primal problems. 




§ PROOF THAT THE CONVEX RELAXATION YIELDS A FINITE BOUND
 

The convex relaxation approach outlined in Section <ref> gives a bound written as the solution to maximizing an objective function of the form  log( + ^†) subject to the PSD constraint ≽0 and various linear trace constraints on . 
Upon first sight, it may be unclear whether such an optimization has a finite maximum: one can imagine that the eigenvalues of ^†≽0 may increase without bound and lead to a diverging log, and indeed this is what happens without any linear trace constraints. 
The key for a finite maximum of (<ref>) is the existence of a linear trace constraint of the form [A] ≤ P_A where A≻ 0 and P_A>0, i.e., the program

    
    max_     log( + ^†) 
    s.t.     [A] ≤ P_A 
       ≽ 0

    

has a finite maximum. 
To see this, observe that given strictly definite A, it must be the case that [A]>0 for all positive semidefinite ≠0. 
Therefore, in the vector space of Hermitian matrices, the intersection between the cone of semidefinite matrices and the halfspace given by [AJ] ≤ P_A is bounded and the log objective cannot diverge. 

In problem (<ref>), we can use a simple sum of the drive power constraint (<ref>) and resistive power conservation over the enter design domain [the imaginary part of (<ref>)] to construct such a linear trace inequality, with A given by

    [ω R _S, (1/χ^*) _D] + [(_0,SS)    (_0,SD) 
    (_0,DS)    (_0,DD)].

The first matrix is a diagonal operator with positive entries and is positive definite. The second operator is positive semidefinite since it is a sub-operator of the positive semidefinite operator

    [(_0,AA)    (_0,AA) 
    (_0,AA)    (_0,AA)]

where the domain A = S ⋃ D. To see that (<ref>) is positive semidefinite, we start from the fact that (_0,AA) is positive semidefinite due to passivity <cit.>, and therefore it has a Hermitian positive semidefinite matrix square root √((_0,AA)). It suffices now to note that 

    [x
    y]^†[(_0,AA)    (_0,AA) 
    (_0,AA)    (_0,AA)] [x
    y]^†
       = √((_0,AA)) (x + y)_2^2 ≥ 0.


It is also worth pointing out that if the constraint matrix  in (<ref>) is indefinite, then the optimum diverges. To see this, consider any vector h outside the kernel of , i.e., h≠0. If h^†h < 0, then we are done: = βhh^† will lead to a divergence of the objective as β→∞ while [] = h^†h→ -∞ < P_A. If h^†h≥ 0, simply pick a vector v such that v^†v < 0; this is always possible given indefinite . Then by constructing = β_1(hh^† + β_2 vv^†) and choosing β_2>0 such that h^†h + β_2 v^†v < 0 the same argument as before works: as β_1 →∞ the objective diverges thanks to the hh^† contribution while ≻ 0 remains within the constraint halfspace. 



§ CONVEX RELAXATION SOLUTION VIA WATERFILLING DUALITY
 

This appendix describes the details of solving the direct convex relaxation (<ref>) using an approach we dub the “waterfilling dual”, which is a generalization of the solution technique used in <cit.>. Specifically, we seek the solution of a convex optimization problem over the n × n PSD matrix :

    
    C = max_     f() ≡log( + ^†) 
    s.t.     [] ≤ P 
       [_l ] = 0    l ∈ [1⋯ m] 
       ≽0

    

where ≻ 0 and we assume that the solution  is non-zero, i.e., the _l and PSD constraints do not restrict =0 and the PD  constraint is needed for C to be finite. While (<ref>) is not explicitly of this form, Appendix <ref> shows how one can straightforwardly transform (<ref>) into this form by taking linear combinations of constraints to obtain a PD constraint. The waterfilling dual approach to (<ref>) starts from the fact that a simplified version of (<ref>) without any [_l ]=0 constraints admits an analytical solution, i.e., waterfilling. Unfortunately, waterfilling does not directly generalize to multiple linear constraints, so we define the waterfilling dual function

    
    𝒟(b) = max_     log( + ^†) 
    s.t.     [( + ∑_l=1^m b_l _l) ] ≤ P, 
       ≽0.

    

The trace constraint in (<ref>) is a linear superposition of the trace constraint in (<ref>), so 𝒟(b) ≥ C for any combination of multipliers b; furthermore, the value of 𝒟(b) can be evaluated using waterfilling. 𝒟 can be understood as a partial dual function of (<ref>) where we have not dualized the ≽0 constraint; 𝒟(b) ≥ C is thus a statement of weak duality. We now show explicitly that, as expected for a convex problem, strong duality holds: the solution to the waterfilling dual problem

    min_ b  𝒟( b)

is exactly equal to C. To see this, consider the objective gradient at the primal optimal point f(). First order local optimality conditions imply that

    f() = (∑_l=1^m b̅_l _l) + a  + ∑ s_k S_k

where a, s_k ≥ 0, and S_k are outward normal matrices of the PSD cone 𝒮^n_+ locally at : such a local description of 𝒮^n_+ is possible due to the fact that any non-zero matrix on the boundary of 𝒮^n is an interior point of a face of 𝒮^n_+ <cit.>. Since by assumption the primal problem has a non-zero solution, the PD constraint must be active and therefore a > 0. Now, consider the dual value 𝒟(b̅/a) as the optimum of (<ref>): it is clear that  is a feasible point of (<ref>), and furthermore by construction f() satisfies the first order local optimality condition

    f() = a ( + ∑_l=1^m b̅_l/a_l ) + ∑ s_k S_k.

Conclude that 𝒟(b̅/a) = f() = C and therefore strong duality holds. Note also that for 𝒟( b) to remain finite, + ∑ b_l _l must be positive semidefinite per Appendix <ref>, so we can always do waterfilling to evaluate finite 𝒟 values. Finally, given the analytic form of the waterfilling solution of (<ref>) one can also derive dual gradients ∂𝒟 / ∂ b using matrix perturbation theory, and thus the dual problem (<ref>) can be solved using gradient-based optimization. 



§ GENERAL SOLUTIONS TO SHANNON CAPACITY RELAXATION
 
In this section, we provide a general solution of (<ref>).
The problem can be written

    x,  ymax         f( x,  y) = ∑_i=1^nlog_2 ( 1 + γ x_i y_i ) 
    s.t.        1^T  x =  1^T  y = 1 
        x,  y ≥ 0

where the c and prime in n_c, γ have been dropped, respectively. Let x^*, y^* be optimal solutions to this problem, which exist because the function is continuous over the compact set defined by the constraints. 
First, we consider asymptotic solutions in γ.

 

γ≪ 1:
∑_i log_2 (1+γ x_i y_i) ≈γ/log 2 ∑_i x_i y_i. Given the constraints 1^T x = 1 and 1^T y = 1 (easy to see the optimum must have equality), and that log(1+x) ≤ x, we get the bound

    C ≤γ / log 2

using the method of Lagrange multipliers, first on the maximization with respect to x_i and then y_i. This is a bound for all γ but is most accurate when γ is small. 

 

γ≫ 1:
In this case, we find that 
∑_i log_2 (1+γ x_i y_i) ≈∑_i log_2 x_i + log_2 y_i + log_2 γ.
This maximization problem can also be analytically solved given the same constraints by the method of Lagrange multipliers to find

    C →  n_c log_2 1n_c + n_c log_2 1n_c + n_c log_2γ = n_c log_2 γn_c^2.

This corresponds to equal allocation 1/n_c and 1/n_c for all channels. 
This value approaches a bound on C as γ→∞. 
Now, we show how this problem can be solved exactly, analytically.
In the following analysis, we take log_2 →log for mathematical convenience, understanding that solutions are scaled by a constant factor. 


    Let S_x^* = { i ∈ [n]: x_i^* > 0} and S_y^* = {i∈ [n] : y_i^* > 0} be the support sets of x^* and y^* respectively. Then S_x^* = S_y^*≡ S.
    


    Pick two indices i≠ j. If x_i = 0, y_i = c, then taking y_i → 0 and y_j → y_j + c strictly increases the objective. 
    Intuition: allocating mass to some y_i if x_i=0 contributes 0 to the objective, whereas allocating that same mass to another y_j contributes positively. Thus the support sets of x^* and y^* are equal.  


 
    For any pair of indices i, j ∈ S,  x_i^* y_i^* x_j^* y_j^* = 1/γ^2 (Condition 1) or x_i^* = x_j^* (Condition 2).


Fix y^* and consider the following optimization problem 

    max_x       f^γ,y^*( x) = ∑_i =1^n log( 1 + γ x_i y_i^* ) 
    s.t.        1^T  x = 1, 
        x ≥ 0.

We know that x^* is a KKT point in the above optimization problem.
The KKT conditions are, for some λ > 0, γ_1, …ν_n ≥ 0,

    λ 1 - ∑_i=1^n ν_i  1_x_i^* = 0   = ∇ f_γ,y^*(x^*),

where 1_x_i^* = 0 is an indicator function for x_i^* = 0.
Now x_i^* = 0 if and only if y_i^* = 0 by Lemma <ref>. 
If y_i^* = 0 then it is easy to see that ∇ f_i^γ,y^*(x^*) = 0, and therefore ν_i = λ. Otherwise, 1_x_i^* = 0 is 0 everywhere. 

 

Therefore, for each pair of i,j∈ S, 

    y_i^* γ1+γ x^*_i y_i^* = γ y_j^* 1+γ x^*_j y_j^*.

Solving for x^*_j - x^*_i, we obtain

    x^*_j - x^*_i = (y_j^* - y_i^*)γ y_j^* y_i^*.

Making the same argument fixing x^*, we find 

    y_j^* - y_i^*    = (x_j^* - x_i^*)γ x_j^* x_i^*

Combining these, we find that 

    x_j^* - x_i^* = 1/γ^2x_i^* y_i^* x_j^* y_j^* (x_j^* - x_i^*),

from which it follows that either x_i^* y_i^* x_j^* y_j^* = 1/γ^2 or x_i^* = x_j^*.



    For every i, x_i^* = y_i^*.




Enumerate the unique nonzero values of x_i^* for i=1… n by X_1, … X_K. Suppose K ≥ 2 (if K=1 the result is clear). For i=1 … k let S_k = {i ∈{1 … n} such that  x_i=X_k} for k ∈{1 … K}.
By (<ref>) if x_i^* = x_j^* then also y_i^* = y_j^*. Denote the unique value of y_i^* in the set S_k by Y_k.

 

Now choose i ∈ S_k  for some k ≥ 2 and j ∈ S_1. From (<ref>) combined with the fact that x_i^* y_i^* x_j^* y_j^* = 1/γ^2, we observe

    γ y_i^*    = γ y_j^* ( 1 + 1/(γ x_j^* y_j^*)/1 + γ x_j^* y_j^*)

and symmetrically in x,

    γ x_i^*    = γ x_j^* ( 1 + 1/(γ x_j^* y_j^*)/1 + γ x_j^* y_j^*).

Now divide the above two equations to see that for every k,

    X_k/Y_k   = X_1/Y_1.

Since 1=∑_i=1^n x_i^* = ∑_i=1^n y_i^*, this implies that x_i^* = y_i^* for each i such that x_i^* > 0.

Now we know that for each i≠ j, X_i^2 X_j^2 = 1/γ^2. This is impossible if K ≥ 3 (xy=c,yz=c,xz=c  y=z) so we know K=1 or K=2. 
First, we consider K=1.
 
    Let g(x) = x log (1 + 1/x^2). There is a unique local maximizer of g over the real numbers, which is also the global maximizer, and is between 0 and 1. Let this real number be r.



    g'(k→ 0) →∞ and g'(1) = log(2) - 2/3 ≤ 0. For x ≥ 1, g'(x) = log(1 + 1/x^2) - 2/(x^2 + 1) ≤ 1/x^2 - 2/(x^2 + 1) = (1 - x^2) / (x^2(x^2 + 1)) ≤ 0. Also,
    
    g”(x)    = 2x/x^2 + 1 - 2x/x^2 + 4x/(x^2 + 1)^2 = 2x(x^2 - 1)/x^2(x^2 + 1)^2

    so g”(x) ≤ 0 for x ∈ [0,1] and therefore g is concave. So there can be only one local maximizer in this interval.


 
    Let h(k) = k log (1 + γ/(k^2)). The unique maximizer of h over the integers is in the interval [⌊ r √(γ)⌋, ⌈ r √(γ)⌉ ]. 
    



    h(k/√(N)) = √(γ) g(k) which is maximized over ℝ at r√(γ). Since the function is concave on the relevant interval, the solution is one of ⌊ r √(γ)⌋, ⌈ r √(γ)⌉.


Since h(k) is concave on [0,n], if r√(γ)≥ n then the optimal k is n. Therefore we show that the K=1 solution to (<ref>) is

    n_optlog(1 + γn_opt^2) 
       n_opt = 
    r√(γ) or r√(γ)     r√(γ)≤ n 
     
         n      r√(γ)≥ n

where n_opt≥ 1, r = max_x∈ [0, 1][ x log(1+1x^2) ] ≈ 0.505 and the floor and ceiling functions are a consequence of the number of populated channels being an integer allocation.

 

Now we consider K=2 solutions.
Suppose that |S_1| = n_1 and |S_2| = n_2 such that n_1 + n_2 ≤ n. 
Observe that X_1 = λ/n_1 and X_2 = (1-λ)/n_2 for some λ∈ [0,1]. 
Consider the function,

    f(λ) = n_1 log( 1 + λ^2 γ/n_1^2) + n_2 log( 1 + γ(1 - λ)^2/n_2^2).

Any maximizer in the K=2 case must be of this form. 
We can explicitly find the maximizing λ for this function in terms of n_1 and n_2.


    The values of λ which can be a maximum of f in [0,1] are λ= 0, λ = 1,
    
    λ_0 = n_1/n_2 + n_1      or   λ_± = 1 ±√(1 - 4 n_1 n_2 / γ)/2.




    These are the values of λ which have zero derivative, and those at the extremes of λ∈ [0, 1].


These are the possible candidates for being a maximum of our function when K=2. The first and second cases correspond to K=1. In the third case, X_1 = 1/(n_1 + n_2) = X_2, so in fact K = 1. 
In the fourth case, we can assume without loss of generality that n_1 ≥ n_2, meaning λ_+ is the optimal K=2 point (as it must be between local minima of λ_0 and λ_-).
In practice, we numerically compare the K=1 with all K=2 solutions to get the global optima of this problem. 
In theory, we believe it is possible to prove that for all K=2 solutions, there exists a K=1 solution with a higher objective value. 

We now prove this for most values of N. 
Towards a contradiction, we assume that λ_+ ≡λ, n_1, n_2 in (<ref>) represent a global maximum of the function (<ref>).
We will now show that λ_+ representing a global maximum leads to a contradiction. That is, we can either change λ, n_1, or n_2 to a different value to get a better objective value.

 

Case 1: 1 ≥max{ (1/√(γ)) n_1 /λ, (1/√(γ)) n_2 / (1 - λ) } 

By the concavity of h (on the interval [0, 1]), we find that for λ such that 1 ≥n_1 / √(γ)λ and 1 ≥n_2 / √(γ)1 - λ,


    f(λ)    = λ h( n_1/λ) + (1 - λ) h( n_2/1 - λ)
       ≤ h(n_1 + n_2)

 meaning we can get a better optimal value by summing n_1 and n_2 and considering a point in the K=1 case, giving a contradiction.

 

Case 2: 1 ≥ 1/√(γ) n_1 /λ, but 1 ≤ (1/√(γ)) n_2 / (1 - λ), and n_2 ≥ 2. First, we prove the following lemma:
 
    For x ≥ 1, g(x/2) ≥ g(3x / 4) ≥ g(x) and g is strictly decreasing on [3x / 4, x].



    We already showed that for x ≥ r g is decreasing, and r ≤ 3/4 so the second inequality is trivial. 

    For the first inequality, for x ≥ 2 r again it is clearly trivial. We just need to show that it is true for x ∈ [1, 2 r].
    Observe that g(1/2) = 1/2 log(5) ≥ 0.7 and for x ≤ r, g is decreasing. Therefore for any x ∈ [1/2, r], g(x) ≥ 0.7. Also for any x ≥ 3/4, again using the decreasing property of g, g(x) ≤ g(3/4) ≤ 0.6.

Now, if n_2/√(γ)1 - λ≥ 1 then (n_2 - 1)/√(γ)1 - λ≥ 1/2 so, using Lemma <ref>, λ h(n_1 / λ) + (1 - λ) h (n_2 / (1 - λ)) ≤λ h(n_1 / λ) + (1 - λ) h ((n_2 - 1) / (1 - λ)).
Therefore, n_2 → n_2 - 1 increases the objective value, violating the assumption of optimality. 

 

Case 3: 1 ≥  n_1 / (√(γ)λ), but 1 ≤ (n_2 / √(γ)) / (1 - λ), and n_2 = 1. 

 

Note that in this case we must have n_1 = ⌊ r √(γ)⌋ / √(γ) because otherwise n_1  n_1 + 1 would yield an objective improvement. 

 

First we must observe that n_2 / (√(γ)λ) ≥ 2. This can be done by numerically solving the following bounded three-dimensional program and observing that the optimal objective value is less than 0,

    x ≤ 1, y ∈ [1, 2], λ≥ 0.5maximizeλ g(x) + (1 - λ) g(y) - g(λ x + (1 - λ) y ).

Now using this fact, plus the fact that g is concave on [0, 1],

    λ g ( n_1/√(γ)/λ) + (1 - λ) g ( 1/√(γ)/1 - λ)  
       ≤λ g ( n_1/λ√(γ)) + (1 - λ) g (2) 
       = λ[ g(n_1 /√(γ)) + g ( n_1/λ√(γ)) - g(n_1√(γ)) ] + (1 - λ) g(2) 
       ≤λ[ g(n_1 /√(γ)) + g'(n_1 /√(γ)) ( n_1/λ√(γ) - n_1√(γ)) ] + (1 - λ) g(2) 
        = g(n_1√(γ)) λ + (1 - λ) n_1√(γ) g'(n_1 /√(γ)) + (1 - λ) g(2).

So an optimal point of the K=2 type is only possible if

    mmax g(m/√(γ)) 
       ≤ g(n_1 /√(γ)) γ + (1 - γ) (n_1 /√(γ)) g'(n_1 /√(γ)) + (1 - γ) g(2)

or,

    g'(n_1 /√(γ))    ≥(mmaxg(m /√(γ)) - λ g(n_1 /√(γ))/1 - λ- g(2)) √(γ)/n_1
       ≥ (g(n_1 /√(γ))- g(2)) √(γ)/n_1.

We can now look at x for which,

    g'(x) ≥ (g(x) - g(2)) 1/x.

These are x ≤ 0.235 and x ≥ 4.245 (solved via a numerical solver). This means that to get an optimal point of this form we need ⌊ r / √(N)⌋√(N) = n_1 √(N)≤ 0.235 which, since r > 0.5, is impossible.
 


In conclusion, if both n_1√(γ)λ and n_2 / √(γ)1-λ≤ 1, by case 1 there exists a value with K=1 with a better objective value. If either n_1λ√(γ) or n_2 / √(γ)1 - λ are greater than 1 (WLOG the second case), with n_2 ≥ 2, by case 2 this allocation is not optimal. 
Finally, if both n_1/√(γ)λ and n_2 / √(γ)/1 - λ are greater than 1, then numerically optimal K=2 points are not possible. 
K=2 points are nevertheless numerically checked to not be optimal in the examples presented in the main text. 



§ ENERGY-CONSERVING QCQP CONSTRAINTS
 

It is possible to derive the energy conservation constraints presented in the main text via operator constraints (which, like Maxwell's equations are dependent on the spatial profile of permittivity) relaxed to structure-agnostic scalar constraints <cit.>. 
Although these constraints encode less information than Maxwell's equations in full generality, they allow the relaxed problem to be formulated as a QCQP from which techniques in the main text (i.e., Lagrange duality) can be employed to calculate structure-agnostic bounds on photonic performance. 
In this section, we show how these constraints can be directly derived from energy conservation as expressed by Poynting's theorem for time-harmonic complex fields <cit.>:

    ∫_∂ Vσ· (EH^*) = iω∫_V (H^*·μ·H -  E·ϵ^*·E^*)   V 
    
        - ∫_V E·J^*   V.
 
For simplicity, we will assume a non-magnetic material μ=μ_0=1 and scalar isotropic permittivity and susceptibility ϵ = 1 + χ, though the derivation is valid for anisotropic ϵ as well <cit.>. 

Consider a scattering theory picture where a free current source J_v generates the fields E_v, H_v in vacuum and E_t, H_t in the presence of a structure with material distribution ϵ(r) = 1 + χ_s(r) where _s(r) is an indicator function. There is an induced polarization current J_s in the material which produces scattered fields E_s and H_s that combine with the vacuum fields to give the total field: E_t = E_v+E_s, H_t = H_v + H_s. The complex Poynting theorem thus (<ref>) applies to three sets of currents, fields, and environments: (J_v, E_v, H_v) in vacuum, (J_s, E_s, H_s) in vacuum, and (J_v, E_t, H_t) over the structure, giving

    ∫_∂ V_k σ·(E_vH^*_v) = iω∫_V_kH_v^* ·H_v   V 
    - iω∫_V_kE_v ·E_v^* - ∫_V_kE_v ·J_v^*   V .



    ∫_∂ V_k σ·(E_sH^*_s) = iω∫_V_kH_s^* ·H_s   V 
    - iω∫_V_kE_s ·E_s^*   V - ∫_V_kE_s ·J_s^*   V .



    ∫_∂ V_k σ·(E_t H^*_t) 
        = iω∫_V_kH_t^* ·H_t   V 
    - iω∫_V_k (1 + χ^*_s) E_t ·E_t^*   V - ∫_V_kE_t ·J_v^*   V .

Subtracting (<ref>) and (<ref>) from (<ref>) gives

    { i ω∫_V_kH_v^* ·H_s   V - ∫_∂ V_k σ·(E_s H_v^*) 
     - iω∫_V_kE_s ·E_v^*   V  }
    
        +{i ω∫_V_kH_s^* ·H_v   V - ∫_∂ V_k σ·(E_v H_s^*) 
     - iω∫_V_kE_v ·E_s^*   V  }
    
        = ∫_V_kE_s ·J_v^*   V - ∫_V_kE_s ·J_s^*   V + iω∫_V χ^* E_t ·E_t^*   V .

Now, using vector calculus identities along with the Maxwell wave equations E_v - ω^2 E_v = iωJ_v and E_s - ω^2 E_s = iωJ_s, the two curly brackets in (<ref>) can be shown to be equal to ∫_V_kE_s ·J_v^*   V and ∫_V_kE_v ·J_s^*   V respectively. Finally, the induced current J_s can be swapped out by the polarization p via J_s = -iωp, and the scattered field E_s = Gp, to give

    ∫_V_kE_v^* ·p  V = ∫_V χ^-1*p^*·p  V - ∫_V_kp^* · ( G^†p)   V.

This can be written in a compact operator notation

    E_v^†_V_kp = p^† (χ^-† - G^†) _v_kp,

giving a form of the energy conservation constraints in the main text, where E_v → S. 
From this derivation it is clear that the constraint (<ref>) encodes conservation of power during the electromagnetic scattering process for every region V_k. 
In the case of many sources (as in this paper), each source defines a different scattering problem, with individual source-polarization pairs (S_j, p_j) satisfying constraints of the form (<ref>). There are also additional “cross-constraints” that capture the fact that the same structured media generates the p_j induced in each case:

    S_j^† p_k -  p_j^†(χ^-† - G^†)  p_k = 0,   ∀ j,k.

For a more detailed discussion of these constraints, we refer the reader to <cit.>.
For computational simplicity, only j=k constraints are enforced in the bounds calculated in this paper, although constraints are enforced for all computational voxels V_k; future incorporation of these constraints is expected to tighten bounds. 







