††thanks: Contributed equally††thanks: Contributed equally

# Spatial wavefront shaping with a multipolar-resonant metasurface for
structured illumination microscopy

Tamal Roy  Department of Mechanical and Process Engineering, ETH Zurich,
Zurich, Switzerland  Peter T. Brown  Center for Biological Physics and
Department of Physics, Arizona State University, Tempe, AZ, USA  Douglas P.
Shepherd<EMAIL_ADDRESS>Center for Biological Physics and
Department of Physics, Arizona State University, Tempe, AZ, USA  Lisa V.
Poulikakos<EMAIL_ADDRESS>Department of Mechanical and Aerospace
Engineering, UC San Diego, La Jolla, CA, USA

###### Abstract

Structured illumination microscopy (SIM) achieves superresolution in
fluorescence imaging through patterned illumination and computational image
reconstruction, yet current methods require bulky, costly modulation optics
and high-precision optical alignment. This work demonstrates how nano-optical
metasurfaces, rationally designed to tailor the optical wavefront at sub-
wavelength dimensions, hold great potential as ultrathin, single-surface, all-
optical wavefront modulators for SIM. We computationally demonstrate this
principle with a multipolar-resonant metasurface composed of silicon
nanostructures which generate versatile optical wavefronts in the far field
upon variation of the polarization or angle of incident light. Algorithmic
optimization is performed to identify the seven most suitable illumination
patterns for SIM generated by the metasurface based on three key criteria. We
find that multipolar-resonant metasurface SIM (mrm-SIM) achieves resolution
comparable to conventional methods by applying the seven optimal metasurface-
generated wavefronts to simulated fluorescent objects and reconstructing the
objects using proximal gradient descent. The work presented here paves the way
for a metasurface-enabled experimental simplification of structured
illumination microscopy.

## Introduction

Structured illumination microscopy (SIM) is a fluorescence superresolution
microscopy technique which can reach twice the diffraction limit by exploiting
patterned excitation light combined with spatial frequency downmixing that
occurs in incoherent imaging [1]. SIM is increasingly popular because it is
compatible with live-cell imaging, does not require specific fluorophore
photophysics, and can capture images at rates approaching standard widefield
microscopy [2, 3, 4]. Most modern SIM instruments rely on illumination
patterning optics, typically diffraction gratings or spatial light modulators
(SLM’s), placed in a conjugate imaging plane. Consequently, these experiments
need precise optical alignment [5, 6], require considerable optical expertise,
and are expensive to build.

Recently, a variety of approaches have emerged to simplify, democratize, and
improve the speed, stability, and efficiency of SIM experiments. These include
fiber-based [7, 8], digital micromirror device (DMD)-based [9, 10, 11], and
photonic chip-based [12, 13] SIM schemes. Such efforts compliment parallel
computational imaging approaches to limit reconstruction artifacts and achieve
high-quality imaging at low signal-to-noise ratios [14, 15]. While each of
these approaches is promising, none fully solves the problems discussed above.
For example, fiber-based approaches simplify alignment but require expensive
or complex methods to phase shifting such as Pockels cells or thermal phase-
shifters. DMD-based approaches significantly reduce the cost of SIM setups but
the diffraction physics of the DMD introduces new complexity and places
constraints on multicolor operation. Photonic chip approaches offer the most
compact and simplified systems, but are either limited to near-field
excitations or have only been demonstrated for 1D resolution enhancement.

Progress in nanoscale optics has enabled emerging technologies for generating
structured light, which offer solutions to the abovementioned challenges.
Nano-optical metasurfaces are composed of sub-wavelength-periodic optical
elements that control light–matter interactions at the nanoscale. The high-
spatial-frequency evanescent near fields generated by plasmonic metasurfaces –
those composed of appropriate metallic nanostructures [16] – have been applied
to SIM to enhance spatial resolution [17, 18, 19, 20]. However, such
resolution enhancement can only be achieved at nm-scale proximity to the
metasurface, thus necessitating close contact to the sample of interest and
hindering applications such as depth-dependent confocal imaging. Moreover,
this method still requires complex piezoelectric stages to achieve
illumination patterns, maintaining the experimental complexity of conventional
SIM experiments. Beyond near-field effects, metasurfaces enable versatile
avenues for light manipulation in the far field. [21] Metasurfaces have been
rationally designed for full phase and polarization control of incident light,
[22] leading to numerous applications encompassing flat optical elements,[23,
24, 25, 26] beam steering,[27, 28, 29] holography,[30, 31, 32] or optical
computation,[33, 34, 35] as has been reviewed in prior work. [36, 37, 38, 39,
40, 41]

Notably, metasurfaces have unique capabilities to control the optical
wavefront at subwavelength dimensions through engineered manipulation of the
phase of incident light [22, 42]. However, the potential of metasurfaces to
act as ultrathin, single-surface, all-optical, generators of tunable SIM
illumination patterns, remains largely unexplored. Wavefront-shaping
capabilities of metasurfaces present an exciting opportunity to drastically
simplify the experimental configuration of SIM by leveraging the ability of
metasurfaces to generate far-field patterns in the optical wavefront.
Moreover, all-optical tunability of the metasurface-generated wavefronts are
enabled by tailoring the interaction of the metasurface with, e.g., varying
polarization states or incident angles of the illumination source.
Superresolution images can then be obtained by tailored computational
reconstruction algorithms.

In this work we introduce multipolar-resonant metasurface SIM (mrm-SIM), an
approach which retains the compact system found in photonic chip-based
approaches but uses far-field patterns and achieves 2D resolution enhancement.
mrm-SIM generates structured light by resonant interactions with nano-optical
structures and allows the illumination patterns to be tuned by changing the
incident polarization and angle. With mrm-SIM we achieve resolutions up to 240
for 1004 excitations, comparable to those achieved with conventional hexagonal
SIM. By embedding such a device in a glass coverslip we propose a nearly
alignment-free SIM approach. Since mrm-SIM uses far-field patterns, we
anticipate it can be extended to 3D SIM by appropriate pattern engineering.
Furthermore, unlike most other SIM approaches the achievable field-of-view
(FOV) in this approach is limited by the size of the metasurface instead of
the FOV of the detection system. Therefore, we expect mrm-SIM to be more
scaleable and achieve larger spatial-bandwidth products than other approaches.
This research paves the way toward significantly simplifying the experimental
configuration required for SIM and the democratization and broad dissemination
of superresolution microscopy techniques.

## Results and discussion

Figure 1:  Nanoresonator characterization.(a) Nanoresonator design and
geometry composed of a Si nanoresonator on a SiO2 substrate in air. The
nanoresonator has a height h = 715 and an elliptical cross sectional area with
width w = 200 and length l = 400 and is arranged periodically in an
equilateral rhombohedral lattice with side length L = 670. (b) Electric (left)
and magnetic (right) field norms of the nanoresonator at resonant wavelengths
of $\lambda$ = 1004 (top) and $\lambda$ = 1194 (bottom), indicating both
electric and magnetic multipolar resonant behavior. The incident plane wave
intensity is 1000/ ^2. (c-e) Transmittance spectra for the nanoresonator for
varying parameters. (c) Varying long l axis of elliptical cross section for
$\theta$ = 0, $\varphi$ = 0, $w=\qty{200}{\nano}$. (d) Varying polarization
$\varphi$ for l = 400, $w=\qty{200}{\nano}$, $\theta$ = 0. (e) Varying
incident angle $\theta$ for $\varphi$ = 0, l = 400, $w=\qty{200}{\nano}$.

