Schrödinger's Red Beyond 65,000 Pixel‐Per‐Inch by Multipolar Interaction in Freeform Meta‐Atom through Efficient Neural Optimizer

Abstract Freeform nanostructures have the potential to support complex resonances and their interactions, which are crucial for achieving desired spectral responses. However, the design optimization of such structures is nontrivial and computationally intensive. Furthermore, the current “black box” design approaches for freeform nanostructures often neglect the underlying physics. Here, a hybrid data‐efficient neural optimizer for resonant nanostructures by combining a reinforcement learning algorithm and Powell's local optimization technique is presented. As a case study, silicon nanostructures with a highly‐saturated red color are designed and experimentally demonstrated. Specifically, color coordinates of (0.677, 0.304) in the International Commission on Illumination (CIE) chromaticity diagram – close to the ideal Schrödinger's red, with polarization independence, high reflectance (>85%), and a large viewing angle (i.e., up to ± 25°) is achieved. The remarkable performance is attributed to underlying generalized multipolar interferences within each nanostructure rather than the collective array effects. Based on that, pixel size down to ≈400 nm, corresponding to a printing resolution of 65000 pixels per inch is demonstrated. Moreover, the proposed design model requires only ≈300 iterations to effectively search a thirteen‐dimensional (13D) design space – an order of magnitude more efficient than the previously reported approaches. The work significantly extends the free‐form optical design toolbox for high‐performance flat‐optical components and metadevices.


Supplementary Text 1.
Parameters of the VAE

Figure. S1. The architecture of the VAE
The detailed parameter of the VAE network is shown in Figure S1.The input of the encoder is a 64 by 64 binary image of 0 or 1 representing the unit cell design.The encoder consists of 4 layers of 2D convolutional (conv.)layers and 2 layers of linear layers.The filter size, the number of channels, the padding size, and the stride are represented by "F", "c", "p", "s" respectively.The letter "I" and "O" represents the input and output size of the linear layer.The input of the decoder is a 10-dimensional array.The input is passed through two linear layers and 4 transposed convolutional layers to reconstruct the original image with 64 by 64 pixels.
During the training, We constructed a geometrical library with 20000 unit cell designs as the training data.The VAE can pick up the pattern and produce unit cell designs with the same C4v symmetry.After training, the decoder can be used to map array inputs into geometrical shapes.

The global optimization module
Here we give a detailed description of the global optimization module of the proposed algorithm.As shown in Figure S2, The global optimization module consists of two submodules, namely the MCTS module and the local Bayesian module.The MCTS takes in all the samples taken so far and runs a search algorithm to find the best partition of the samples so that the selected samples have the highest average FOM : . The selected samples represent the most promising area to search.After the partition, a local Bayesian process provides the next sample to evaluate within the promising area.
The MCTS follows the standard procedure: selection, expansion, rollout, and backpropagation.In each iteration, the MCTS selects a node based on the upper confidence bound (UCB), which is a standard way of choosing the best node.In the expansion process, each node is split into a predefined number of children nodes by using the K-means clustering algorithm implemented by the open-source scikit-learn library [1,2] .Each node represents a different way of splitting the existing samples.During the rollout, the splitting algorithm is called progressively until the number of samples in a node reaches a certain limit.are backpropagated to the selected node to select the best partition strategy.Since all the data points are not clustered in high dimensions, different random states are used to get different divisions.The number of samples is used as a threshold to determine whether a node is splittable.It also controls the ratio between exploration and exploitation.If the threshold is big, the search area is also big.In contrast, if this number is small, the area is contracted to a smaller area.After the MCTS picks out a set of samples.An SVM can be trained to determine the boundary between the selected sample and the rest of the sample.Bayesian optimization process is used subsequently to suggest a new sample for the next evaluation in the promising area.We first train a Gaussian process regressor with the full set of samples to build up a surrogate model, after which we maximize the acquisition function to give the next sample to evaluate.To sample for the acquisition function, the pre-trained SVM is used as a reference for rejection samples to get samples distributed in the picked volume.In this way, the BO is confined to a promising region and thus improves the optimization efficiency.The acquisition function is either expected improvement or probability of improvement acquisition function.The main control takes the suggested sample and calls the numerical simulation module to evaluate the FOM.Afterwards, the FOM and the sample are added to the existing sample pool and a new iteration of optimization is started.This global optimization iterates until the FOM or the number of iterations reaches a threshold value.Figure .S3 compares our algorithm with previous reports [3][4][5][6][7][8][9] using RL and heuristics optimization algorithms.The typical number of function calls are in the 10 3 range with a much smaller design parameter size for previous reports, as opposed to around 300 function call in 13-dimensional space in this work.The high efficiency in data usage and convergence speed is attributed to the combination of the efficient global algorithm and the local optimization algorithm.

