Inverse design of plasma metamaterial devices for optical computing
Jesse A. Rodriguez, Ahmed I. Abdalla, Benjamin Wang, Beicheng Lou, Shanhui Fan, Mark A. Cappelli
IInverse design of plasma metamaterial devices for optical computing
Jesse A. Rodr´ıguez, ∗ Ahmed I. Abdalla, Benjamin Wang, Beicheng Lou, Shanhui Fan, and Mark A. Cappelli Department of Mechanical Engineering, Stanford University, Stanford, CA 94305 Department of Physics, Stanford University, Stanford, CA 94305 Department of Electrical Engineering, Stanford University, Stanford, CA 94305 (Dated: February 11, 2021)We apply inverse design methods to produce two-dimensional plasma metamaterial (PMM) de-vices. Backpropagated finite difference frequency domain (FDFD) simulations are used to designwaveguides and demultiplexers operating under both transverse electric (TE) and transverse mag-netic (TM) modes. Multiplexing and waveguiding are demonstrated for devices composed of plasmaelements with reasonable plasma densities ∼ in-situ training and experi-mental realization of these designs. We also explore the possible applicability of PMMs to nonlinearboolean operations for use in optical computing. Functionally complete logical connectives (NORand NAND) are achieved in the TM mode. I. INTRODUCTION
Inverse design methodology is an algorithmic techniqueby which optimal solutions are iteratively approximatedfrom an a priori set of performance metrics. When ap-plied to electromagnetically-active systems [1–18], de-vices are built by designating a design region whereinvarious numerical methods are used to solve Maxwell’sequations for a given distribution of sources and materialdomain. The constitutive properties (e.g., relative per-mittivity or permeability) of the design region are thenmapped to trainable parameters which are modified tomold the propagation of source fields through the des-ignated region in a manner that fulfills the user’s per-formance criteria. Typically, the performance criteria isencoded as a objective function that needs to maximizedor minimized. In many problems the objective functionsinvolve an inner product between the realized and the de-sired electromagnetic fields in regions of interest. GivenMaxwell’s equations as a constraint, automatic differenti-ation allows the calculation of exact numerical gradientsrelating the chosen objective function to the parametersthat encode the permittivity structure of the training re-gion. The algorithm then iteratively adjusts these pa-rameters to maximize the overlap between the realizedand desired modes of propagation [1–3]. The result isa novel device geometry specialized to the user’s perfor-mance criteria.Unfortunately, the set of physically realizable devicesfor inverse design methods is a small subset of the con-figuration space, leading to intensive investigation of de-vices which are amenable to these constraints such asphotonic crystal-style devices [1, 4–7] and metasurfaces[8–11]. In general these limitations require that certainrestrictions be applied during the training process. Forexample, utilizing a continuous range of permittivitiesthroughout the simulation domain greatly improves theability of the algorithm to extremize the objective func- ∗ [email protected] tion, but continuity is not easily realizable with conven-tional manufacturing techniques and materials. To re-solve this issue, a map from the parameterized designregion to the set of actually manufacturable device con-figurations is used. These maps encode restrictions onthe range and spatial distribution of the domain permit-tivity. A common technique in the literature is the useof nonlinear projection to design less accurate, but moreeasily manufacturable, binarized photonic devices [12].Such a device can be printed or created via lithographyfor uses in one-shot computing, but would be otherwiselimited in general computing applications.Rather than restricting the range of permittivities, wefocus instead on limiting the spatial configuration by pa-rameterizing the design region as a plasma metamate-rial (PMM). Our PMM is a periodic array of cylindri-cal gaseous plasma elements (rods) which can be experi-mentally realized via the use of discharge lamps or laser-generated plasmas. PMMs are ideal candidates for in-versely designed devices since plasma density, and hencepermittivity of each element within the PMM, can bedynamically and precisely tuned through a continuousrange of values ( ε ∈ ( −∞ ,
1) for a collisonless plasma).This yields an infinite configuration space for trainingpurposes [19] and allows the possibility of a single devicewhich serves multiple functions. Past work with inversedesign electromagnetic devices achieved reconfigurabilityvia methods such as refractive index changes [8] in ma-terials like liquid crystals [10]. In our case, the relativepermittivities of elements within the PMM device can bedynamically tuned by varying the gas discharge current.In addition, when the PMM is composed of elements thatare small compared to the operating wavelength, we canaccess localized surface plasmon (LSP) modes along theboundary between a positive permittivity backgroundand a negative permittivity plasma element. Such sur-face modes can yield more complex and efficient powertransfer [20, 21], allowing for a particularly high degree ofreconfigurability. Prior work by authors Wang and Cap-pelli with plasma photonic crystals (PPCs) has alreadyhighlighted the richness of this geometry in waveguide[22] and bandgap [20, 23, 24] devices. The rich and var- a r X i v : . [ phy s i c s . a pp - ph ] F e b ied physics of PMMs is well matched to inverse designmethods which promise to allow more holistic and effi-cient exploration of the configuration space.The design processes for binarized and PMM geome-tries are contrasted in Fig. 1. Binarized geometries (Fig.1a) are challenging to optimize because of the requisitediscontinuities. The jumps in permittivity from one ma-terial to another are difficult to handle in gradient-basedmethods that require the permittivity to vary continu-ously. Binarization can instead be approximated by anynumber of techniques, including gaussian blurring andhighly nonlinear projections. In these configurations, thetrainable parameters ρ (of dimension n x ∗ n y ) maintain aone-to-one correspondence with the individual pixels ofthe training region. In contrast, the number of trainingparameters in our PMM (Fig. 1b) is determined by thedimensions of the plasma rod array, where there is a one-to-one correspondence between each plasma element andthe n x ∗ n y elements of ρ (Fig. 1c). In the end, bothparameterization schemes result in a relative permittiv-ity matrix ε r (of size N x × N y ) containing permittivityvalues at each pixel in the simulation domain. FIG. 1. Flow chart describing the algorithmic design of ei-ther (a) binarized material with penalized material use or (b)array of plasma columns either in the TM mode ( E x ) or thethe TE mode ( E z ). Though the parameterization stages aredifferent, both configurations can be given identical sourcesand objective criteria. J represents the modal source for theFDFD simulation. Examples of the training matrix ρ , per-mittivity matrix (cid:15) r and simulated E fields for a simple PMMdevice are provided in (c). II. METHODS
Fig. 2 provides a diagram of our PMM device in adirectional waveguide configuration. Subsequent configu-rations exhibit the exact same array of tunable plasma el-ements but with different entrance/exit waveguides. Thedevice is simply a 10 ×
10 array of tunable plasma rodssuspended in air and spaced according to the limits ofexisting experimental facilites [22, 23]. Modal sourcesare introduced at the input waveguide(s) and allowed toscatter through the training region before again collect-ing at the desired output waveguide(s). In this work,the fields are propagated through the domain via finitedifference frequency domain (FDFD) simulations com-puted with Ceviche, an autograd compliant electromag-netic simulation tool [2]. Note that the use of FDFDimplies that all computed devices represent the steadystate solution achieved after some characteristic time.The simulation domain was discretized using a resolutionof 50 pixels per lattice constant a and in each case pre-sented N x = N y = 1000 pixels. The only exceptions arethe logic gate devices which were optimized at a resolu-tion of 30 pixels /a due to the high cost of the simulation.In those cases, the final device configurations were testedat the higher resolution of 50 pixels /a and yielded resultsthat were qualitatively identical. A perfectly-matchedboundary layer (PML) 2 a in width was applied alongthe domain boundaries. The polarization of the inputsource has a strong effect in devices of this nature, ei-ther E z ( E out of the page), which we call the TM case,or E x ( H out of the page), which we call the TE case.Previous work with PMMs indicates that both the TMand TE responses are tunable, with the latter benefit-ting from the presence of LSP modes [20, 25], while theformer makes more direct use of dispersive and refrac-tive effects [21–23]. The simulated E x and E z fields aremasked to compute the field intensity and mode over-lap integrals along planes of interest within the problemgeometry. These integrals are used to calculate the ob-jective L ( ρ ). Ceviche then computes numerical gradientsof the objective with respect to the training parametersvia forward-mode differentiation [2] and uses the Adamoptimization algorithm [26] to iteratively adjust ρ andthereby maximize L . Optimization was conducted withlearning rates ranging from 0 . − . β = 0 .
9, and β = 0 . ρ L ( E ) given ∇ × µ ∇ × E − ω ε ( ω, ρ ) E = − iωJ. where L is the objective composed of a set of mode over-lap integrals ( L inner product of simulated field withdesired propagation mode), ρ is an n x ∗ n y -dimensionalvector that contains the permittivity values of each ofthe PMM elements, µ is the vacuum permeability, E isthe electric field, J is a current density used to define a r a PML BC
Source Probes ra = 0.433 Metal Waveguide ε rod ε bg = 1 ˜ ω = 1 FIG. 2. Schematic diagram for the PMM waveguide deviceconsisting of a 10 ×
10 array of gaseous plasma rods suspendedin air with metal waveguides functioning as the entrance andexit(s). The rod radius is chosen to be commensurate withcurrent experimental facilities. The ’probe’ slices representwhere the mode overlap and intensity integrals that constitutethe objective function are calculated. ε bg is the backgroundpermittivity. ˜ ω is the nondimensionalized operating frequency ω/ (2 πc/a ). fundamental modal source at the input waveguide, ω isthe field frequency, and ε ( ω, ρ ) is the spatially-dependentpermittivity that is encoded by ρ . In practice, this per-mittivity distribution is achieved by varying the plasmadensity through control of the discharge current in eachof the PMM elements according to the Drude model,which is where the dependence on ω arises. The non-dimensionalized plasma frequencies in ρ are mapped toelement permittivities via the Drude model with no col-lisionality; ε = 1 − ω p ω where ω p = n e e ε m e is the plasma frequency squared, n e isthe electron density, e is the electron charge, m e is theelectron mass, and ε o is the free-space permittivity. III. RESULTS AND DISCUSSIONA. Directional waveguide
Electromagnetic waveguides are fairly common in in-verse design devices either as the primary function oras important building blocks to more complicated de-vices [5, 18]. Presented first in Fig. 3 are field simu-lations and permittivities of straight waveguide devicesdesigned for either the TM or TE polarization. Operat-ing at ω = 1 × (2 πc/a ), or a non-dimensionalized sourcefrequency ˜ ω = ω/ (2 πc/a )=1, these represent the sim-plest type of optimization problem for this configuration. In both polarization cases, the domain was initializedwith r/a = 0 . ε bg = 1, and ε rod = 0 .
75 for every plasma element in the array. Theplasma frequency of each of the rods was unconstrainedduring training ( ω p ∈ [0 , ∞ ]). The rod permittivitieswere initialized at 0.75 to allow the source to reach theexit waveguides in the first training epoch which consis-tently led to higher field intensities at the desired exitafter training. In this and all other simulations, the per-mittivity of the input and output waveguide walls was setto ε = − L inner product of the simulated field withan m = 1 propagation mode at the location of the probein the desired exit waveguide, L wvg = (cid:90) E · E ∗ m =1 dl desired exit . No penalty was needed to achieve a decent standard ofobjective functionality. We can see from Fig. 3 that boththe TE and TM modes yield good performance, with aslightly better result in the TM mode straight waveguide.The optimized distribution of plasma densities (hence rel-ative permittivities) in the plasma rods for the TM case isnot surprising, representing a rectangular bridge of lowerrefractive index than its surroundings. However, the TEcase settled on a distribution that was unintuitive. TE ⃗ E w TM ⃗ B w (a)(b) (a.i) (a.ii)(b.ii)(b.i) FIG. 3. (a) TE FDFD simulation showing (a.i) the field mag-nitude | H z | (ˆ z out of the page) for an optimized straightwaveguide and (a.ii) the relative permittivity domain thatgives rise to this behavior, where the maximum plasma fre-quency among the plasma elements is ˜ ω p = 1 .
17. (b) FDFDsimulation for the TM case showing (b.i) | E z | and (b.ii) therelative permittivity with a maximum plasma frequency of˜ ω p = 1 . Next, we defined an objective intended to simulatea waveguide with a 90-degree bend. The domain wasinitialized almost identically to the straight waveguide,the only difference being that the objective overlap in-tegral was computed at the other exit waveguide. Nopenalty was imposed at the incorrect exit waveguide asit was deemed unnecessary for satisfactory performance.The results for this case are found in Fig. 4. Again wesee that, while both cases result in good device perfor-mance, the TM case again slightly outperforms the TEcase which had some field leakage into to the undesiredoutput waveguide. In both modes, the distributions set-tled on by the training process were wholly unintuitive. TE ⃗ E w TM ⃗ B w (a)(b) (a.i) (a.ii)(b.ii)(b.i) FIG. 4. (a) TE FDFD simulation showing (a.i) the field mag-nitude | H z | (ˆ z out of the page) for an optimized bent waveg-uide and (a.ii) the relative permittivity domain that givesrise to this behavior, where the maximum plasma frequencyamong the plasma elements is ˜ ω p = 1 .
27. (b) FDFD simula-tion for the TM case showing (b.i) | E z | and (b.ii) the relativepermittivity with a maximum plasma frequency of ˜ ω p = 1 . The devices presented in Figs. 3 and 4 appear to makeexpert use of the underlying electromagnetics. Inversedesign allows a nuanced search of the configurationspace, yielding more effective designs that preservethe input mode. Again, we emphasize that the finaltraining results are not intuitive configurations. Ahuman-reasoned design for this device may more closelyresemble that of Wang et al. in ref. [22], where theelements along the desired path of propagation aremade equivalent to the background permittivity (i.e.turned off) and the remainder of the rods are activatedat the lowest permittivity (highest plasma frequency)possible to approximate a metallic waveguide. Instead,the algorithm settled on a configuration where a highcontrast in permittivity of neighboring plasma columnsis only present near the edges of the device to preventleaking. The remaining elements are given a seeminglyrandom distribution of plasma frequencies, keeping inline with typical results of inverse design schemes wherestrange structures often arise [12, 15]. Fig. 5 graphsthe evolution of the objective while training for thesewaveguide devices. The curves suggest that a localmaximum was achieved in each optimization scheme.The spurious changes in the TE objectives are likelydue to the excitation of LSPs. Once the permittivity ofthe plasma elements is pushed to negative values by theparameter evolution, small changes in plasma frequency can excite localized surface plasmon resonances thatstrongly affect the overall wave propagation.
FIG. 5. Evolution of the objective functions for (a) the TEwaveguides in the (a.i) straight waveguide case and (a.ii) thebent waveguide case. The same is presented for (b) the TMwaveguides. The y-values of these plots are arbitrary as theyare normalized by the initial (epoch 0) value.
Overall, the results depicted in Figs. 3 and 4 are en-couraging. Unlike binarized photonic devices, the PMMgeometry allows the same device configuration to achieveboth propagation modes. Quickly switching between thestraight and bent waveguide behavior, which is easily en-abled by the controllable plasma current, opens up thepossibility of transistor-like optical switches.
B. Demultiplexer
Next, we present a frequency demultiplexer designed todistinguish between frequencies of ˜ ω = 1 and ˜ ω = 1 . m = 1 mode. This means we cannotuse the same mode integral as before. Instead we simplypenalize the simulated field intensity at the incorrect ex-its: (cid:82) | E | dl . With all this considered, the demultiplexerobjective is: L mp = (cid:18)(cid:90) E ω · E ∗ m =1 dl ω (cid:19) (cid:18)(cid:90) E ω . · E ∗ m =1 dl ω . (cid:19) − (cid:18)(cid:90) | E ω | dl ω . (cid:19) (cid:18)(cid:90) | E ω . | dl ω (cid:19) , where E ω is the simulated field for the ˜ ω = 1 source and E ω . is the simulated field for the ˜ ω = 1 . ω = 1 . ω p > ω src and there-fore (cid:15) rod <
0, the algorithm did not push rod permit-tivities deep into negative values. Perhaps bringing thesource frequencies closer together might allow more ef-fective usage of LSP excitations to differentiate the twofrequencies. Fig. 7 shows the evolution of the objectivefor the two devices. It appears as before that a suitablelocal maximum is obtained in both cases. Once again, thefinal distribution of plasma frequency among the PMMelements for the demultiplexers provides little suggestionof their intended function. ˜ ω = 1 ˜ ω = 1.1 TE ⃗ E w TM ⃗ B w (a)(b) (b.ii)(b.i) (b.iii)(a.ii)(a.i) (a.iii) ˜ ω = 1 FIG. 6. (a) TE FDFD simulation showing (a.i) the field mag-nitude | H z | (ˆ z out of the page) for ˜ ω = 1 and (a.ii) ˜ ω = 1 . ω = 1 where the maximum plasma frequency amongthe plasma elements is ˜ ω p = 1 .
17. (b) FDFD simulations and˜ ω = 1 permittivities for the TM case where the maximumplasma frequency is ˜ ω p = 1 . C. Low plasma frequency designs
Now, while these cases do exhibit compelling results,there is much to be said about their experimental re-alizability. Using the dimensions of plasma dischargetubes like those in refs. [22] and [23], the dimensional-ized lattice frequency is estimated to be ∼
20 GHz. Thus
FIG. 7. Evolution of the objective functions for (a) the TEdemultiplexer and (b) the TM demultiplexer. Once again, they-axes on these plots are effectively arbitrary since they arenormalized by the initial (epoch 0) value. these devices call for plasma frequencies as high as ∼ ∼ × m − .These conditions are only possible through pulsed opera-tion, which introduces complex transient behavior. In aneffort to examine how these devices may be realized withcurrently existing plasma sources, we consider the casewhere the lowest frequency source allowed by the waveg-uides (width 2 a yielding ˜ ω = 0 .
25) is used and the plasmafrequency is limited to 7 GHz (an attainable quasi steady-state operation condition). These limits are enforced viaa reparameterization of the optimization algorithm. Theresults for these waveguide cases along with the evolutionof the objective function are presented below in Figs. 8and 9 respectively.We observe that for both objective functionalities, theTM-polarization case exhibits very strong performancewhile the TE case struggles to match its higher frequencyalternatives. The reason for this is likely enhanced cou-pling to LSP modes. Since the plasma rods are now muchsmaller in radius than the wavelength of the source, LSPsare much more likely to appear. This is evident in the TEobjective training curves which are far more erratic thantheir TM counterparts. Since the algorithm adjusts thepermittivities of many rods at each step after calculatingthe gradients with respect to each parameter, the pres-ence of LSP excitation can cause unexpected spuriouslosses that cause the objective to abruptly drop. We seethe same performance discrepancies in the demultiplexer,shown in Fig. 10, and the demultiplexer objective, plot-ted in Fig. 11. Despite the losses in the TE-polarizationcases, these lower frequency simulations suggest that areconfigurable device operating in the TM mode couldbe constructed with experimental discharge tubes suchas those described in our prior studies [22, 23]. Morework needs to be done to identify objective functionali-ties that can reap the benefits of coupling into TE LSPs, TE ⃗ E w TM ⃗ B w (a)(b) (a.i) (a.ii)(b.ii)(b.i) (c.i)(d.i) (d.ii)(a.iv)(b.iv) TE ⃗ E w TM ⃗ B w (a)(b) (a.i) (a.ii)(b.ii)(b.i) (a.iii)(b.iii) FIG. 8. (a) TE FDFD simulation showing (a.i) the field magnitude | H z | (ˆ z out of the page) for an optimized straight waveguideand (a.ii) the relative permittivity domain that gives rise to this behavior, along with (a.iii) the field magnitude for the bentwaveguide and (a.iv) its permittivity domain. (b) FDFD simulations and domains for the TM polarization. Among all of thesecases, the highest plasma frequency is ˜ ω p = 1 .
35 which corresponds to ∼ .
75 GHz within existing experimental facilities.
TETM (a) (b) (a.i) (a.ii)(b.ii)(b.i)
Training EpochTraining Epoch
Straight
Training EpochTraining Epoch
Bent
FIG. 9. Evolution of the objective functions for (a) the low-frequency TE waveguides in the (a.i) straight waveguide caseand (a.ii) the bent waveguide case. The same is presented for(b) the TM waveguides. The y-axes on these plots are againarbitrary. perhaps with devices which require an extremely sensi-tive frequency response.
D. Optical logic gates
The success of these preliminary devices encouragesexperimentation with more complicated objectives anddevice geometries. Past work has established a connec-tion between the wave dynamics inherent to electromag-netic phenomena and the complicated computations per-formed by recurrent neural networks (RNNs) [27]. Thistype of analog computing appears to be the natural ap-plication of these PMM devices. The reconfigurability of ˜ ω = 0.25 ˜ ω = 0.27 TE ⃗ E w TM ⃗ B w (a)(b) (b.ii)(b.i) (b.iii)(a.ii)(a.i) (a.iii) ˜ ω = 0.25 FIG. 10. (a) TE FDFD simulation showing (a.i) the fieldmagnitude | H z | (ˆ z out of the page) for ˜ ω = 0 .
25 and (a.ii)˜ ω = 0 .
27, along with (a.iii) the relative permittivity of theplasma rods when ˜ ω = 0 .
25 where the maximum plasmafrequency among the plasma elements is ˜ ω p = 1 .
15 whichcorresponds to ∼ .
75 GHz within our existing experimentalfacilities. (b) FDFD simulations and ˜ ω = 0 .
25 permittivi-ties for the TM case where the maximum plasma frequencyis ˜ ω p = 1 .
43 or ∼ .
13 GHz experimentally. (a) (b)
Training Epoch TE Training Epoch TM FIG. 11. Evolution of the objective functions in the low-frequency case for (a) the TE demultiplexer and (b) the TMdemultiplexer. Once again, the y-axes on these plots are ar-bitrary. gaseous plasmas means that the complicated classifica-tion tasks which are the staple of RNNs can be readilytrained and implemented in-situ . As a natural depar-ture from traditional Von Neumann architectures, theseoptical computing platforms might yield significant ad-vantages in bandwidth and throughput in an age whereelectronic computing appears to be reaching physical lim-its. Here we explore the application of inverse design tothe fabrication of devices crucial to the development ofmore general optical computing platforms. Specifically,we consider the basic building block of traditional com-puting platforms: the electronic logic circuit. As physicalrealizations of boolean algebra, logic gates serve as thebackbone for the core arithmetic, storage, and processingoperations used in modern computing architectures. Op-tical logic gates could serve as the basic components ofoptical computing platforms [28], and have been realizedusing photonic crystal systems [29–31] as well as othertechniques like semiconductor optical amplifiers [32]. Wepropose the application of inverse design PMM devicesas a novel avenue for the fabrication of reconfigurableoptical logic devices.These devices represent the most complicated objec-tives considered thus far, particularly since the booleanarithmetic is an inherently nonlinear algebra. Becauseof this, we remove the plasma frequency constraints im- posed for the simpler devices and use a TM polarizedsource frequency of ˜ ω = 2 to allow for more complexpropagation structures. As with the demultiplexers, thesame device must be optimized to handle multiple casesdepending on the source. Figs. 12 and 13 illustrateour platform for the realization of boolean logic throughplasma arrays. Logical 1s and 0s are captured by thepresence or absence of a modal source at each of the twoinput data waveguides. To account for situations where alogical 1 is required while both databits are off (as is thecase for a NOR or NAND gate), a constant field source isprovided via the bottom waveguide. This in turn neces-sitates the addition of a ”ground” sink in the case whena logical 0 is required. The inclusion of this sink has theadded benefit of dissuading reflection back into the inputwaveguides, which would hinder any efforts to link con-nectives and perform more complicated computations.While the previous cases were run entirely by the algo-rithm, the complex nature of this problem required theaddition of a set of variable weights ( w - w and q - q )that could be manually tuned after each run. With addi-tional computing power, one could design another shellof optimization around the Ceviche algorithm to obtainthese weights automatically. The objectives used to trainthe permittivity domain for the gates were: L AND = w (cid:90) E Q ,P · E ∗ m =1 dl NAND − w (cid:90) | E Q ,P | dl AND + w (cid:90) E Q ,P · E ∗ m =1 dl NAND − w (cid:90) | E Q ,P | dl AND + w (cid:90) E Q ,P · E ∗ m =1 dl NAND − w (cid:90) | E Q ,P | dl AND + w (cid:90) E Q ,P · E ∗ m =1 dl AND − w (cid:90) | E Q ,P | dl NAND , L OR = q (cid:90) E Q ,P · E ∗ m =1 dl NOR − q (cid:90) | E Q ,P | dl OR + q (cid:90) E Q ,P · E ∗ m =1 dl OR − q (cid:90) | E Q ,P | dl NOR + q (cid:90) E Q ,P · E ∗ m =1 dl OR − q (cid:90) | E Q ,P | dl NOR + q (cid:90) E Q ,P · E ∗ m =1 dl OR − q (cid:90) | E Q ,P | dl NOR , where the weights are (5 , , , , , , . , .
5) for theAND case and (4 , , , . , , . , ,
1) for the OR case.The results were found to be highly dependent on thechoice of these weighting constants, suggesting that theshell optimization mentioned above might lead to betterperformance as opposed to using the lens maker’s algo-rithm. Experimenting with different configurations of theobjective determined that the simple linear combinationof the field integrals above was the easiest to coax intogood performance. Attempting to maximize the ratio be-tween field intensities at the output and ground waveg-uides resulted in an optimization domain with many sin-gularities, making convergence exceedingly difficult, andan objective consisting of a product of differences wouldyield the same objective for both the AND and OR gatesince the parity of these objectives is the same.The results shown in Figs. 12 and 13 are very com- pelling. In all cases, the field intensity at the correct out-put waveguide is substantially higher than at the incor-rect output. Though some operating cases exhibit sub-stantial loss, we understand that these losses are neces-sary when trying to achieve a nonlinear objective throughlinear dynamics. For example, the AND gate exhibits alarge amount of leakage in the upper right corner of thedevice. This odd mode of propagation (halfway betweenthe two exit guides) allows small changes in the sourcesto change which output waveguide sees a higher field in-tensity. Overall, these results show us that simple PMMdevice architectures can be reconfigured to achieve a widearray of highly non-trivial functions when inverse designis utilized. QP Q NOR P Q OR PConst. src. I n p u t Training Epoch
Q = 0, P = 0 Q = 0, P = 1 Q = 1, P = 0 Q = 1, P = 1 (a)(b) (c)
FIG. 12. (a) E-field magnitudes for the OR gate in each of theinput cases. (b) The evolution of the objective throughouttraining. (c) The relative permittivity domain along withlabels for the device input and exit waveguides where themaximum plasma frequency is ˜ ω p = 0 . QP Q NAND P Q AND PConst. src. I n p u t Training Epoch
Q = 0, P = 0 Q = 0, P = 1 Q = 1, P = 0 Q = 1, P = 1 (a)(b) (c)
FIG. 13. (a) E-field magnitudes for the AND gate in each ofthe input cases. (b) The evolution of the objective throughouttraining. (c) The relative permittivity domain along withlabels for the device input and exit waveguides where themaximum plasma frequency is ˜ ω p = 0 . IV. CONCLUSION
There are several potential opportunities for furtherwork based on the results discussed so far. The waveg-uide and demultiplexer devices designed here can be read-ily implemented and compared against simulated func-tionality using existing 2D PPC devices. Those parame-ters which are here taken as initial conditions ( r/a, ε rod )can themselves be used as training parameters. The lat-tice constant in particular was noted to have a large im-pact on the trade-off between transmission to the outputwaveguide and control of propagation in the design re-gion. Though in this study we chose the lattice constantto be commensurate with existing experimental equip-ment, more tightly-packed rods would likely lead to bet-ter performance. Consider also the prospect of a hexago-nal or a 3D PMM arrangement. The non-dimensionalized operation frequency ( ω in terms of c/a ) may also proveto be an important optimization parameter as the disper-sive physics which are the operational basis for these de-vices are highly frequency-dependent. In this study, asidefrom choosing the lowest frequency possible to reduce theplasma density requirements, other operating frequen-cies were determined without any substantial amount ofreasoning. Of particular interest in the future are the ef-fects of magnetization on the landscape of possible PMMdevices. When used in conjunction with circularly po-larized light, gaseous plasmas are capable of achieving ε r >
1, opening the door to even more complex logi-cal operations. While several tools exist for the simu-lation and inverse design of photonic devices [33], moredevelopment is required on the application of inverse de-sign to magnetized PMMs to incorporate such physics.Even these non-magnetized simulations can be furtherrefined by adding collisionality (loss) and non-uniformplasma density profiles to the simulation. If the objectivefunctionalities shown in this report can be achieved withthese two non-ideal factors considered, then experimentalrealization would likely prove to be much more attain-able. Finally, should experimental realization directlyfrom simulation results prove difficult, the device train-ing could actually be carried out in-situ , with plasmaelement currents adjusted while the objective (transmis-sion at a certain point, for example) is measured in real-time. In effect, the Ceviche forward-mode scheme couldbe carried out in the laboratory.In conclusion, we consider the application of inversedesign methods to the creation of highly optimized 2DPMM devices – including waveguides, demultiplexers,and novel photonic logic gates – for both TE and TMpropagation modes. Objective functions of varying com-plexities are considered depending on the desired device.Optimization is shown to take advantage of the richphysics inherent to PMMs including tunable reflection,refraction, and the presence of localized surface plasmonmodes. In the waveguiding and demultiplexer cases, sim-ulated devices are shown with experimentally realizableparameters with few necessary changes to existing exper-imental infrastructure. However, further work is needed,both in simulations and experimental development, toenable device development for more complex objectiveslike boolean logic as well as other devices with complexwave propagation structures.
ACKNOWLEDGMENTS
A.A. would like to thank the Tomkat Center for sup-port via the Energy Impact Fellowship. This research isalso partially supported by the Air Force Office of Scien-tific Research through a Multi-University Research Ini-tiative (MURI) with Dr. Mitat Birkan as Program Man-ager. J.A.R. acknowledges support by the U.S. Depart-ment of Energy, Office of Science, Office of Advanced Sci-entific Computing Research, Department of Energy Com-putational Science Graduate Fellowship under Award Number [DE-SC0019323]. [1] M. Minkov, I. A. Williamson, L. C. Andreani, D. Gerace,B. Lou, A. Y. Song, T. W. Hughes, and S. Fan, Inversedesign of photonic crystals through automatic differenti-ation, ACS Photonics , 1729 (2020).[2] T. W. Hughes, I. A. D. Williamson, M. Minkov, andS. Fan, Forward-Mode Differentiation of Maxwell’s Equa-tions, ACS Photonics , 3010 (2019).[3] L. Su, D. Vercruysse, J. Skarda, N. V. Sapra, J. A.Petykiewicz, and J. Vuˇckovi´c, Nanophotonic inverse de-sign with SPINS: Software architecture and practical con-siderations, Applied Physics Reviews , 11407 (2020).[4] M. Burger, S. J. Osher, and E. Yablonovitch, Inverseproblem techniques for the design of photonic crystals,IEICE transactions on electronics , 258 (2004).[5] P. I. Borel, A. Harpøth, L. H. Frandsen, M. Kristensen,P. Shi, J. S. Jensen, and O. Sigmund, Topology optimiza-tion and fabrication of photonic crystal structures, Opt.Express , 1996 (2004).[6] Z. Lin, A. Pick, M. Lonˇcar, and A. W. Rodriguez, En-hanced spontaneous emission at third-order dirac excep-tional points in inverse-designed photonic crystals, Phys.Rev. Lett. , 107402 (2016).[7] O. D. Miller, Photonic design: From fundamental so-lar cell physics to computational inverse design, arXivpreprint arXiv:1308.0212 (2013).[8] R. E. Christiansen, Z. Lin, C. Roques-Carmes,Y. Salamin, S. E. Kooi, J. D. Joannopoulos, M. Soljaˇci´c,and S. G. Johnson, Fullwave maxwell inverse design ofaxisymmetric, tunable, and multi-scale multi-wavelengthmetalenses, Opt. Express , 33854 (2020).[9] R. Pestourie, C. P´erez-Arancibia, Z. Lin, W. Shin, F. Ca-passo, and S. G. Johnson, Inverse design of large-areametasurfaces, Opt. Express , 33732 (2018).[10] H. Chung and O. D. Miller, Tunable metasurface in-verse design for 80% switching efficiencies and 144 ° angular deflection, ACS Photonics , 2236 (2020),https://doi.org/10.1021/acsphotonics.0c00787.[11] M. Meem, S. Banerji, C. Pies, T. Oberbiermann, A. Ma-jumder, B. Sensale-Rodriguez, and R. Menon, Large-area, high-numerical-aperture multi-level diffractive lensvia inverse design, Optica , 252 (2020).[12] A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz,T. M. Babinec, and J. Vuˇckovi´c, Inverse design anddemonstration of a compact and broadband on-chipwavelength demultiplexer, Nature Photonics , 374(2015).[13] D. Liu, L. H. Gabrielli, M. Lipson, and S. G. Johnson,Transformation inverse design, Optics express , 14223(2013).[14] T. W. Hughes, M. Minkov, I. A. Williamson, andS. Fan, Adjoint method and inverse design for nonlinearnanophotonic devices, ACS Photonics , 4781 (2018).[15] S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vuckovi´c,and A. W. Rodriguez, Inverse design in nanophotonics,Nature Photonics , 659 (2018).[16] R. E. Christiansen and O. Sigmund, Inverse design inphotonics by topology optimization: tutorial, JOSA B , 496 (2021). [17] J. Jiang, M.-K. Chen, and J. A. Fan, Deep neural net-works for the evaluation and design of photonic devices,ArXiv abs/2007.0 (2020).[18] N. M. Andrade, S. Hooten, S. A. Fortuna, K. Han,E. Yablonovitch, and M. C. Wu, Inverse design optimiza-tion for efficient coupling of an electrically injected opti-cal antenna-led to a single-mode waveguide, Opt. Express , 19802 (2019).[19] A. Iwai, F. Righetti, B. Wang, O. Sakai, and M. A. Cap-pelli, A tunable double negative device consisting of aplasma array and a negative-permeability metamaterial,Physics of Plasmas , 23511 (2020).[20] B. Wang, J. A. Rodr´ıguez, and M. A. Cappelli, 3D wood-pile structure tunable plasma photonic crystal, PlasmaSources Science and Technology , 02LT01 (2019).[21] B. Wang, J. Rodriguez, and M. Cappelli, 3D topolog-ical plasma photonic crystal with surface plasmon andFano-resonance modes, in APS Annual Gaseous Electron-ics Meeting Abstracts , APS Meeting Abstracts (2019) p.FT1.010.[22] B. Wang and M. A. Cappelli, Waveguiding and bendingmodes in a plasma photonic crystal bandgap device, AIPAdvances , 65015 (2016).[23] B. Wang and M. A. Cappelli, A plasma photonic crys-tal bandgap device, Applied Physics Letters , 161101(2016).[24] B. Wang, O. Miller, and M. Cappelli, A 3D hexagonal-packed photonic crystal with a tunable plasma-filled de-fect, in APS Annual Gaseous Electronics Meeting Ab-stracts , APS Meeting Abstracts (2019) p. MW1.079.[25] F. Righetti, B. Wang, and M. Cappelli, Enhanced at-tenuation due to lattice resonances in a two-dimensionalplasma photonic crystal, Physics of Plasmas , 124502(2018).[26] D. P. Kingma and J. Ba, Adam: A method for stochasticoptimization (2014), arXiv:1412.6980 [cs.LG].[27] T. W. Hughes, I. A. D. Williamson, M. Minkov, andS. Fan, Wave physics as an analog recurrent neural net-work, Science Advances , eaay6946 (2019).[28] L. De Marinis, M. Cococcioni, P. Castoldi, and N. Andri-olli, Photonic Neural Networks: A Survey, IEEE Access , 175827 (2019).[29] Y. Fu, X. Hu, and Q. Gong, Silicon photonic crystal all-optical logic gates, Physics Letters A , 329 (2013).[30] P. Rani, Y. Kalra, and R. Sinha, Realization of and gatein y shaped photonic crystal waveguide, Optics Commu-nications , 227 (2013).[31] P. Rani, Y. Kalra, and R. Sinha, Design and analy-sis of polarization independent all-optical logic gates insilicon-on-insulator photonic crystal, Optics Communica-tions , 148 (2016).[32] G. Berrettini, A. Simi, A. Malacarne, A. Bogoni, andL. Poti, Ultrafast integrable and reconfigurable xnor,and, nor, and not photonic logic gate, IEEE PhotonicsTechnology Letters , 917 (2006).[33] J. Jiang, R. Lupoiu, E. W. Wang, D. Sell, J. P. Hugonin,P. Lalanne, and J. A. Fan, MetaNet: a new paradigmfor data sharing in photonics research, Opt. Express ,0