On the design of random metasurface based devices

Metasurfaces are generally designed by placing scatterers in periodic or pseudo-periodic grids. We propose and discuss design rules for functional metasurfaces with randomly placed anisotropic elements that randomly sample a well-defined phase function. By analyzing the focusing performance of random metasurface lenses as a function of their density and the density of the phase-maps used to design them, we find that the performance of 1D metasurfaces is mostly governed by their density while 2D metasurfaces strongly depend on both the density and the near-field coupling configuration of the surface. The proposed approach is used to design all-polarization random metalenses at near infrared frequencies. Challenges, as well as opportunities of random metasurfaces compared to periodic ones are discussed. Our results pave the way to new approaches in the design of nanophotonic structures and devices from lenses to solar energy concentrators.

Recent advances on the control of light in complex media 36 have motivated the study of random or disordered metasurfaces for specific applications such as decreasing the radar cross-section [37][38][39][40] , improving SERS enhancement 41 , reducing laser coherence 42 , designing wide band-gaps 43 , or increasing light-matter interaction and the absorption of solar cells 44,45 . One of the advantages of random media is the very high number of degrees of freedom that they support and which can be harnessed to control waves on scales smaller than the wavelength, or to multiplex more information for communication 36,46 . This has recently led to the design of random metasurfaces for wave front shaping [47][48][49] . However, the design of such devices still remains elusive due to the disordered distances between neighboring elements, the near-field coupling, and variations of the local density of elements. Some theoretical approaches can address the homogenization problem of homogeneous random polarizability materials in periodic arrays of resonators 50 , or for identical polarizabilities in disordered arrays of scatterers 51,52 .
The relation providing the phase-shift of the elements constituting a metasurface as a function of their dimensions is determined either analytically, when possible, or with numerical simulations for a single element or for periodic arrays of identical elements. However, metasurfaces are generally made of elements of different sizes to provide a phase-shift that varies spatially. Hence, the previous approaches may fail 53 as near-field coupling introduces errors in the phase-shifts provided by the elements. Important questions are thus whether this periodic arrangement is always the best solution and whether it is possible to design functional metasurfaces within a random framework with general guidelines.
While random and disordered metasurfaces can be complex to design, they also have potential advantages. For instance, the random design process optimizes the area of the metasurface. In a periodic metasurface, small and large elements have the same footprint. On the contrary, in a random metasurface, the random design or the pseudo random algorithm finds more easily a spot for a small element than for a larger one. This optimizes the local density and the footprint of the elements. Furthermore, the absence of periodicity eliminates any spurious diffraction orders that arise from large periods 54 , due for example to large resonators made of low index materials. The circular symmetry of the elements is also statistically restored by the randomness 55 , which enables the design of polarization independent metalenses with anisotropic elements. This contrasts with the current works implementing polarization independent lenses using circular or fourfold symmetric cross-section elements 19,31,[56][57][58] .
Here, using anisotropic gold nano-elements as resonators, we design random metasurface lenses at the wavelength of λ 0 = 1.5 μm. Such metasurfaces are designed using randomly sampled elements, which lengths are computed in order to provide the required phase law at each of the sampled position. To establish general rules and guidelines for such designs, the resonators are first considered in a periodic framework to numerically obtain their phase maps φ(l), i.e. the phase-shift provided by the element as a function of a geometrical parameter-the length in our case-for different periods of the array. These periods, corresponding to a density of the phase-map, are then used as references from which a resonator can be chosen. Then, using different phase-maps, one-and two dimensional metasurfaces of periodic or random elements of various densities can be constructed. The performances of the designs are then discussed.