The wavefront-shaping metasurface presented in this work is based on
multipolar-resonant silicon nanoresonators,[43, 44, 45, 46, 47, 42] with
elliptical cross sections of short axis w, long axis l and height h. For
initial characterization, the nanostructures are arranged in an equilateral
rhombohedral lattice with side length L = 670 (Fig. 1a) and periodic boundary
conditions. For our purposes, effective multipolar wavefront shaping is
achieved at infrared frequencies where silicon absorption becomes negligible
and electric and magnetic multipolar resonances overlap, thus enabling total
transmission blocking at resonant frequencies.[42] For a nanostructure of w =
200, l = 400 and h = 715, Fig. 1b shows cross sectional near-field plots of
the electric and magnetic field norms, demonstrating the simultaneous
occurrence of electric and magnetic multipolar resonances at 1004 and 1194
upon $y$-polarized incident light.

The aspect ratio of the elliptical cross sectional area dictates the resonant
behavior and optical tunability of the metasurface response. Figure 1c
demonstrates variation of the nanoresonator aspect ratio by changing the long
axis l of the structure from a circular to elliptical cross sectional area in
100 steps. Based on this parameter study, we select an ellipse with w = 200
and l = 400 for further investigation, as the multipolar resonances along the
long and short ellipse axes will be spaced sufficiently far apart to ensure
their independent behavior, while maintaining narrow-bandwidth resonant
properties.

The asymmetrical elliptical nanoresonator cross section and the rhombohedral
unit cell of the periodic boundaries result in an optically tunable
metasurface response. Figure 1d,e demonstrates how varying the polarization
($\varphi$) and incoming angle ($\theta$) can controllably tune the
transmittance at a given wavelength of interest. We define the polarization
angle $\varphi$ from the $\hat{\mathbf{y}}$ axis, i.e. the electric field
points along
$\hat{\mathbf{\epsilon}}=-\sin\varphi\hat{\mathbf{x}}+\cos\varphi\hat{\mathbf{y}}$
and the incidence angle is in the $xz$ plane so that the beam propagates along
$\hat{\mathbf{k}}=-\sin\theta\hat{\mathbf{x}}+\cos\theta\hat{\mathbf{z}}$.
Figure 1d ($0^{\circ}\leq\varphi\leq 160^{\circ}$ with $20^{\circ}$
increments, $\theta$ = 0∘), shows that varying the incident polarization state
$\varphi$ leads to complete tunability of the transmittance at the multipolar
resonant wavelengths $\lambda=\qty{1004}{\nano}$ and
$\lambda=\qty{1194}{\nano}$. In particular, transmittance is completely
blocked for excitation along the long axis of the ellipse
($\varphi=0^{\circ}$). Figure 1e ($\varphi$=0∘, $-10^{\circ}\leq\theta\leq
10^{\circ}$ with $2^{\circ}$ increments), demonstrates that varying the
incident angle $\theta$ does not significantly alter the metasurface
multipolar resonant transmittance. The optical properties observed in Fig.
1d,e indicate all-optical tunability of the studied multipolar resonant
metasurface and will now be leveraged to rationally design a wavefront-shaping
metasurface for SIM.

Figure 2:  Design and characterization of wavefront shaping metasurface for
SIM. (a) Metasurface geometry is composed of an equilateral rhombohedral
superlattice of 9 cells with side length 3$L$ = 2010 (left). Each cell
contains a nanoresonator as described in Fig. 1 at rotations of 0, 120 and 240
(right). (b) Electric field norm of the far-field patterns generated by the
metasurface at z = 5000 for varying incident angles of
$\theta=\\{-6^{\circ},0^{\circ},6^{\circ}\\}$ (bottom to top) and polarization
states $\varphi=\\{0^{\circ},120^{\circ},240^{\circ}\\}$.

To generate optically tunable far-field wavefronts which are suitable for SIM,
we designed a metasurface where the elliptical nanostructures described in
Fig. 1 are arranged in a periodically repeating geometry which generates
intensity patterns with high-contrast and large spatial frequency support.
These patterns are formed by constructive and destructive interference of
light from the multipolar nanoresonators characterized in Fig. 1. Using
multipolar resonances ensures the patterns have regions of low intensity, and
hence high contrast. The lattice geometry with rotational symmetry generates
the large spatial frequencies and ensures coverage in different spatial
directions, which is vital for achieving isotropic resolution enhancement. SIM
requires a diversity of different intensity patterns, and these can be
generated by varying the incident angle $\theta$ and polarization angle
$\varphi$ while retaining the high contrast and high spatial frequency content
imposed by the metasurface. Therefore, our metasurface generates all-optically
tunable, SIM-appropriate wavefront patterns.

Figure 2a shows a three-dimensional rendering (left) and top-down schematic
(right) of our SIM wavefront shaping metasurface. The metasurface consists of
a periodically arranged equilateral rhombohedral supercell with side length 3L
= 2010. The supercell contains 9 nanoresonators, with dimensions described in
Fig. 1, each located at the center of an equilateral rhombohedral sub-unit-
cell with side length L = 670. As shown in the right panel, the nanoresonators
are arranged at varying rotations of 0 (yellow), 120 (green) and 240 (blue),
respectively. Due to the rotational symmetry of the nanoresonator arrangement,
varying polarization states $\varphi$ will result in varying transmittance
values for each nanoresonator. For example for
$\varphi=\\{0^{\circ},120^{\circ},240^{\circ}\\}$, the nanoresonators marked
in yellow, green or blue, respectively, will be excited along their long axis,
resulting in transmittance blocking (see Fig. 1d). Moreover, varying incident
angles $\theta$ do not significantly alter the multipolar resonant wavelength
(Fig. 1e), yet varying $\theta$ will alter the propagation of transmitted
light to the far field for mutually rotated nanoresonators. This polarization
and angle-dependent excitation of the multipolar nanoresonators enables the
generation of versatile, SIM-appropriate wavefront patterns in the far field
through constructive and destructive superposition of the optical response of
each nanoresonator.

Figure 2b demonstrates examples of far field wavefront patterns generated by
the SIM wavefront shaping metasurface at height $z=\qty{5000}{\nano}$ from the
metasurface substrate for varying incident polarizations
$\varphi=\\{0^{\circ},120^{\circ},240^{\circ}\\}$ and angles
$\theta=\\{-6^{\circ},0^{\circ},6^{\circ}\\}$ (see Figures S2 and S3 in the
Supporting Information for further far-field analysis). Dashed lines in the
top left panel indicate the position of the nanoresonators at $z=0$. These
examples demonstrate the suitability of the generated wavefronts for SIM as
they obey the above-noted criteria of high contrast and large spatial
frequency support. Moreover, they exhibit similar wavefront patterns for a
given $\theta$ which rotate and translate at varying $\varphi$.

We can better understand the structure of the mrm-SIM patterns by considering
the underlying symmetry of the device. The metasurface consists of a periodic
array of pillars arranged on a hexagonal (triangular) lattice with primitive
lattice vectors $\mathbf{a}_{1}=\qty{670}{\nano}\times(3,0)$ and
$\mathbf{a}_{2}=\qty{670}{\nano}\times(3\cos\qty{60}{},\sin\qty{60}{})$,
together with a three element basis describing the position of the three non-
equivalent posts. The underlying lattice has $D_{6}$ dihedral symmetry while
the basis has $D_{3}$ dihedral symmetry. The rhombohedral supercell the
simulations are performed on contains three primitive cells (Fig. 2a). This
supercell has less symmetry than the full lattice, as its only symmetry is a
reflection plane 30 counterclockwise from the $x$-axis. This symmetry
guarantees that the patterns for $\varphi=\qty{0}{},\qty{240}{}$ polarization
are reflected copies of each other for $\theta=\qty{0}{}$ as illustrated in
Fig. 2b because these polarization directions are ±60 from the reflection
axis.

Figure 3:  SIM pattern optimization. (a) We select $7$ patterns from an
initial set of $21$ incidence angles $\theta$ and $3$ polarization angles
$\varphi$ by choosing a pattern set which optimizes the individual pattern
contrast (top right), flat sum intensity for the pattern set (center right),
and intensity balance between the various patterns (bottom right). Here the
green and blue curves are cartoons of 1D intensity line cuts of individual
patterns, and magenta of the sum pattern. (b) An optimal set of $7$ intensity
patterns and the incident angles and polarization angles $\theta,\varphi$.
Lattice vectors are shown in black (left, inset). (c) Simulated images of the
patterns from part b applied to a set of fluorescent line pairs with variable
spacing. The number of photons for each pixel is indicated by the color scale
(right panel).