Multipole expansion
The multipole expansion is carried out in Cartesian coordinate based on the displacement currents  = iω 0 (  − 1) , where  is the angular frequency,  0 is the permittivity in vacuum, and   is the complex relative permittivity.The multipole moments are calculated as follows: Electric dipole moment:  = The total radiation power of each multipole moment is: Here,  is the speed of light.α, β, γ = x, y, z are the Cartesian coordinates.

Size dependent reflection
Figure S4(a) depicts the reflection spectra for different array size.The suppression of reflection at wavelength below 600 nm is consistent for all the array size.The reflection peak in the red band due to the response is also present for the smallest array simulated.The dropped reflection when the array gets smaller is due to the limitation of the monitor size during the simulation, which can only capture part of the reflected light.As the array size get smaller, the scattered field resembles more spherical wave, and only a small portion of the scattered field can be captured by the field monitor.However, it is clear from the measured optical microscope data, even a single unit cell shows high saturation red color.Figure S4(b) shows the corresponding CIE colors for various array size.As the array size gets bigger, the coordinates move towards the red region.It is estimated that a color saturation comparable to Cadmium red can be obtained with an array size of 6 by 6-unit cells.Figure S4(c-d) shows the multipole decomposition of the scattering cross section of a single scatterer.The dashed line in Figure S4(c) is the scattering cross section, it shows a pronounced peak at 600-700 nm range, which coincides with the enhancement of reflection in the array structures.Moreover, we also observe the cancellation of the quadrupole components in the lower wavelength range, similar to that of the array.This indicates that the single unit cell is able to reach certain cancellation and the array effect is less important in our structure.

Comparison of the optimized structure with sub-optimized structures
In this section, we show the comparison of the optimized structure with various unoptimized structure to illustrate the role of the multipolar cancellation.In all the structures shown in Figure S5, the volume is conserved while only the arrangement of the material is altered.In all these structures, the reflection peak at the red region remains high.However, there is substantial reflection in the shorter wavelength region for the unoptimized structures, as shown in Figure S5(a, b).The reflection at shorter wavelength results in the CIE coordinate moving toward the blue and yellow region and therefore lower red saturation.The effective cancellation of the reflection at lower wavelength is due to the minute variation of the multipole moments and multipole phase as shown in Figure S5, last two columns to the right.The electric field intensity distribution is very similar at longer wavelength (680 nm) for the optimized structures and a slightly altered one (Figure S6).As the wavelength get shorter, the difference become more significant, which is also consistent with the previous observation.This further confirms the variation of the geometries can effectively modulate the optical response in a broadband wavelength.

Figure. S2 .
Figure.S2.The flowchart for the global optimization module

Figure. S4 .
Figure.S4.(a) The simulated reflection spectra for different array size.(b) The green dots show the CIE coordinates for the array of different size.The Cadmium red (C-red), are also shown for reference.The dashed line is the boundary of sRGB standard.And the green arrow indicates the increasing array size.(c) The multipole components and the scattering cross section of a single isolated unit cell.(d) The multipole phase of the scattering cross section of a single isolated cell.

Figure. S5 .
Figure.S5.Comparison of different geometries.(a) Circular geometry with the same volume as the optimized geometry.(b) Complex geometry with slight variation.(c) The optimized geometry.The rows in each panel are simulated reflection spectra, CIE color coordinate, multipole moments and multipole phase respectively.The electric field intensity distribution is very similar at longer wavelength (680 nm) for the optimized structures and a slightly altered one (FigureS6).As the wavelength get shorter, the difference become more significant, which is also consistent with the previous observation.This further confirms the variation of the geometries can effectively modulate the optical response in a broadband wavelength.

Figure. S6 .
Figure.S6.Comparison of the E field distribution at different wavelength between (a) suboptimized structure and (b) optimized structure.FigureS7 (a, d) shows the comparison of the SEM images of the optimized structure and the sub-optimized one, respectively.The differences in the geometry are subtle, but the optical response is quite different, especially at shorter wavelength.FigureS7(e) shows the measured

Figure. S7 .
Figure.S7.SEM image, reflection spectrum, and CIE coordinate of the optimized structure (ac) and the sub-optimized structure (d-f), respectively.