Results
Phase-maps and unit cell simulations. To design random metalenses such as the one represented in Fig. 1(a), phase-shifting elements need to be used. Figure 1(b) presents a possible implementation using anisotropic gold nano-resonators. The elements has a width of 50 nm, a height of 40 nm, and the length is tuned between 150 nm and 500 nm. Gold is modeled using a Drude model with a plasma frequency ω p = 1.367 × 10 16 rad/s and a collision frequency ω c = 6.478 × 10 13 rad/s 59 . This plasmonic particle is supported by a dielectric SU8 spacer with a refractive index n SU8 = 1.59, optimized to a thickness of 70 nm, on top of a metallic ground plane (more details in supplementary information).
The building block of the metasurfaces is first investigated using in-plane periodic boundary conditions, as shown in Fig. 1(b). The period in the direction of the width of the element (period p y in the y direction) will be swept but is kept to 100 nm in Fig. 1. The period along the longer dimension of the element (p x in the x direction) is set to 900 nm and is kept constant all along the paper. Using the frequency domain solver of the commercial software CST, the complex reflection coefficient of our structure is computed. The illuminating plane wave has a frequency varying from 50 THz to 350 THz and is polarized along x (i.e. the long axis of the particle). The particle is transparent to the orthogonal polarization (see supplementary information). Varying the length of the particle from 150 nm to 500 nm shifts its fundamental resonance frequency as shown in Fig. 1(c). The phase of the wave reflected by the particle and the metallic plane can thus be controlled. It is worth noting that there is no transmission because of the metallic ground plane and the absorption is 1-R. Figure 1(d) shows the phase shift of an element around 200 THz (λ 0 = 1500 nm) for different particle lengths. The shortest element is taken as phase reference. Figure 2(a) shows the phase shift as a function of the length of the elements for different periods p y at 200 THz. For a single resonance, the complete 2π phase shift is only obtainable asymptotically far away from the resonance. The SU8 spacer thickness sets the quality factor Q of the resonances which in turn, controls the maximum value of the phase-shift-the thicker the SU8 layer, the lower the Q factor and the smaller the maximum phase-shift. However, the higher the Q factor, the sharper the slope of the reflected phase. Hence, a compromise has to be made between the maximum value of the phase shift and the slope of the reflected phase around 200 THz ( Fig. 1(d)). A very steep change of phase introduces phase errors 22,53 . A thickness of 70 nm (the red curve on Fig. 2(a)) appears to be a good compromise as the difference between 2π and the maximum phase-shift is smaller than 37°.
The phase-shift required to design a parabolic lens in reflection or concentrator with a focal length f is given by the parabolic law 22,27 : The phase-shift required to design a metalens of 30 μm width with a focal spot of 20 μm is shown in Fig. 2(b). Knowing the phase-shift required at any position of the metasurface ( Fig. 2(b)) and the phase shift induced at reflection on a periodic array as a function of the element length ( Fig. 2(a)) (reference phase-map), leads to choose the length of the elements as a function of their position on the metasurface. Figure 2(c) presents the length required at a given position to realize the phase shift plotted in Fig. 2(b), for different phase-maps represented by different periods p y (from 150 nm to 500 nm) in the periodic array. Changing the period shifts the phase of the field reflected by an array of identical elements. This originates from two reasons: the near-field coupling that becomes stronger as the distance between the elements is decreased, and the density of elements itself in the limit of negligible near-field coupling. Indeed, the denser the array, the more field will be phase-shifted by the elements compared to the field which is only reflected by the ground plane. The total reflected field which is the sum of the field reflected by the mirror and the field scattered by the elements has therefore different phases for different densities. The two insets of Fig. 2(c) show the "z" component of the electric field in the near field of the elements of 300 nm length and 150 nm at their resonant frequency (respectively 200 and 300 THz).
One-dimensional random metalenses. For a periodic metasurface, on one hand, choosing a period p y of a phase-map sets the density of elements per unit area of the periodic array: the relation between the period and the density is: ρ = p p 1/ x y . On the other hand, the main questions that arise are which density optimizes the focusing of a random metalens, and which phase-map should be chosen to design a metasurface at this density. A naïve response would be to select the phase-map with the same density as the random metasurface to be designed, but, as will be seen later, the response is not straightforward. These questions are all the more important that random metasurfaces have fluctuations of the near-field coupling that may affect the efficiency if the density becomes too high.
We simulated 16 different random metasurfaces corresponding to four densities of phase-maps (four periods in the y direction (p y ) of 100 nm, 150 nm, 250 nm and 500 nm) and four equivalent densities in the random metasurface of 25, 17, 10 and 5 elements per squared wavelength. On a matrix with the rows representing the density of the phase-maps and the columns representing the density of the random metasurface, elements on the diagonal are thus those for which the two densities are equal. Figure 3 represents 16 one-dimensional metasurfaces with a width (along "y") of 30 μm and a focal length (along "z") of 20 μm. Elements are set along "x" perpendicularly to the 1D metasurface ( Fig. 3(a), out of plane of the Fig. 3(b)). Periodic boundary conditions define the out of plane direction, with a period of 900 nm.
Our algorithm to design a random metasurface corresponds to random loose packing 60 and consists in the following steps. First, we randomly select a position (for 2D metasurfaces, we also choose a random orientation for the element). For 1D metasurfaces, all elements are parallel to their longest length. We compute, using the phase maps of Fig. 2(a-c) the length that a particle at this position should have. We then check if it is overlapping or is too close to previously placed elements. A minimum distance between the surfaces of the elements of 10 nm is set in order to put the elements as close as possible to get the maximum density but this value can be varied. If the minimum distance between elements is too small or if they overlap, we remove it and select another random position. If they do not overlap, we approve the change and move to the next particle. The process is repeated until we manage to place a defined number of elements (from 300 for a density of 25 λ 0 −2 to 60 for a density of 5λ 0 −2 ) or until we have failed to place a given element. This would mean that the maximal density of the random metasurface has been reached. Using CST simulations in time domain for better computational efficiency, we computed the density of energy of the reflected field normalized by the density of energy of the incident plane wave for the sixteen 1D random metasurfaces. Each design is randomly repeated, simulated ten times, and averaged to ensure that the results are statistically meaningful. The average results for the 16 (4 by 4) metasurfaces are displayed on Fig. 3. Unaveraged results are very similar as the maximum value of the standard deviation is found to be about 10% of the average value (see supplementary information). Figure 3 presents the density of reflected energy normalized to the density of incident energy for the 16 metasurfaces. Metasurfaces on the diagonal of the figure and below the diagonal have better focusing performance, i.e. higher energy at the focal spot, than those above the diagonal. This impression is confirmed by the measured values in Table 1. This figure and its table show that for the designed 1D random metasurfaces, the density of elements in the metasurface plays a more important role than the density of the particular phase-map chosen to design the metasurface. This is an interesting conclusion as near-field coupling fluctuation in the random system, when the phase-map is no longer representative of the metasurface, does not seem to alter the focusing ability of the metasurface. Two-dimensional random metalenses. We now consider 2D metalenses. As in the 1D case, metasurfaces are constructed using the same phase-maps corresponding periods p y of the reference of 100 nm, 150 nm, 250 and 500 nm. The size of the 2D metasurfaces is chosen to be 10 by 10 μm, with a focal length of 10 μm to limit the computational volume. The phase shift provided by the elements is now a function of x and y. The nominal frequency is still 200 THz (λ 0 = 1.5 μm). We again simulate 10 sets of 16 metasurfaces with densities from 25 λ 0 −2 to 5 λ 0 −2 and present the averaged results. The process to design the random metalenses is the same as for the 1D case, but with randomly placed and oriented elements as shown in Fig. 1(a). A minimum distance between the surfaces of the elements of 10 nm between adjacent elements is enforced and makes the structures realistic to fabricate.
The focusing results, i.e., the reflected density of energy in the plane y = 0, which corresponds to a cross section of the central volume, are shown in Fig. 4, while Table 2 displays the main characteristics of the focal spots. Metasurfaces produce well defined focal spots of 1.9 μm by 8 μm, i.e. slightly superior to Abbe's limit λ 0 /2NA = 1.6 μm and λ 0 /NA 2 = 6 μm. Very interestingly, the best results are not achieved anymore for metasurfaces below the diagonal. This contrasts with the 1D case, where metasurfaces of the lower-left part of the figure, i.e. which are designed for equal or higher densities than their corresponding phase-maps, provided good results. Here the three best results are obtained with the densest phase-map with p y = 100 nm. With this phase-map density, even the random metasurface with a density of 5/λ 0 2 leads to a visible focusing spot. The density of the phase-maps thus seems to play a more important role than the density of the random metasurface itself even if the two parameters obviously play a role. The phase shift in denser phase-maps accounts for a more important near-field coupling between elements. Near-field coupling thus seems to play a more important role in the designed 2D random metasurfaces than in the design of 1D metasurfaces. In 2D metasurfaces the elements are randomly oriented whereas they were all aligned in the 1D random metasurfaces. Near-field interactions in 2D metasurfaces are different from what is modeled in the phase-maps and thus raises fundamental questions on the optimal design of random metasurfaces. As noted in a recent paper 53 , near-field interactions between adjacent elements in a gradient metasurface do not account for the fact that elements are different, because phase-maps use identical elements. This problem is exacerbated in metasurfaces in which the cross-talk between elements is not negligible. While ideas to address this problem and improve the efficiency were proposed for gradient metasurfaces using periodic grids 53 , the question still remains open for random structures that thus operate away from their optimum.
It is worth underlining that near-field coupling is present in both 1D and 2D metasurfaces. However, in 1D metasurfaces, near-field coupling is well accounted for by phase-maps from simulations with periodic and aligned structures. Randomly oriented elements in 2D metasurfaces lead to near-field coupling that deviates more from the phase-maps simulations. Even if 2D near-field coupling is in average weaker than 1D near-field coupling (the coupling being stronger for collinear dipoles), the use of the proposed phase-maps leads to a stronger dependence of the 2D metasurface on the scatterers configuration.