Once the metasurface geometry is fixed, we select a set of illumination
patterns which are suitable for superresolution SIM from the range of
polarization and incidence angles simulated previously. We assume the electric
field above the metasurface is relayed to the sample, and defer discussion of
proposed experimental geometries. To choose an optimal set of patterns, we
developed an algorithm to select $N$ patterns which have similar
characteristics to the classical patterns used in sinusoidal and hexagonal
SIM. Our algorithm, which is illustrated in Fig. 3a, relies on the following
heuristic features of high-quality SIM patterns: (i) each individual pattern
should have high contrast, with regions where there is no optical intensity
(ii) the sum intensity of all the patterns should be spatially uniform (iii)
all patterns should have the same peak intensity (iv) all patterns should have
large spatial-frequency support. To quantify these features, we define a cost
function which formalizes the first three requirements as separate terms

$\mathcal{C}\left(\\{I_{i}\\}_{i=1}^{N}\right)=\frac{\text{sd}(\bar{I})}{\bar{I}}+\frac{1}{N}\sum_{i}\frac{\min
I_{i}}{\max I_{i}}+\frac{1}{N}\sum_{i}\left(1-\frac{\max
I_{i}}{I_{\max}}\right).$ (1)

Here $I_{i}$ is the $i$th intensity pattern and $\bar{I}$ and $I_{\max}$ are
the the average and maximum of the $N$ intensity patterns respectively. While
requirement (iv) could also be included in the cost function, we find all
patterns considered here have high-spatial frequency content so this is not
necessary. For $M$ simulated patterns there are $\binom{M}{N}$ pattern sets,
which is too large a space to exhaustively search. We randomly sample
$150\,000$ combinations and choose the pattern set with the lowest cost. Then,
we take this optimal set and iteratively replace each pattern by any patterns
not contained in the set. If this lowers the cost, we retain the new set of
patterns.

The structure of our mrm-SIM patterns more closely resemble the hexagonal SIM
patterns than they do sinusoidal SIM patterns. To perform a fair comparison
with standard hexagonal SIM we chose $N=$7$$ and set the maximum spatial
frequency of the hexagonal SIM patterns to match that in the mrm-SIM patterns
(section S7). The selected patterns, shown in Fig. 3b, fulfill requirements
(i)–(iv). The total cost from eq. 1 is $0.36$ with (i)–(iii) contributing
$0.14$, $0.09$, and $0.13$ respectively. We display the images which result
from applying these patterns to a synthetic fluorescent object in Fig. 3c.

Figure 4:  SIM pattern spatial Fourier structure. (a) Amplitude of the Fourier
transform of the phasor electric fields corresponding to the intensity
patterns shown in Fig. 3b. The maximum free-space propagating frequency
(green), microscope electric field band pass frequency (red), and vectors
$k_{\perp}(n_{1},n_{2})$ (purple) are shown. The electric field magnitude is
illustrated using the color scale, $\hat{\mathbf{\epsilon}}(\mathbf{k})_{1,2}$
(cyan, magenta) and $\hat{\mathbf{k}}$ (yellow). (b) Amplitude of the Fourier
transform of the intensity patterns shown in Fig. 3b (gamma = $0.3$). The
(intensity) band pass frequency of the microscope is shown in red. The
reciprocal lattice basis vectors are shown in blue, and other reciprocal
lattice frequencies within the band pass frequency are shown in purple. The
mrm-SIM patterns are primarily composed of $9$ Fourier modes.

We can better understand the spatial structure of the mrm-SIM patterns by
exploiting the periodic structure of the metasurface. Due to Bloch’s theorem,
the spatial frequency content of the mrm-SIM electric field patterns is
constrained by the primitive reciprocal lattice vectors, defined such that
$\mathbf{a}_{i}\cdot\mathbf{b}_{j}=\delta_{ij}$. In this case,
$\mathbf{b}_{1}=(1/3,-1/\sqrt{3})\times 1/\qty{670}{\nano}$ and
$\mathbf{b}_{2}=(0,2/\sqrt{3})\times 1/\qty{670}{\nano}$. A plane wave
incident through the glass with wavevector
$\mathbf{k}=k_{\perp}^{\text{in}}+k_{z}\hat{\mathbf{z}}$, where
$|\mathbf{k}|=2\pi n/\lambda$, interacts with the metasurface and couples into
various modes which have 2D wavevectors

$k_{\perp}(n_{1},n_{2})=k_{\perp}^{\text{in}}+2\pi\left(n_{1}\mathbf{b}_{1}+n_{2}\mathbf{b}_{2}\right)$
(2)

for $n_{1},n_{2}\in\mathbb{Z}$. On the other side of the metasurface, the
light is diffracted into the set of plane waves with wavevectors
$\mathbf{k}=k_{\perp}(n_{1},n_{2})+\sqrt{\left(2\pi/\lambda\right)^{2}-|k_{\perp}|^{2}}\hat{\mathbf{z}}$
and polarizations
$\hat{\mathbf{\epsilon}}_{2}(\mathbf{k})=\hat{\mathbf{k}}\times\hat{\mathbf{z}}/\left|\hat{\mathbf{k}}\times\hat{\mathbf{z}}\right|$
and
$\hat{\mathbf{\epsilon}}_{1}(\mathbf{k})=\hat{\mathbf{\epsilon}}_{2}(\mathbf{k})\times\hat{\mathbf{k}}$.
Where $|k_{\perp}(n_{1},n_{2})|\leq k$, these modes are propagating plane
waves which describe the far-field intensity pattern produced by the
metasurface. The intensity of the light that couples into each plane wave is
controlled by the physics of the light–metasurface interaction and varies with
incident angle and polarization.

We first consider the electric field structure produced by the metasurface for
orthogonal excitation ($\theta=\qty{0}{}$). Due to our choice of reciprocal
lattice vectors, only the shortest $6$ non-zero reciprocal lattice vectors
correspond with propagating modes in vacuum. In particular,
$\mathbf{b_{2}}\lambda\sim$0.999$$, and there are $5$ other lattice vectors of
the same length due to the lattice symmetry. Similar to a Bessel beam [48],
the plane-wave components at the symmetric spatial frequencies experience the
same phase evolution as they propagate and hence they form an invariant
pattern. The only pattern variation with $z$-position comes from interference
between this fixed pattern and the DC term. Introducing a finite-sized beam
will broaden these peaks and lead to additional divergence and pattern
evolution.

When a tilt angle is introduced, this shifts the supported frequency by
$k_{\perp}^{\text{in}}=-(2\pi n/\lambda)\sin\theta\hat{\mathbf{x}}$. Due to
the aggressive choice of reciprocal lattice vectors even a 1 tilt converts one
of the two supported reciprocal lattice vectors along the $x$-axis into an
evanescent wave. This effect is clear in the Fourier transform of the phasor
electric field shown illustrated in Fig. 4a. For small tilts, the spatial
frequencies maintain approximately the same magnitude and hence the patterns
evolve little over a significant $z$ propagation distance (see section S2). To
display these Fourier peaks clearly we apodize the real-space pattern with a
2D Hann window and Fourier broaden them by considering a ∼15 x 15 (500 x 500
pixel) area.