Discussion
In the previous section, we have proposed a strategy to design 2D random metasurface lenses. A 2D random metasurface using anisotropic elements such as rectangle elements are expected to be polarization independent, an important property for many applications. The theoretical focusing power of the 2D random metasurfaces is expected to be half of the focusing power of the corresponding periodic 2D metasurface using the same anisotropic element. However, there is not a one to one map between the periodic and the random structure of the same density as near-field coupling between elements in the two surfaces is different. We present in supplementary information the polarization dependent periodic metasurfaces with the same elements. Figure 5 shows the cross section of energy at the focal spot at a distance of 10 μm from the metasurfaces for two incident polarizations for 2D random metasurfaces and for the metasurface with a periodic grid. The random 2D metalens is polarization independent as expected while the metasurface with periodic grid and aligned elements is not. However, as previously discussed, the power at the focal spot in Fig. 5(b) is smaller than half the power at the focal spot in Fig. 5(a) and this stems from the differences in near-field interaction configurations.
To further optimize the efficiency of random metasurface, a possible solution would consist in refining the references and use phase-maps representative of the configuration of randomness and accounting for the near-field coupling fluctuation in the metasurface. A local phase method that optimizes the length of elements for quasi periodic metasurfaces has been proven to be effective 53 , but it would be very computationally intensive for 2D random metasurfaces. Designing metasurfaces with isotropic elements in the unit cell also leads to polarization independence. Whether randomizing such systems can lead to better performance is still an open question. Indeed, the approach presented here could be improved by computing phase-maps of the random structures with randomly placed and oriented identical elements. Finally, an interesting feature of random metalenses can be seen in Fig. 1(a), the density of elements is not homogeneous. At locations where elements are larger, the density is smaller, while at positions where elements  are shorter, the density is higher. Such a distribution can be expected to compensate the fact that smaller elements have a smaller scattering cross-section. Hence, using disordered metalens provides additional degrees of freedom in the design of devices. Figure 5(b) presents the distribution of the length of the elements for a random and a periodic metasurface lens of 100 by 100 μm with a focal length of 250 μm and a density of 25 elements per squared wavelength (111,000 elements on the surface). We can see that a simple random loose packing 60 algorithm favors smaller elements. Such feature may be engineered by tailoring the algorithm used to design the metasurface, for instance using close packing algorithms or hyperuniform media [61][62][63][64] . Such algorithms may also be used to explore higher densities.