The structure of the mrm-SIM intensity patterns considered in Fig. 3b,c can be
understood from similar considerations. The Fourier transform of the intensity
pattern is related to the auto-correlation of the Fourier transform of the
electric field components, hence all supported spatial frequencies are the
difference of two electric field spatial frequencies, i.e.
$\mathbf{k}_{\text{int}}=\mathbf{k}_{\perp}(n_{1},n_{2})-\mathbf{k}_{\perp}(m_{1},m_{2})$
which are necessarily reciprocal lattice frequencies. For the intensity
patterns discussed above, we find only $9$ peaks out of the $19$ reciprocal
lattice vectors within the band pass contribute strongly (Fig. 4b), as these
are the only peaks formed by auto-correlation of the $4$ peaks in Fig. 4a. The
relative weight of the various Fourier peaks is shown in Fig. 4b. The most
important peaks for frequency resolution enhancement in SIM are the
$\pm\mathbf{b}_{2}$ peaks found at normalized frequency $\approx 0.66$ and
$\pm(2\mathbf{b}_{1}+\mathbf{b}_{2})$ at $\approx 0.38$. The absence of modes
at the highest spatial frequencies along directions other than
$\hat{\mathbf{y}}$ is due to the beam tilt effect discussed above.

Figure 5:  SIM reconstruction. (a) Widefield (WF), hexagonal SIM
reconstruction (H), and mrm-SIM reconstruction (M) for sets of line pairs with
spacings of 60450 in 30 steps. Arrows indicate the closest resolvable line
pairs which are 180 for hexagonal SIM (left arrow) and 240 for mrm-SIM. (b)
Widefield (purple), hexagonal SIM reconstruction (green), and mrm-SIM
reconstruction (orange) line cuts along the path in part a (cyan line). This
shows line pairs with spacings of 180330. (c) Widefield, hexagonal SIM, and
mrm-SIM reconstructions for synthetic filaments. (d) Line cuts through the
synthetic filament images along the path shown in part c (cyan lines).

To quantify the quality of our mrm-SIM approach, we simulated images of two
types of synthetic fluorescent samples: simulated line pairs with variable
spacing and simulated filamentous structures (see section S4 and S5). We then
reconstructed the underlying fluorescent object using an iterative
reconstruction approach described below. In all cases we assume a numerical
aperture of 1.3 and consider images with up to $1000$ photons per pixel, a
reasonable signal level for biological samples.

Due to the complexity of our patterns, applying standard SIM reconstruction
approaches cannot effectively infer superresolution information. Instead, we
adopt an iterative optimization-based reconstruction approach which relies on
our prior knowledge of the pattern structure. This approach infers the sample
structure using the fast iterative shrinkage-thresholding algorithm (FISTA)
[49], an accelerated proximal gradient-descent approach which incorporates
constraints and regularization (see section S6). In particular, we apply a
positivity constraint and use total-variation regularization to enforce
smoothness. A similar approach was recently applied to 3D SIM [50]. For
realistic instruments, blind-SIM style pattern recovery may be employed to
optimize the patterns [51].

We find mrm-SIM patterns produce high quality reconstructions of the synthetic
line pair images, as shown in Fig. 5a. For the line pairs we chose spacings of
60450 in steps of 30. mrm-SIM resolves the line pair separated by 240, while
the hexagonal SIM resolves the line pair separated by 180. Here we rely on the
Sparrow criterion. Both resolution values should be compared with the maximum
expected frequency enhancement in SIM, which is
$1+\frac{f_{\text{sim}}}{2\text{NA}/\lambda}$ based on the maximum spatial
frequency imaged by the microscope. For the patterns chosen here this is
$\approx 1.67$, which corresponds to a resolution of 230. We also compare the
widefield, hexagonal SIM reconstructed object and mrm-SIM reconstructed object
by taking a line cut (Fig. 5b). mrm-SIM clearly resolves line pairs which are
not resolved in the widefield image. However, it does not achieve as high
resolution or contrast as hexagonal SIM. We also considered a more realistic
model of a biological sample: simulated filament networks modelling
microtubules (Fig. 5c,d). This simulation shows similar results to the
synthetic line-pairs.

The modest difference in the resolution and contrast of mrm-SIM and hexagonal
SIM reconstructions can be understood by considering the distribution of
pattern weight at different spatial frequencies. The hexagonal SIM patterns
have the advantage that all of the Fourier weight is at either DC or
$k_{\text{sim}}$, and $I(k_{\text{sim}})/I(0)=1/3$. The mrm-SIM patterns have
significant weight at intermediate frequencies (Fig. 4b) and thus
$I(k_{\text{sim}})/I(0)\leq 0.1$. Therefore, for fixed intensity mrm-SIM
patterns contain less signal related to the highest frequency improvements.

Figure 6:  Proposed experimental realizations. (a) A proposed experimental
approach to realizing mrm-SIM. The excitation light (green) impinges on the
metasurface which is placed below the coverslip. Metasurface patterns excite
fluorescence (orange) in a sample. The fluorescence light is collected in
transmission by an imaging system and detected on a camera. (b) An
epifluorescence microscopy approach for mrm-SIM. The metasurface is placed in
a conjugate imaging plane which is optically relayed to the sample. A dichroic
mirror reflects the excitation light and transmits the fluorescence.

Although we have so far assumed the mrm-SIM patterns are relayed to the
sample, mrm-SIM is compatible with several different experimental geometries.
In the simplest scheme, the metasurface could be placed immediately underneath
the sample and coverslip (Fig. 6a) and be illuminated by a collimated beam
with angle control via a galvenometric mirror and polarization control via a
polarization rotator or electric-optic modulator. To avoid fluorescence light
interacting with the metasurface, fluorescence light should be collected in
transmission and excitation light rejected with a blocking filter. Embedding
the metasurface in the glass coverslip 5 from the surface would result in the
simplest alignment. Otherwise, the metasurface could be placed below the
coverslip, in which case the evolution of the mrm-SIM patterns through ∼100 of
glass must be considered. Alternatively, if the metasurface geometry is chosen
such that the excitation light is diffracted by the metasurface but, due to
the Stokes shift the fluorescence light is not, such a configuration might be
compatible with an epifluorescence arrangement, albeit for only a single
wavelength.

mrm-SIM can be applied to an epifluorescence geometry more generally by
placing the metasurface in a conjugate imaging plane and relaying the patterns
to the sample, as shown in Fig. 6b. The beam angle and polarization can be
controlled similar to the approach discussed above. The epifluorescence
approach makes mrm-SIM compatible with many commercial microscopes, but at the
cost of introducing a more complex optical train, precise alignment, and
limiting the achievable SIM pattern numerical aperture to $1$.

The dependence of the mrm-SIM patterns on the lattice parameters suggests the
natural way to tune their frequency is to change the lattice constant. For
example, decreasing the lattice constant by $\sim 2/3$ would presumably push
the SIM component at $\mathbf{b}_{2}$ to the band pass frequency. However, in
our current design the metasurface already operates at nearly the maximum
spatial frequency allowed for propagating the pattern in vacuum. Nevertheless,
modern high-numerical aperture microscope objectives support considerably
larger spatial frequencies, enabled by immersing the objective in water or
oil. Immersion increases the maximum propagating spatial frequency by a factor
of the refractive index of the immersion medium. Using an analogous approach,
immersing the metasurface in a liquid or solid medium with $n\sim$1.5$$ would
allow higher-frequency SIM patterns supporting larger resolution gains.

## Conclusion

In conclusion, we introduced the concept of mrm-SIM, where a wavefront-
shaping, multipolar-resonant metasurface acts as an ultrathin, single surface,
all-optical generator of SIM illumination patterns. The metasurface is
composed of silicon nanoresonators with elliptical cross-sectional area to
enable an incident angle and polarization tunable optical response. SIM-
appropriate illumination patterns were achieved by leveraging geometric
anisotropy of the nanoresonators, resulting in their mutually rotated assembly
in a periodic supercell. An optimization algorithm was developed and employed
to identify seven incident angles $\theta$ and polarization states $\varphi$
to generate SIM illumination patterns. Two types of fluorescent samples were
simulated to assess the quality of mrm-SIM, and the results were directly
compared to conventional hexagonal SIM. High quality image reconstructions
were achieved with our mrm-SIM approach, where resolutions of 240 were
observed, thus reaching similar resolutions to conventional hexagonal SIM
(180). Modest resolution differences were explained by the pattern weight
distribution at varying frequencies for each technique.

The proof of concept mrm-SIM demonstrated in this work operates at near-
infrared (NIR) wavelengths due to the optical properties of the silicon
nanoresonator material system. Such NIR illumination systems are suitable for
SIM research with upconverting dyes [52, 53] which, e.g., can enable new
applications in deep-tissue imaging. For future work, analogous SIM-compatible
metasurfaces which operate at visible frequencies can be developed with
alternative material systems such as SiN_x and TiO2 which exhibit low losses
and can thus support narrow bandwidth multipolar resonances such as those
presented in this work [25, 54]. Furthermore, it would be advantageous to
develop geometries which provide sufficient pattern diversity utilizing only
polarization rotation or metasurfaces embedded in high-refractive-index media
which can support SIM pattern frequencies with $\text{NA}\geq 1$. The
metasurface SIM platform demonstrated here makes steps toward significant
experimental simplification and democratization of SIM. The proposed all-
optical tunability of the metasurface-generated far-field wavefront on a
single, ultrathin metasurface paves the way toward next-generation SIM
experiments with larger fields of view, decreased alignment sensitivity,
miniaturized experimental footprint and minimized experimental complexity.

## acknowledgement

The authors thank Prachi Thureja and Zaid Haddadin for fruitful discussions.
All authors acknowledge funding from grant number 2021- 236170 from the Chan
Zuckerberg Initiative DAF, an advised fund of Silicon Valley Community
Foundation. PTB and DPS acknowledge funding from Scialog, Research Corporation
for Science Advancement, and Frederick Gardner Cottrell Foundation 28041. LVP
acknowledges funding from the Air Force Office of Scientific Research Young
Investigator Research Program Grant (FA9550-23-1-0263).

## Data availability

The electric and magnetic field data associated with figures 2–6, the
synthetic image generation code, and the SIM reconstruction code are available
online [55].