Conclusion
We proposed a method to design 1D and 2D random metasurface lenses. Using extensive numerical simulations, we demonstrated successful focusing by 1D and 2D random metasurfaces. By implementing random metalenses of various densities using phase-maps of same density (but periodic), we found that the main metric affecting the performance of random 1D metasurfaces is the density of the metasurface, while, in 2D random metasurfaces, the density of the phase-maps or the near-field coupling between elements seems to play a more important role than the density of elements in the metasurface itself. Randomness statistically restores the circular symmetry of the devices and enables polarization independent lenses. We have also demonstrated that random metasurfaces contain a larger number of small scatterers than their periodic counterpart and this may favor higher intensity at the focus if the optimal near-field couplings between random structures is obtained. Further investigations need to be performed to understand the role of the orientation disorder and the strength of the near-field coupling to optimize 2D random metalenses. Our results pave the way to the design of random metasurfaces for devices as diverse as lenses and concentrators. We also believe that random metasurfaces may overcome limitations on the diffraction efficiency of periodic systems, especially for dielectric metasurfaces that are made with larger elements. Random structures are also more amenable to self-assembly fabrication for large scale systems.

Methods
All simulations are performed using the commercially available software CST, version 2017. 1D and 2D metasurface simulations are performed in the time domain for better efficiency. Each layer of materials has an optimized maximum mesh size. It is λ/10 in air (the domain of the propagation), 15 nm in the gold elements and 30 nm in the SU8 spacer layer. This leads to about 35,000 mesh cells in unit cell simulations, about 16 million for 1D metasurfaces and up to 60 million for 2D metasurfaces. Convergence is ensured by changing the maximum size of the cells to 10 or 20 nm and making sure the results are identical, and, by making sure that energy in the computational volume is monotonously decreasing after the plane wave excitation.