## References

  * Saxena  _et al._ [2015] M. Saxena, G. Eluru, and S. S. Gorthi, Structured illumination microscopy, Advances in Optics and Photonics 7, 241 (2015).
  * Heintzmann and Cremer [1999] R. Heintzmann and C. G. Cremer, Laterally modulated excitation microscopy: improvement of resolution by using a diffraction grating, in  _Optical biopsies and microscopic techniques III_, edited by I. J. Bigio, H. Schneckenburger, J. Slavik, K. Svanberg, and P. M. Viallet (SPIE, 1999).
  * Gustafsson [2000] M. G. L. Gustafsson, Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy., Journal of Microscopy 198, 82 (2000).
  * Prakash  _et al._ [2021] K. Prakash, B. Diederich, S. Reichelt, R. Heintzmann, and L. Schermelleh, Super-resolution structured illumination microscopy: past, present and future, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 379, 10.1098/rsta.2020.0143 (2021).
  * Ball  _et al._ [2015] G. Ball, J. Demmerle, R. Kaufmann, I. Davis, I. M. Dobbie, and L. Schermelleh, SIMcheck: A toolbox for successful super-resolution structured illumination microscopy, Scientific Reports 5, 10.1038/srep15915 (2015).
  * Demmerle  _et al._ [2017] J. Demmerle, C. Innocent, A. J. North, G. Ball, M. Müller, E. Miron, A. Matsuda, I. M. Dobbie, Y. Markaki, and L. Schermelleh, Strategic and practical guidelines for successful structured illumination microscopy, Nature Protocols 12, 988 (2017).
  * Hinsdale  _et al._ [2021] T. A. Hinsdale, S. Stallinga, and B. Rieger, High-speed multicolor structured illumination microscopy using a hexagonal single mode fiber array, Biomedical Optics Express 12, 1181 (2021).
  * Calvarese  _et al._ [2023] M. Calvarese, P. Paiè, A. Candeo, G. Calisesi, F. Ceccarelli, G. Valentini, R. Osellame, H. Gong, M. Neil, F. Bragheri, and A. Bassi, Integrated optical device for structured illumination microscopy, Optics Express 30, 30246 (2023).
  * Sandmeyer  _et al._ [2021] A. Sandmeyer, M. Lachetta, H. Sandmeyer, W. Hübner, T. Huser, and M. Müller, Cost-effective live cell structured illumination microscopy with video-rate imaging, ACS Photonics 8, 1639 (2021).
  * Brown  _et al._ [2021] P. T. Brown, R. Kruithoff, G. J. Seedorf, and D. P. Shepherd, Multicolor structured illumination microscopy and quantitative control of polychromatic light with a digital micromirror device, Biomedical Optics Express 12, 3700 (2021).
  * Li  _et al._ [2020] M. Li, Y. Li, W. Liu, A. Lal, S. Jiang, D. Jin, H. Yang, S. Wang, K. Zhanghao, and P. Xi, Structured illumination microscopy using digital micro-mirror device and coherent light source, Applied Physics Letters 116, 233702 (2020).
  * Helle  _et al._ [2020] Ø. I. Helle, F. T. Dullo, M. Lahrberg, J.-C. Tinguely, O. G. Hellesø, and B. S. Ahluwalia, Structured illumination microscopy using a photonic chip, Nature Photonics 14, 431 (2020).
  * Lin  _et al._ [2023] C. Lin, D. Schaubroeck, R. Baets, N. Boon, and N. L. Thomas, UV photonic-integrated-circuits-based structured illumination microscopy with a field of view larger than 100 $\mu$m2, IEEE Journal of Selected Topics in Quantum Electronics , 1 (2023).
  * Huang  _et al._ [2018] X. Huang, J. Fan, L. Li, H. Liu, R. Wu, Y. Wu, L. Wei, H. Mao, A. Lal, P. Xi, L. Tang, Y. Zhang, Y. Liu, S. Tan, and L. Chen, Fast, long-term, super-resolution imaging with Hessian structured illumination microscopy, Nature Biotechnology 36, 451 (2018).
  * Qiao  _et al._ [2022] C. Qiao, D. Li, Y. Liu, S. Zhang, K. Liu, C. Liu, Y. Guo, T. Jiang, C. Fang, N. Li, Y. Zeng, K. He, X. Zhu, J. Lippincott-Schwartz, Q. Dai, and D. Li, Rationalized deep learning super-resolution microscopy for sustained live imaging of rapid subcellular processes, Nature Biotechnology 10.1038/s41587-022-01471-3 (2022).
  * Maier  _et al._ [2007] S. A. Maier _et al._ ,  _Plasmonics: fundamentals and applications_ , Vol. 1 (Springer, 2007).
  * Wei and Liu [2010] F. Wei and Z. Liu, Plasmonic structured illumination microscopy, Nano Letters 10, 2531 (2010).
  * Bezryadina  _et al._ [2017] A. Bezryadina, J. Li, J. Zhao, A. Kothambawala, J. Ponsetto, E. Huang, J. Wang, and Z. Liu, Localized plasmonic structured illumination microscopy with an optically trapped microlens, Nanoscale 9, 14907 (2017).
  * Ponsetto  _et al._ [2017] J. L. Ponsetto, A. Bezryadina, F. Wei, K. Onishi, H. Shen, E. Huang, L. Ferrari, Q. Ma, Y. Zou, and Z. Liu, Experimental Demonstration of Localized Plasmonic Structured Illumination Microscopy, ACS Nano 11, 5344 (2017).
  * Bezryadina  _et al._ [2018] A. Bezryadina, J. Zhao, Y. Xia, X. Zhang, and Z. Liu, High Spatiotemporal Resolution Imaging with Localized Plasmonic Structured Illumination Microscopy, ACS Nano 12, 8248 (2018).
  * Yu  _et al._ [2011] N. Yu, P. Genevet, M. A. Kats, F. Aieta, J.-P. Tetienne, F. Capasso, and Z. Gaburro, Light propagation with phase discontinuities: generalized laws of reflection and refraction, Science 334, 333 (2011).
  * Arbabi  _et al._ [2015] A. Arbabi, Y. Horie, M. Bagheri, and A. Faraon, Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission, Nature nanotechnology 10, 937 (2015).
  * Lin  _et al._ [2014] D. Lin, P. Fan, E. Hasman, and M. L. Brongersma, Dielectric gradient metasurface optical elements, Science 345, 298 (2014).
  * Lindfors  _et al._ [2016] K. Lindfors, D. Dregely, M. Lippitz, N. Engheta, M. Totzeck, and H. Giessen, Imaging and steering unidirectional emission from nanoantenna array metasurfaces, ACS Photonics 3, 286 (2016).
  * Khorasaninejad  _et al._ [2016] M. Khorasaninejad, W. T. Chen, R. C. Devlin, J. Oh, A. Y. Zhu, and F. Capasso, Metalenses at visible wavelengths: Diffraction-limited focusing and subwavelength resolution imaging, Science 352, 1190 (2016).
  * Klopfer  _et al._ [2020] E. Klopfer, M. Lawrence, D. R. Barton III, J. Dixon, and J. A. Dionne, Dynamic focusing with high-quality-factor metalenses, Nano Letters 20, 5127 (2020).
  * Holsteen  _et al._ [2019] A. L. Holsteen, A. F. Cihan, and M. L. Brongersma, Temporal color mixing and dynamic beam shaping with silicon metasurfaces, Science 365, 257 (2019).
  * Thureja  _et al._ [2020] P. Thureja, G. K. Shirmanesh, K. T. Fountaine, R. Sokhoyan, M. Grajower, and H. A. Atwater, Array-level inverse design of beam steering active metasurfaces, ACS Nano 14, 15042 (2020).
  * Lawrence  _et al._ [2020] M. Lawrence, D. R. Barton III, J. Dixon, J.-H. Song, J. van de Groep, M. L. Brongersma, and J. A. Dionne, High quality factor phase gradient metasurfaces, Nature Nanotechnology 15, 956 (2020).
  * Huang  _et al._ [2013] L. Huang, X. Chen, H. Mühlenbernd, H. Zhang, S. Chen, B. Bai, Q. Tan, G. Jin, K.-W. Cheah, C.-W. Qiu, _et al._ , Three-dimensional optical holography using a plasmonic metasurface, Nature communications 4, 2808 (2013).
  * Ni  _et al._ [2013] X. Ni, A. V. Kildishev, and V. M. Shalaev, Metasurface holograms for visible light, Nature communications 4, 2807 (2013).
  * Zheng  _et al._ [2015] G. Zheng, H. Mühlenbernd, M. Kenney, G. Li, T. Zentgraf, and S. Zhang, Metasurface holograms reaching 80% efficiency, Nature nanotechnology 10, 308 (2015).
  * Silva  _et al._ [2014] A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, Performing mathematical operations with metamaterials, Science 343, 160 (2014).
  * Mohammadi Estakhri  _et al._ [2019] N. Mohammadi Estakhri, B. Edwards, and N. Engheta, Inverse-designed metastructures that solve equations, Science 363, 1333 (2019).
  * Cordaro  _et al._ [2023] A. Cordaro, B. Edwards, V. Nikkhah, A. Alù, N. Engheta, and A. Polman, Solving integral equations in free space with inverse-designed ultrathin optical metagratings, Nature Nanotechnology , 1 (2023).
  * Yu and Capasso [2014] N. Yu and F. Capasso, Flat optics with designer metasurfaces, Nature materials 13, 139 (2014).
  * Genevet  _et al._ [2017] P. Genevet, F. Capasso, F. Aieta, M. Khorasaninejad, and R. Devlin, Recent advances in planar optics: from plasmonic to dielectric metasurfaces, Optica 4, 139 (2017).
  * Kamali  _et al._ [2018] S. M. Kamali, E. Arbabi, A. Arbabi, and A. Faraon, A review of dielectric optical metasurfaces for wavefront control, Nanophotonics 7, 1041 (2018).
  * Shaltout  _et al._ [2019] A. M. Shaltout, V. M. Shalaev, and M. L. Brongersma, Spatiotemporal light control with active metasurfaces, Science 364, eaat3100 (2019).
  * Barton III  _et al._ [2020] D. Barton III, J. Hu, J. Dixon, E. Klopfer, S. Dagli, M. Lawrence, and J. Dionne, High-q nanophotonics: sculpting wavefronts with slow light, Nanophotonics 10, 83 (2020).
  * Thureja  _et al._ [2022] P. Thureja, R. Sokhoyan, C. U. Hail, J. Sisler, M. Foley, M. Y. Grajower, and H. A. Atwater, Toward a universal metasurface for optical imaging, communication, and computation, Nanophotonics 11, 3745 (2022).
  * Hail  _et al._ [2022] C. U. Hail, M. Foley, R. Sokhoyan, L. Michaeli, and H. A. Atwater, High quality factor metasurfaces for two-dimensional wavefront manipulation, arXiv preprint arXiv:2212.05647 (2022).
  * Schuller  _et al._ [2007] J. A. Schuller, R. Zia, T. Taubner, and M. L. Brongersma, Dielectric metamaterials based on electric and magnetic resonances of silicon carbide particles, Physical review letters 99, 107401 (2007).
  * Van de Groep and Polman [2013] J. Van de Groep and A. Polman, Designing dielectric resonators on substrates: Combining magnetic and electric resonances, Optics express 21, 26285 (2013).
  * Evlyukhin  _et al._ [2013] A. B. Evlyukhin, C. Reinhardt, E. Evlyukhin, and B. N. Chichkov, Multipole analysis of light scattering by arbitrary-shaped nanoparticles on a plane surface, JOSA B 30, 2589 (2013).
  * Miroshnichenko  _et al._ [2015] A. E. Miroshnichenko, A. B. Evlyukhin, Y. S. Kivshar, and B. N. Chichkov, Substrate-induced resonant magnetoelectric effects for dielectric nanoparticles, ACS Photonics 2, 1423 (2015).
  * Butakov and Schuller [2016] N. A. Butakov and J. A. Schuller, Designing multipolar resonances in dielectric metamaterials, Scientific reports 6, 38487 (2016).
  * Céspedes Vicente and Caloz [2021] O. Céspedes Vicente and C. Caloz, Bessel beams: A unified and extended perspective, Optica 8, 451 (2021).
  * Beck and Teboulle [2009] A. Beck and M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM Journal on Imaging Sciences 2, 183 (2009).
  * Cai  _et al._ [2022] M. Cai, H. Zhu, Y. Sun, L. Yin, F. Xu, H. Wu, X. Hao, R. Zhou, C. Kuang, and X. Liu, Total variation and spatial iteration-based 3d structured illumination microscopy, Optics Express 30, 7938 (2022).
  * Jost  _et al._ [2017] A. Jost, E. Tolstik, P. Feldmann, K. Wicker, A. Sentenac, and R. Heintzmann, Optical Sectioning and High Resolution in Single-Slice Structured Illumination Microscopy by Thick Slice Blind-SIM Reconstruction, PLOS ONE 10, e0132174 (2017).
  * Mettenbrink  _et al._ [2022] E. M. Mettenbrink, W. Yang, and S. Wilhelm, Bioimaging with Upconversion Nanoparticles, Advanced Photonics Research 3, 2200098 (2022).
  * Malhotra  _et al._ [2023] K. Malhotra, D. Hrovat, B. Kumar, G. Qu, J. V. Houten, R. Ahmed, P. A. E. Piunno, P. T. Gunning, and U. J. Krull, Lanthanide-doped upconversion nanoparticles: Exploring a treasure trove of nir-mediated emerging applications, ACS Applied Materials & Interfaces 15, 2499 (2023).
  * Yang  _et al._ [2020] J.-H. Yang, V. E. Babicheva, M.-W. Yu, T.-C. Lu, T.-R. Lin, and K.-P. Chen, Structural colors enabled by lattice resonance on silicon nitride metasurfaces, ACS nano 14, 5678 (2020).
  * Roy  _et al._ [2023] T. Roy, P. T. Brown, D. P. Shepherd, and L. V. Poulikakos, Data and code associated with "spatial wavefront shaping with a nanostructured metasurface for structured illumination microscopy" (2023).
  * Pierce and Spicer [1972] D. T. Pierce and W. E. Spicer, Electronic structure of amorphous Si from photoemission and optical studies, Physical Review B 5, 3017 (1972).
  * Brown  _et al._ [2023] P. T. Brown, S. Sheppard, A. Coullomb, and D. P. Shepherd, localize-psf (2023), 10.5281/zenodo.7662600.
  * Brown and Shepherd [2023] P. T. Brown and D. P. Shepherd, mcsim (2023), 10.5281/zenodo.7872728.
  * [59] T. Lambert, psfmodels, https://github.com/tlambert03/psfmodels.
  * Venkataramani  _et al._ [2016] V. Venkataramani, F. Herrmannsdörfer, M. Heilemann, and T. Kuner, SuReSim: Simulating localization microscopy experiments from ground truth models, Nature Methods 13, 319 (2016).
  * Chambolle [2004] A. Chambolle, An algorithm for total variation minimization and applications, Journal of Mathematical Imaging and Vision 20, 89 (2004).

## Supporting Information

##  S1 Finite Element Simulations of Metasurface

Full-field numerical computations in the frequency domain were performed with
COMSOL Multiphysics 6.0, resulting in time-harmonic electromagnetic fields
governed by Maxwell’s Equations. Phasor quantitites are assumed to carry time
dependence $\exp(i\omega t)$, implying phasor
$\exp(-i\mathbf{k}\cdot\mathbf{r})$ represents a plane wave propagating along
the $\hat{\mathbf{k}}$ direction. The optical properties of silicon [56] and
glass (Corning HPFS 7980 Fused Silica) were adapted from the COMSOL material
library, with the refractive index of the glass and the silicon domains at the
resonant wavelength of 1004 as $1.4508$ and $3.6356$, respectively. The glass
and the air domain spanned 5 each, with silicon nanostructures on the glass-
air interface (Fig. S1a).

Figure S1:  Schematic of the simulation domain. (a) Side view of the domain
indicating the dimensions, materials, and ports. (b) Top view (z-normal) of
the domain indicating the periodic boundaries. (c) Representative view of the
mesh.

Port boundary conditions were used on the glass and the air domains for
excitation and measuring the transmitted and reflected waves. A plane wave of
intensity 1000/ ^2 was applied along the positive z-axis on the port on the
glass domain. The port on the air domain was set to measure the transmitted
wave. Periodic boundary conditions were applied on the side walls of the glass
and the air domains (Fig. S1b). The mesh size was $\lambda/6n$ in the glass
and the air domains, where $\lambda$ is the wavelength of the incident light
and $n$ is the refractive index of the respective medium. The maximum element
size was 40 in the silicon domain (Fig. S1c). Analogous simulation parameters
were employed for the mrm-SIM metasurface supercell (Fig. 2).

##  S2 Far Field Analysis of Illumination Patterns

Figures S2 and S3 present additional information analyzing the mrm-SIM
patterns in the far field.

Figure S2 shows example electric field norms of the illumination patterns
generated by the multipolar resonant metasurface presented in this work for
varying incident polarization states
$\varphi=\\{0^{\circ},120^{\circ},240^{\circ}\\}$, on the top, middle and
bottom rows, respectively. The incident angle $\theta=0^{\circ}$ is studied in
all cases considered. Far-field illumination patterns are shown at varying
far-field heights of
$z=\\{\qty{4000}{\nano},\qty{4500}{\nano},\qty{5000}{\nano}\\}$ on the left,
center and right, respectively. Fig. S2 demonstrates how the illumination
patterns remain consistent at these varying heights. Moreover, while the
illumination patterns at $\varphi=0^{\circ}$ and $\varphi=240^{\circ}$ are
rotations of one another, the differences in the illumination patterns for
$\varphi=120^{\circ}$ arise from the symmetry of the illumination source with
respect to the metasurface unit cell. Because COMSOL uses a single plane wave
illumination with the studied port condition, symmetry of the metasurface unit
cell with respect to the incident plane wave is analogous but rotated for
$\varphi=0^{\circ}$ and $\varphi=240^{\circ}$, while $\varphi=120^{\circ}$ has
a different symmetry with respect to the metasurface unit cell as the incident
plane wave traverses the long axis of all three ellipses with axes along the
long diagonal of the rhombohedral simulation domain. In an experiment, these
differences in symmetry would average out due to a finite-sized illumination
beam.

Figure S3 shows the electric field norm for the studied simulation domains at
$\theta=0^{\circ}$ for varying incident polarizations
$\varphi=\\{0^{\circ},120^{\circ},240^{\circ}\\}$ from left to right. A cross
section of the simulation domain along the long diagonal of the rhombohedral
unit cell shows the full z-range. This view of the electric field norm
illustrates the consistence of the multipolar metasurface illumination
patterns observed in the far field.

Figure S2:  Far-field characterization of wavefront-shaping metasurface.
Electric field norm of wavefront patterns at varying heights from
$z=\\{\qty{4000}{\nano},\qty{4500}{\nano},\qty{5000}{\nano}\\}$ (left to
right) for varying incident polarizations
$\varphi=\\{0^{\circ},120^{\circ},240^{\circ}\\}$ (top to bottom) and incident
angle $\theta=0^{\circ}$.  Figure S3:  Far-field characterization of
wavefront-shaping metasurface. Electric field norm of wavefront patterns along
the long diagonal of the metasurface supercell at varying incident
polarizations $\varphi=\\{0^{\circ},120^{\circ},240^{\circ}\\}$ (left to
right) for incident angle $\theta=0^{\circ}$.

##  S3 Obtaining illumination patterns

The electromagnetic fields obtained from the COMSOL simulations must be
processed to determine the illumination intensity patterns on a square grid
suitable for simulating eq. S2. The COMSOL simulations are performed for a
single unit cell of the metasurface, but we would like to simulate
fluorescence densities covering many unit cells. To achieve this, we first
duplicate the electric field patterns at coordinates shifted by
$n_{1}\mathbf{c}_{1}+n_{2}\mathbf{c}_{2}$ for enough $n_{1},n_{2}$ to cover
our desired field of view. Here $\mathbf{c}_{1}=\mathbf{a}_{1}$ and
$\mathbf{c_{2}}=3\mathbf{a}_{2}-\mathbf{a}_{1}$ are the vectors describing the
periodicity of the simulation cell. Each shifted copy of the electric field
must be multiplied by a phase factor
$\exp\left[-i\mathbf{k}_{\text{in}}\cdot\left(n_{1}\mathbf{c}_{1}+n_{2}\mathbf{c}_{2}\right)\right]$
to account for the Bloch (Floquet) boundary conditions used in the simulation,
$\mathbf{k}_{\text{in}}$ is the wavevector of the plane wave illuminating the
metasurface and $|\mathbf{k}_{\text{in}}|=2\pi n/\lambda$ where $n$ is the
refractive index of the glass. We the interpolate the electric field
components onto a square grid and take the intensity pattern to be the modulus
squared. This effectively assumes the fluorescence object consists of emitters
with randomly oriented dipoles.

##  S4 Simulating realistic fluorescence images

We simulate realistic images using custom Python software, specifically our
localize_psf [57] and mcSIM [58] packages. To simulate a realistic microscope
image, we begin with a ground truth fluorophore density $S(\mathbf{r})$ and
illumination pattern $P(\mathbf{r})$. Then we apply the discretized microscope
forward model given in eq. S2, using a vectorial model of the microscope
point-spread function generated using psfmodels [59]. For our simulations we
assumed a numerical aperture of $1.3$, an object grid size 1368x1368 pixels
and pixel size of 29.25, and a raw image grid pixel size of 58.5. This grid is
oversampled with respect to Nyquist sampling. Next, we determine each pixel’s
“noisy” value by drawing from a Poisson distribution with mean value given by
the pixel value after applying the forward model. We model the action of an
sCMOS camera by assuming a gain of 2/ adding an offset of 100 and Gaussian
noise with standard deviation 2.5. We do not include pixel-to-pixel variation
in the camera parameters although this is expected to be present in real sCMOS
images.

##  S5 Synthetic fluorescent objects

The filament structures are generated from a filament network file provided by
Super Resolution Simulation (SuReSim) [60]. The positions given in the file
are shifted and scaled to match the desired field of view. The filament
network file provides vertices and the filaments themselves are drawn as one-
pixel wide lines.

Synthetic line pair images are generated with spacings of 60450 in steps of
30. The line patterns are rotated by angle of 15 from the pixel grid. Each
line has a flattop profile over 30 which then decays exponentially to zero
with decay length half a pixel. To mitigate pixelation effects, the ground
truth value at each pixel is calculated by subsampling the pixel 3x3 times.

##  S6 Iterative SIM reconstruction algorithm

Figure S4:  SIM reconstruction algorithm convergence. (a) The cost function
evaluated for each angle. Averaging these functions produces the cost function
of eq. 1 Results are shown for the synthetic line pairs for mrm-SIM (right
panel) and hexagonal SIM (left panel) over $100$ iterations of the algorithm.
(b) Mean-squared error between the reconstructed fluorophore density and the
ground truth values for mrm-SIM (purple) and hexagonal SIM (yellow).

Suppose we have (2D) fluorophore density $S(\mathbf{r})$ and SIM patterns
$P^{(i)}$ with ($i=1,...N$). The fluorescence intensity in the imaging plane
is given by

$\displaystyle I^{(i)}(\mathbf{r})$
$\displaystyle=\left[P^{(i)}(\mathbf{r})S(\mathbf{r})\right]*h(\mathbf{r})$
(S1)

where $h$ is the point-spread function of the imaging system and $*$ is the
convolution operator.

The camera detects the integrated irradiance over each pixel, so we can
rewrite this in discrete form

$\displaystyle I^{(i)}$ $\displaystyle=Bh\left[P^{(i)}\odot S\right]$ (S2)

where we regard $I$, $P^{(i)}$, and $S$ as vectors and $B$ and $h$ as linear
operators. Here $\odot$ represents elementwise multiplication and $B$ is the
binning operator which sums each $2\times 2$ pixel area. Including binning
allows us to define $S$ on a finer grid than the camera images, which is
necessary to infer superresolution information as the camera pixel size is
chosen to be Nyquist sampled with respect to the diffraction limit.

We determine the fluorophore density by solving

$\displaystyle S^{*}$
$\displaystyle=\mathop{\text{argmin}}_{S\in\mathcal{P}}\left[\mathcal{C}(S)+\tau\text{TV}(S)\right]$
(S3) $\displaystyle\mathcal{C}(S)$
$\displaystyle=\frac{1}{2Nn_{x}n_{y}}\sum_{i}\left\|I^{(i)}-D^{(i)}\right\|^{2}_{2}$
(S4)

where $D^{(i)}$ is the measured intensity for the $i$th pattern, $\mathcal{P}$
is the set of positive valued vectors $S$, TV is the total variation operator,
and $\tau$ is a hyperparameter controlling the influence of the regularization
term compared with the data penalty term, and $n_{x}n_{y}$ is the total number
of pixels in the image $D^{(i)}$.

We solve eq. S3 using the FISTA approach. Specifically, if $S_{k}$ is the
estimated fluorophore density at the $k$th iteration, then we generate the
next estimate according to

$\displaystyle y_{k+1}$
$\displaystyle=\text{prox}\left[S_{k}-\gamma\nabla_{S}\mathcal{C}(S)\right]$
(S5) $\displaystyle S_{k+1}$
$\displaystyle=y_{k+1}+\frac{q_{k}-1}{q_{k+1}}\left(y_{k+1}-y_{k}\right)$ (S6)
$\displaystyle q_{k+1}$
$\displaystyle=\frac{1}{2}\left(1+\sqrt{1+4q_{k}^{2}}\right).$ (S7)

Here we first apply gradient descent to our differentiable cost function, then
incorporate the regularization terms using their proximal operators (prox),
and finally accelerate the convergence of the scheme using a momentum type
approach. We apply proximal operators to (1) project $S$ onto the space of
positive vectors by setting negative values to zero and (2) enforce the TV
regularization [61]. Here $\gamma$ is the step-size, and when $\mathcal{C}$ is
convex and Lipschitz continuous with Lipschitz constant $L$ we can guarantee
convergence by choosing $\gamma<2/L$. We estimate $L$ using a power iteration
scheme. We initialize the scheme by setting $S_{0}$ equal to the average of
$D_{i}$ and $q_{0}=1$.

The gradient of the forward model required in eq. S5 is given by

$\displaystyle\nabla\mathcal{C}_{i}$
$\displaystyle=\frac{1}{n_{x}n_{y}}M^{(i){\dagger}}\left(M^{(i)}-D^{(i)}\right)$
(S8)

where $M^{(i)}=BhP^{(i)}$ and $M^{(i){\dagger}}$ is the Hermitian adjoint of
$M^{(i)}$. We find that better convergence is achieved by replacing the
gradient update with a stochastic gradient update where for each iteration we
randomly select a subset $\sigma_{k}\subseteq\\{1,...,N\\}$ and replace
$\nabla\mathcal{C}\to\frac{1}{|\sigma_{k}|}\sum_{i\in\sigma_{k}}\nabla\mathcal{C}_{i}$.

We illustrate the convergence of this algorithm for the synthetic line-pair
data presented in the main text in Fig. S4a for $100$ iterations. As expected
for an appropriately chosen step size, the cost function is always decreasing.
In Fig. S4b we show the mean-squared error defined by

$\displaystyle\epsilon_{k}$
$\displaystyle=\frac{1}{n_{x}n_{y}}\left\|\left(S_{k}-S_{\text{true}}\right)*h\right\|^{2}_{2}.$
(S9)

We find the error decreases rapidly for both sets of SIM patterns. The error
asymptotically approaches a limiting value after $100$ iterations, and the
hexagonal SIM achieves moderately smaller error than mrm-SIM due its greater
high-spatial frequency content.

##  S7 Hexagonal SIM patterns

In the main text we compare mrm-SIM patterns with hexagonal SIM patterns. The
hexagonal SIM patterns are given by

$\displaystyle I_{j}(\mathbf{r})$
$\displaystyle=\left|e^{i\mathbf{k}_{0}\cdot\mathbf{r}}+e^{i\mathbf{k}_{1}\cdot\mathbf{r}}e^{i\eta_{j}}+e^{i\mathbf{k}_{2}\cdot\mathbf{r}}e^{-i\eta_{j}}\right|^{2}$
(S10) $\displaystyle\mathbf{k_{i}}$
$\displaystyle=\frac{k_{\text{sim}}}{\sqrt{3}}\left(\cos\phi_{i},\sin\phi_{i}\right)$
(S11)

with $\phi_{i}=\frac{2\pi}{3}i$ for $i=0,1,2$ and $\eta_{j}=\frac{2\pi}{7}j$
for $j=0,...6$. In the main text we use
$k_{\text{sim}}=\frac{2}{3}k_{\max}=\frac{2}{3}\frac{\text{NA}}{\lambda}$