1. Introduction

Dielectric metasurfaces [1–3], due to their low optical losses, chemical stability, nontoxicity, and high refractive index enabling Mie-type resonances, represent a convenient choice for various applications in nanophotonics, such as flat optics [4], structural colors [5–10], nonlinear photonics [11–13], light structuring [14], holography [15] and augmented reality [16,17]. Typical high refractive index materials used in dielectric metasurfaces include silicon [18–20] and germanium [21,22]. These materials are lossy in the blue part of the spectrum, which hinders their use when all colors are needed, such as in achromatic metalenses and color generation for display applications. In this case, materials such as silicon nitride [8,9], gallium nitride [6], lithium niobate [23], titanium dioxide [5,10] and diamond [24,25] are preferred.

Optical metasurfaces are often made to operate at the single-resonant regime, as in the case of dielectric metasurface-based pixels individually designed for red, green and blue (RGB) colors [9,26]. On the other hand, multiresonant metasurfaces (also known as multiwavelength or multiband metasurfaces) can be used to have multiband compact devices, and are desired for applications in, e.g., biosensing [27], solar cells [28], nonlinear optics [29], color generation [25,30], achromatic metalenses [31], optical beam steering [32], and quantum state engineering [33]. Multiresont behavior specifically refers to achieving separate resonances, as opposed to the case of broadband response that can be achieved by means of many adjacent resonances. Multiresonant metasurfaces have been proposed by combining in a unit cell two or more different particles [30,34], or by using a single nanostructures supporting several multipole resonances [35]. In the context of surface lattice resonances (SLRs), multiple resonances with high quality factor have been demonstrated in dielectric [36] and plasmonic [37] metasurfaces. For example, by introducing many periodicities in the metasurface, using a Fourier Transform approach, one can achieve several SLRs at the desired wavelengths [38]. The nanofabrication process for some multiresonant configurations may become complex and costly. Having a relatively simple optical metasurface that exhibits multiple resonances at desired spectral positions with the possibility of tuning or switching them on-off would be advantageous.

In general, lossless materials such as diamond and titanium dioxide enable the manifestation of high-order Mie-type resonances in the visible with amplitude comparable to low-order modes, thus opening new opportunities for the engineering of resonances. Such high-order multipoles, although well understood in single nanostructures, are often suppressed by design in metasurfaces in order to achieve monochromatic response [26], or avoided, e.g., to reduce the analytical complexity associated with the modelling of collective resonances based on the dipole lattice sum theory. Furthermore, nanostructures supporting only low-order Mie-type resonances have a lower aspect ratio, which is easier to obtain in fabrication.

In this paper, we propose and theoretically demonstrate an approach to realizing multiresonant metasurfaces in the visible range, which is based on lattice-induced multipole coupling effects. The proposed metasurface design consists of regularly arranged (i.e., periodic lattice) array of dielectric nanocylinders with fixed diameter and height. Multiresonant response of the metasurface is governed by an energy exchange (coupling) between particular group of Mie-type multipoles. The paper is organized as follows. In Sec. 2, we present a theoretical model of multipole coupling in metasurfaces considering multipoles up to the 16-th order. In Sec. 3, we investigate the multipolar response of a dielectric nanostructure supporting high-order multipoles, and use it as a building block for a metasurface. Via full-wave numerical simulations, we observe the emergence and evolution of multiple resonances as a function of the unit cell size, and we explain which multipoles are contributing to each resonance based on the theoretical model introduced in Sec. 2. Up to three resonances at red, green and blue colors can be created. Control strategies to enable tunability of these resonances are discussed in Sec. 4, where we investigate how several parameters individually impact the multiple resonances, including geometric parameters (e.g., finite size of the metasurface and shape of the unit cell), illumination conditions (e.g., polarization), and surrounding conditions (e.g., the refractive index of a superstrate).

2. Theoretical background

In order to understand how high-order multipoles can be exploited for the creation of multiresonant metasurfaces, we first look at the multipole decomposition of reflection and transmission coefficients for an infinite periodic metasurface, and then present the coupling mechanisms between such multipoles.

2.1 Multipole decomposition of reflection and transmission coefficients

In order to facilitate the multipole decomposition for a metasurface, we consider an infinite periodic array of nanoparticles surrounded by a homogeneous medium with refractive index $n_s$. Furthermore, we assume that all nanoparticles are identical and have inversion symmetry (e.g., sphere, cylinder, cube). The electric field reflection $r$ (reflectance $R=|r|^2$) and transmission $t$ (transmittance $T=|t|^2$) coefficients of such system can be analytically calculated based on the multipole decomposition approach [1,39]. If the system is illuminated by a normal incidence $x$-axis polarized plane wave propagating along the $z$-axis, then due to translational symmetry, all nanoparticles in the metasurface have the same multipole moments and the coefficients $r$ and $t$ can be approximated as [39]

In contrast to previous works [1,39], the expressions for $r$ and $t$ are extended up to 16-th order poles. Higher order multipoles are relevant for this study as they significantly contribute in the visible spectral range at shorter wavelengths. On the other hand, Eqs. (1) are relatively complex since all multipole contributions are present. To better understand the physics described by these equations, we introduce effective electric and magnetic dipole moments, and reformulate Eqs. (1) as

2.2 Multipole coupling mechanism

Our system is invariant under inversion symmetry. However, multipole moments of the system behave differently under inversion operation. For the sake of simplicity, let us consider electric and magnetic dipole moment vectors. In Fig. 1(a), an electric dipole moment vector $\mathbf {p}$ has positive direction in the original black coordinate system. Under the inversion operation, the black coordinate system transforms to the red one. In the new inverted coordinate system, the vector $\mathbf {p}$ has negative coordinates, i.e., $\mathbf {p}\xrightarrow {i} -\mathbf {p}$ ($\xrightarrow {i}$ denotes inversion operation). Such vectors are called "true" (polar) vectors. In this context the multipole moment components also change sign in the inverted coordinate system (e.g., ED, MQ, EO and M16), and are known as odd parity (or simply odd) multipoles [40]. Meanwhile, a magnetic dipole moment vector $\mathbf {m}$ can be defined through a vector product of two "true" vectors (radius vector and $\mathbf {p}$). Since both "true" vectors change their sign under inversion operation, in the new inverted coordinate system, the vector $\mathbf {m}$ has positive coordinates, i.e., $\mathbf {m} \xrightarrow {i} \mathbf {m}$ [Fig. 1(b)]. Such vectors are called pseudovectors (axial vectors). Multipole moments, which retain the sign of their components under inversion operation, (such as MD, EQ, MO and E16) are known as even parity (or simply even) multipoles.

 

Fig. 1. Dipoles and inversion operation. Sketch that shows behavior of (a) electric and (b) magnetic dipole moment vectors under inversion operation (from black to red coordinate).

Download Full Size | PDF

In Fig. 2, 2D radiation diagrams of odd and even parity multipoles are illustrated. To understand the difference between them, let us concentrate on their electric field (black arrows). As we can see in Fig. 2(a), for the odd multipoles, the electric field has the same phase in opposite radiation directions. In the case of even multipoles, the phase of the electric field is opposite in those directions. This means that if fields generated by multipoles with different parities interfere constructively in one direction, then they interfere destructively in the opposite direction [41].

 

Fig. 2. Odd and even (under inversion) multipoles. 2D radiation patterns of (a) odd and (b) even parity multipoles. Arrow, circle (to reader) and cross (from reader) indicate direction of the electric (black) and magnetic (red) fields. (c) Odd-odd, (d) even-even, (e) odd-even and (f) even-odd multipole coupling (energy exchange) diagrams in an array.

Download Full Size | PDF

We now discuss how odd and even multipoles influence each other in periodic 2D arrays with inverse symmetry. First, we consider odd-odd multipole interaction, shown in Fig. 2(c). In the panel, part of a lattice with period $P_x$ and $P_y$ is shown. We assume that, at the lattice nodes, any type of odd multipole can be located (blue asterisk). At the position of the central node, the local electric field generated by surrounding multipoles is $\mathbf {E}_{\rm loc} \neq 0$, since the electric field of odd multipoles has the same phase in opposite directions [see Fig. 2(a)]. This means that there is an energy exchange (coupling) between odd multipoles in the lattice. The coupling between odd multipoles happens through their $\textit {electric field}$, which is a polar (true) vector. The magnetic fields of odd multipoles come to the central point in opposite phase and $\mathbf {H}_{\rm loc} = 0$. Nevertheless, the total magnetic field at that position can be non-zero because of even multipoles.

In a similar way, we can explain even-even multipole interaction, that is illustrated in Fig. 2(d). In this case, coupling between even multipoles happens through their $\textit {magnetic field}$, which is an axial vector. In this case, at the position of the central node, we have $\mathbf {E}_{\rm loc} = 0$ and $\mathbf {H}_{\rm loc} \neq 0$.

Now, let us consider cross, odd-even multipole interaction. It is known that multipoles with different parity do not couple if the system has inversion symmetry [42]. To illustrate it, the corresponding interaction diagram is shown in Fig. 2(e). The local electric field generated by even multipoles (surrounding green asterisks) at the position of the odd multipole (blue asterisk) is $\mathbf {E}_{\rm loc} = 0$, but the local magnetic field $\mathbf {H}_{\rm loc} \neq 0$. For instance, let us take an electric dipole with moment $\mathbf {p}$, and assume that there is a coupling between electric (blue asterisk) and magnetic (green asterisk) dipoles. Hence, one can write $\mathbf {p}=\alpha _{\rm p}\mathbf {E}_{\rm loc} + \gamma \mathbf {H}_{\rm loc}$, where $\alpha _{\rm p}$ is the dipole polarizability, $\gamma$ is some proportionality coefficient (bianisotropic polarizability) depending only on the particle properties. Under the inversion operation, $\mathbf {p}\xrightarrow {i} -\mathbf {p}$, and the fields $\mathbf {E}_{\rm loc}\xrightarrow {i} -\mathbf {E}_{\rm loc}$, $\mathbf {H}_{\rm loc}\xrightarrow {i} \mathbf {H}_{\rm loc}$ because the electric field is a true vector whereas the magnetic field is a pseudovector. Since for particles with inverse symmetry, $\alpha _{\rm p}$ and $\gamma$ are the invariant values under the inversion operation, from the relation $\mathbf {p}\xrightarrow {i} -\mathbf {p}$ it follows that the magnetic field $\mathbf {H}_{\rm loc}$ cannot determine the electric dipole and $\gamma \equiv 0$. Similar analysis with the same result can be done for even-odd multipole interaction shown in Fig. 2(f).

Thus, coupling (energy exchange) between multipoles is possible if they have the same parity. Due to this reason, the effective electric dipole moment $p^{\rm eff}_x$ [see Eqs. (3)] contains only odd multipole contributions, while the effective magnetic dipole moment $m^{\rm eff}_y$ includes only even ones.

3. Multiresonant metasurfaces: results and discussion

In this section, we design and explain multiresonance behavior of metasurfaces in the visible spectral range that exploits coupling effects between high-order multipoles. Hence, it is important that building blocks of the metasurface, i.e., single, isolated particles have to support high-order multipole response in the same spectral range. As a demonstration example, in this work we consider a cylindrical nanostructure with fixed diameter $D=200$ nm and height $H=300$ nm. Furthermore, we assume that surrounding and particle materials are dispersionless and homogeneous, and have refractive indices of $n_{\rm s}=1.45$ (similar to that of fused silica [43] in the visible) and $n_{\rm p}=2.45$ (similar to that of diamond [44], titanium dioxide [45] in the visible), respectively. These parameters allow having the high-order multipoles in the interested spectral range while being the system realistic.

3.1 Single particle resonances

The numerically calculated scattering efficiency of the single nanocylinder for different polarization and direction of plane wave excitation are shown in Fig. 3. Formulas for the scattering cross-section are provided in the Appendix A. From Fig. 3(a), we see that for the considered excitation condition (see inset sketch), the total scattering efficiency is relatively broad and the most of the incoming waves are scattered in the forward direction (see 3D far-field radiation pattern). Multipole content analysis of the total scattering efficiency via multipole decomposition shows that up to the 16-th order electric and magnetic poles can be resonantly excited in the considered wavelength range. Furthermore, individual multipole peaks have also relatively broad spectra and overlap, which makes possible efficient coupling between low (dipole, quadrupole) and high order (octupole, 16-th pole) multipoles. Similar behavior can be observed for the "side electric" Fig. 3(b) and "side magnetic" Fig. 3(c) excitation as well. Although, there are several important differences. First, the spectral positions of ED and MD resonances are flipped in comparison to Fig. 3(a), i.e., now the ED (MD) resonance is at longer (shorter) wavelength. Furthermore, as we can see in Fig. 3(b), there is a notable difference between black dashed (total scattering) and solid (sum of up to 16-th poles) lines below $\lambda = 400$ nm, which means that 32-th pole and higher poles are also excited.

 

Fig. 3. Single particle response. Multipole decomposition of the total scattering efficiency (black dashed line) of a single dielectric cylinder with refractive index $n_{\rm p}=2.45$, diameter $D=200$ nm and height $H=300$ nm in glass ($n_{\rm s}=1.45$) surrounding for different polarizations and directions of a plane wave excitation in the visible spectral range. Multipole sum is presented as a black solid line. The right-hand side of the panels presents the far-field radiation pattern at the spectral positions indicated by triangular, square and circular orange points. The scattering efficiency is $Q=\sigma /\sigma _0$, i.e., scattering cross-section $\sigma$ normalized to the particle’s geometrical cross-section $S$ which is (a) $\sigma _0=\pi D^2/4$, (b, c) $\sigma _0=HD$.

Download Full Size | PDF

On the right-hand side of Fig. 3, it is shown the 3D far-field radiation pattern at the wavelength indicated by triangular, square and circular orange points. As expected, we see more nonuniform far-field scattering at shorter wavelengths due to the dominating high order multipole contributions. While at long wavelengths, we see a typical dipole type radiation behavior. Note that the orange circular point in Fig. 3(c) corresponds to the well-known Kerker effect [41], i.e., no backscattering.

The same calculations were performed for a cuboid particle (height $H=300$ nm and side length of the square base $L=200$ nm), and multipoles up to the 16-th order were observed (not shown). As long as the single nanostructure exhibits such high order multipoles, we will observe multiple resonances in array configuration as presented in this paper.

3.2 Optical resonances of metasurfaces

We now use the cylinder nanostructure supporting high-order multipoles as the building block for a metasurface. We consider an array with a square lattice of period $P$, and investigate its optical response via numerical simulations. A sketch of the simulated metasurface and illumination condition are shown in Fig. 4(a). The metasurface is excited by a $z$-propagating and $x$-polarized plane wave. The refractive indices of materials and geometrical parameters of the cylinder are the same as in Fig. 3. Figure 4(b) shows a 2D color plot of the reflectance $R$ of the metasurface as a function of the lattice period $P$. The white solid line indicates the spectral position corresponding to the main order Rayleigh anomaly (RA, $\lambda ^{\rm RA}=n_{\rm s}P$) where a diffraction order starts to be an evanescent wave propagating in the plane of the array. Hence, the RA position divides the $P\lambda$-plane into two, diffractionless ($\lambda >\lambda ^{\rm RA}$) and diffractive ($\lambda <\lambda ^{\rm RA}$) regions. Here, we are particularly interested on the resonances occurring in the diffractionless region. Depending on the period $P$, the metasurface can have a different number of resonances in the visible and diffractionless regime. We identify four cases, numbered from 0 to 3, based on the number of resonances that occur in the reflectance spectrum. Below, we look at these particular cases more in detail by analyzing their multipole content.

 

Fig. 4. Metasurface resonances for different lattice periods. (a) Sketch of the considered infinite array of cylinders (with diameter $D$, height $H$ and index $n_{\rm p}$) placed in a homogeneous medium with index $n_{\rm s}$ and illuminated with the normal incidence (in $z$ direction) linearly polarized plane wave. (b) The reflectance $R$ of the metasurface [panel (a)] as a function of period $P$. The metasurface parameters are $D=200$ nm, $H=300$ nm, $n_{\rm s}=1.45$ and $n_{\rm p}=2.45$. The white dashed line displays the main Rayleigh anomaly (RA) condition ($\lambda ^{\rm RA}=n_{\rm s}P$), i.e., the area below this line corresponds to the diffractionless regime.

Download Full Size | PDF

We start with the case denoted by dashed line 0 in Fig. 4(b) corresponding to $P=475$ nm. In Fig. 5(a), the solid line correspond to the numerically calculated total reflectance, while the dashed line present the reflectance obtained via multipole approximation, i.e., Eqs. (2). One can note that (in other cases below as well) the full numerical and multipole approximation agree very well in the dispersionless regime ($\lambda >\lambda ^{\rm RA}$, right-hand side of the vertical dashed line), while the latter notably differs in the diffractive regime ($\lambda <\lambda ^{\rm RA}$, left-hand side of the vertical dashed line). The reason is that the current form of Eqs. (2) take into account only normally reflected and transmitted waves, i.e., zero-order diffraction, while full numerical results in our case include all diffraction order contributions. Oblique incidence case within dipole approximation is considered in Ref. [46].

 

Fig. 5. Resonances in cases 0 and 2. (a) Reflectance $R$ (solid line - full numerical, dashed line - multipole approximation) and (b) amplitudes $A$ of individual multipole contributions to the reflection coefficient at the numbered case 0 in Fig. 4 at $P=475$ nm. Panels (c) and (d) visualize the same calculation results for the numbered case 2 in Fig. 4(b) at $P=400$ nm.

Download Full Size | PDF

As observed in Fig. 5(a), we do not see any distinct resonance in the reflectance. This behavior can be explained by looking at the individual lattice multipole contributions to the reflection coefficient. Figure 5(b) visualizes the amplitudes of dominating multipoles contributions over wavelength. The solid (dashed) lines correspond to the amplitude $A$ of odd (even) multipoles. We see that around $\lambda =720~$nm all multipoles are in resonance condition and at the resonance position $R=0$. This is the well-known lattice Kerker effect, due to the suppression of (mainly) ED and MD radiations in backward direction [47]. They have nearly identical amplitude and phase in longer wavelengths as well. Since multipoles with same parity couple, as we discussed in Sec. 2.2, strong ED (MD) lattice resonance "pull up" MQ and EO (EQ) resonances. They also suppress each other in backward direction at the resonance (crossing of the green dashed and black solid lines at $R=0$ position). Using the introduced effective electric and magnetic dipole moments in Eqs. (2), the condition for the Kerker effect can be written as $p_x^{\rm eff}=m_y^{\rm eff}/v_{\rm s}$. Under the conditions of the resonant Kerker effect, electrical and magnetic energy is accumulated in the particles of the metasurface, as can be seen from the insets in Fig. 5(a).

Decreasing the period of the metasurface leads to detuning the spectral position of ED and MD resonances with changes of their Q-factors. Such behavior is a result of multipole coupling in the array. At the beginning, ED and MD contribute to a single resonance in the reflectance, e.g., the case denoted by dashed line 1 in Fig. 4(b) corresponding to $P=450$ nm. Further decreasing the period [e.g., see dashed line 2 in Fig. 4(b)] leads to a split of the single resonance into two resonances as one can also see in Fig. 5(c). This is because, now, the spectral separation between ED and MD resonances is large enough and each of them form a separate resonance in the reflection. Importantly, not only ED and MD, but also other higher order multipoles are also separated according to their parity. This is because they, i.e., the multipoles with the same parity, are coupled and behaves as an effective dipole.

Another interesting point is that at certain wavelengths, the amplitude of a multipole contribution becomes zero. For example, in Fig. 5(d), the amplitudes of the ED and MD are zero at around $\lambda =600$ nm and $\lambda =620$ nm, respectively. At those spectral positions the multipoles are suppressed. It is known that within the dipole approximation, the suppression happens at the RA position (indicated by vertical black dashed line), and is observed as $R=0$. However, in our case, the suppression positions of ED and MD are shifted to the long wavelength side, and hence, their amplitudes are nonzero at the RA. This is due to the coupling of the dipoles with higher order multipoles in an array, as was shown analytically in Ref. [48]. In particular, ED (MD) can disappear due to the vanishing of the local electric (magnetic) field caused by destructive interference between the incident external field and the field generated by odd (even) multipoles of a higher order. Due to the reciprocity principle, such suppression can also occur in relation to multipoles of higher order but at other spectral points.

In Fig. 5(d), we observe that the third resonance starts emerging around $\lambda =580$ nm. At that position, amplitudes of resonant even parity multipoles (indicated by a black dashed arrow) are a bit bigger than the odd ones (indicated by a black solid arrow) and due to partial suppression, there is partial reflection. Further, as we can see in Fig. 6(a), this resonance is "fully realized" for $P=350~$nm at $\lambda =510~$nm by making the metasurface to have three resonances. The reason is that, now, even (dashed arrow) and odd (solid arrow) parity multipole resonances are spectrally separated and provide partially constructive and destructive contributions into the maximum reflection. Note that due to the interference suppression of the dipole contributions, the reflection resonances have the asymmetric line-shape (Fano-type profile) [49].

 

Fig. 6. Resonances in case 3. (a) Reflectance $R$ (solid line - full numerical, dashed line - multipole approximation) and (b) amplitudes $A$ of individual multipole contributions to the reflection coefficient at the numbered case 3 shown in Fig. 4 (at $P=350$ nm). (c) Electric and (d) magnetic near-field distributions at the wavelengths indicated by triangle, asterisk, square and circle orange markers in (a).

Download Full Size | PDF

In Fig. 6(c,d), the near-field distributions on the $xy$-plane at the resonance positions (indicated by orange markers) are shown. For the middle resonance (square marker), the electric field is mostly confined in the space between particles and its intensity is enhanced up to 100 times. For such resonance, the average enhancement of the electric field intensity over a volume $V$ external to the nanoparticle is $1/V\int |{\rm E}/{\rm E}_0|^2dV\approx 25$ ($V$ in the $z$-direction is symmetric with respect to the nanoparticle, and has a height of 600 nm to include the near-field in all directions). This can have potential applications in areas such as sensing and nonlinear optics.

Finally, we would like to demonstrate the impact of the high-order multipoles on the spectrum. Equations (2) for the reflection $r$ and transmission $t$ coefficients include contributions of up to the 16-th order poles. Without those contributions, the calculations will be not accurate enough. This is demonstrated in Fig. 7, which shows the reflectance $R$ and transmittance $T$ (dashed lines) obtained within the dipoles, quadrupoles and octupoles approximations, respectively, for $P=350$ nm. We observe that dipoles and quadrupoles are not enough to qualitatively catch all resonances, especially for transmittance. It is necessary to consider at least up to octupoles to qualitatively reproduce all spectral features. However, even octupoles are not enough to preserve the energy conservation law.

 

Fig. 7. Impact of high order multipoles. Comparison of the full numerical (solid lines) versus multipole approximation based (dashed lines) reflectance and transmittance of a metasurface with period $P=350$ nm, (a) up to dipole, (b) quadrupole and (c) octupole contributions.

Download Full Size | PDF

We highlight that the system considered here can also exhibit four and five resonances, e.g., for periods $P=300$ nm and $P=250$ nm, respectively. However, for these periods the distance between the nanoparticles becomes very small, which makes the fabrication of such metasurfaces problematic. Furthermore, for an accurate description of the resonances occurring at shorter wavelengths the next order multipoles may be required.

4. Tunability of multiresonances: impact of size, lattice and external conditions

In order to explore the tunability of our multiresonant metasurface, in this section, we investigate the influence of inhomogeneous environment, absorption, array size, and rectangular lattice with different periods in different directions on the metasurface resonances.

Since real metasurfaces have finite size, we first investigate the influence of the metasurface size on the multiresonant response. In Fig. 8, the reflectance of the finite size metasurface with period $P=350$ nm for different array sizes is shown. The numerical approach for the estimation of the reflectance for a finite size metasurface is presented in Appendix C. In order to realize all three resonances with average $50{\%}$ of reflectance, the array has to have approximately 50 disks in each direction (in total 2500 disks). This in general holds for relatively broad resonances (e.g., quality factor ${\rm QF} \lesssim 100$). However, for a narrow resonances with ${\rm QF} \gtrsim 1000$, one can estimate that the number of particles has to be $\sim 100$ times larger (approximately 500 disks in each direction). The reason is that radiation loss due to the finite size effects can be considered as an absorption. Accordingly, in such lossy systems, realization of high quality factor resonances, which has long lifetime, is generally problematic.

 

Fig. 8. Impact of the array size. Reflectance of the metasurface with $P=350$ nm for different array size $N_x \times N_y$, where $N_x$ and $N_y$ are the number of particles in $x$ and $y$ directions, respectively.

Download Full Size | PDF

In many cases array of nanoparticles are fabricated on a glass substrate. In that case, a superstrate can be used to tune or control their optical response [20]. In Fig. 9, the reflectance of an infinite metasurface on a substrate with index $n_{\rm sub}=1.45$ and with period $P=350$ nm for different superstrate refractive indices $n_{\rm sup}$ is shown. The resonances are quite sensitive to the superstrate-to-substrate index contrast. For example, a single broad resonance associated with the resonant excitation of even multipole moments is observed in air ($n_{\rm sup}=1$). Immersion of the structure in water ($n_{\rm sup}=1.33$) leads to the appearance of the second relatively narrow resonance associated with the odd multipoles. As we already know, for the homogeneous case $n_{\rm sup}=n_{\rm sub}=1.45$, the metasurface exhibits three resonances. For $n_{\rm sup}>n_{\rm sub}$, there are only two resonances in the diffractionless region (the third short-wavelength resonance occurs in the diffractive region). These resonances become closer to the RA of the superstrate at longer wavelength and eventually disappear. Here, it is important to note that when $n_{\rm sup} \neq n_{\rm sub}$, the inversion symmetry is broken and due to the reflection at the interface, coupling between odd and even multipoles becomes possible [42].

 

Fig. 9. Impact of the superstrate. Reflectance of the metasurface on a substrate ($n_{\rm sub}=1.45$) with $P=350$ nm for different superstrate materials with refractive index $n_{\rm sup}$. The horizontal dashed line corresponds to the metasurface in a homogeneous medium ($n_{\rm sup}=n_{\rm sub}$). The vertical and inclined dashed lines indicate the spectral position of the main Rayleigh anomaly of glass substrate and superstrate, respectively.

Download Full Size | PDF

In Sec. 3.2, we considered the lattice with square elementary cells: $P=P_{x}=P_{y}$. However, as we have seen in Fig. 4(b), by varying the period one can change resonances strongly. Metasurface with a rectangular lattice ($P_{x} \neq P_{y}$) allows to switch between two different states by changing the excitation polarization. For example, as we see in Fig. 5(b), the square lattice metasurfaces with period $P=350$ nm and $P=450$ nm exhibit three and single resonances, respectively. We can switch back and forth between these states in the case of rectangular lattice with $P_{x}=350$ nm and $P_{y}=450$ nm. This is demonstrated numerically in Fig. 10(a). It is possible to "turn off" the first two short wavelength resonances (around $\lambda = 600$ nm and $\lambda = 650$ nm) by changing the polarization from $x$ ($\varphi =0$) to $y$ ($\varphi =90$ deg.), while the third long wavelength resonance ($\lambda = 700$ nm) experiences a small red shift. Similar polarization control of the colors (resonances) can be done other $P_x$ and $P_y$ combinations within or outside of Fig. 5(b). The lattice period can be optimized to tune the resonances around RGB colors, which is important for metasurface-based pixel application. For example, Fig. 10(b) depicts the above simulation for $P_{y}=275$ nm. The green resonance (letter G) position basically does not depend on the polarization, but the red (letter R) and blue (letter B) do. Hence, for the $x$ polarization, we mix green and red colors, while for the $y$ polarization, blue and green. In the middle, i.e., at $\varphi =45$ deg., we have essentially all three RGB colors. There are also two relatively narrow resonances between blue and green, and green and red resonances. In practice, they have less impact on the final spectra because of manufacturing imperfections and spectral resolution.

 

Fig. 10. Impact of the lattice shape. Reflectance of the metasurface with rectangular lattice for different polarization angle $\varphi$ for (a) $P_x=350$ nm and $P_y=450$ nm and (b) $P_x=350$ nm and $P_y=275$ nm. In (b), the top color bar visualizes the wavelength range of the $x$-axis. The left color bar shows the color appearance of the reflectance (using CIE 1931 color space and D65 light source) as a function of the angle $\varphi$.

Download Full Size | PDF

5. Conclusion

We have proposed and numerically demonstrated an approach to realize multiresonance all-dielectric metasurfaces in the visible spectral range based on high-order multipole moments. Our theoretical and numerical analysis shows that the physical mechanism behind the emergence of multiple resonances is the coupling between lattice multipoles with the same parity. In this mechanism, high order multipoles play a crucial role. High-order multipoles were enabled by using dielectric nanostructures made of materials with negligible losses in the visible, such as diamond or titanium dioxide. By engineering the period of a metasurface with a square lattice, i.e., by modulating the coupling between the nanostructures, it is possible to obtain an increasing number of resonances, up to three resonances simultaneously nearly at red, green and blue (RGB) wavelengths. Furthermore, by changing the refractive index of a superstrate (which can be effectively obtained, e.g., by changing the level of water or other index matching liquids), one can dynamically control the position of these resonances. A strategy based on a rectangular unit cell and tuning of the resonances based on the polarization angle of the incident wave was also presented, which may lead to the possibility of mixing colors, i.e., for different polarization angles, the reflectance can include resonances from different color regions. The theoretical details on multipole decomposition provided in this paper leads to a deeper understanding of the metasurface’s response, and may open new opportunities for the engineering of multiresonant metasurfaces, including the development of techniques for the indipendent tuning of each single resonance. All-dielectric multiresonant metasurfaces may find application in coloring technology, displays, multiband sensors, augmented reality, nonlinear nanophotonics (e.g., for sum frequency generation and four-wave mixing), and quantum state engineering.

Appendix A: multipole moments

This section contains formulas for multipole moments that have been used in Eqs. (2). We would like to list them explicitly for the sake of completeness and for convenience of a reader. Essentially, the exact detraced and symmetrized multipole moments up to E16 is provided in Ref. [50]. However, the expression for the latter one is more cumbersome than others. Therefore, we simply use the symmetrized primitive moment for M16 which is relatively easy to implement.

Electric dipole (ED) Vector of the electric dipole moment

(4)$${\bf p} = \frac{{\rm i}}{\omega} \int_{V}j_0(k_{\rm s}r){\bf j}dV+\frac{{\rm i}k_{\rm s}^2}{2\omega} \int_{V} \frac{j_2(k_{\rm s}r)}{(k_{\rm s}r)^2} [3({\bf j} \cdot {\bf r}){\bf r}-{r}^2{\bf j}]dV.$$

The scattering cross-section of the electric dipole radiation is

(5)$$\sigma_{\rm ED} = \frac{k_0^4}{6\pi\varepsilon_0^2|{\bf E}_0|^2} |{\bf p}|^2.$$

The $x$ component of ${\bf p}$ is

(6)$$p_x = \frac{{\rm i}}{\omega} \int_{V}j_0(k_{\rm s}r)j_xdV+\frac{{\rm i}k_{\rm s}^2}{2\omega} \int_{V} \frac{j_2(k_{\rm s}r)}{(k_{\rm s}r)^2} [3({\bf j} \cdot {\bf r})x-r^2j_x]dV.$$

Magnetic dipole (MD) Vector of the magnetic dipole moment

(7)$${\bf m} = \frac{3}{2} \int_{V}\frac{j_1(k_{\rm s}r)}{k_{\rm s}r}({\bf r} \times {\bf j})dV.$$

The scattering cross-section of the magnetic dipole radiation is

(8)$$\sigma_{\rm MD} = \frac{k_0^4\varepsilon_{\rm s}\mu_0}{6\pi\varepsilon_0|{\bf E}_0|^2} |{\bf m}|^2.$$

The $y$ component of ${\bf m}$ is

(9)$$m_y = \frac{3}{2} \int_{V}\frac{j_1(k_{\rm s}r)}{k_{\rm s}r}(zj_x-xj_z)dV.$$

Electric quadrupole (EQ) Tensor of the electric quadrupole moment

(10)$$\begin{aligned} {\hat{Q}} =& \frac{3{\rm i}}{\omega} \int_{V}\frac{j_1(k_{\rm s}r)}{k_{\rm s}r}[3({\bf r} {\bf j}+{\bf j} {\bf r})-2({\bf r} \cdot {\bf j})\hat{U}dV\\ +&\frac{6{\rm i}k_{\rm s}^2}{\omega} \int_{V}\frac{j_3(k_{\rm s}r)}{(k_{\rm s}r)^3}[5({\bf r} \cdot {\bf j}){\bf r}{\bf r}-r^2({\bf j} {\bf r}+{\bf r} {\bf j}) - r^2({\bf r} \cdot {\bf j})\hat{U}dV. \end{aligned}$$

The scattering cross-section of the electric quadrupole radiation is

(11)$$\sigma_{\rm EQ} = \frac{k_0^6\varepsilon_{\rm s}}{720\pi\varepsilon_0^2|{\bf E}_0|^2} \sum_{\alpha\beta}|Q_{\alpha\beta}|^2.$$

The $xz$ component of ${\hat {Q}}$ is

(12)$$\begin{aligned} Q_{xz} = &\frac{3{\rm i}}{\omega} \int_{V}\frac{j_1(k_{\rm s}r)}{k_{\rm s}r}[3(xj_z+j_xz)]dV\\ +&\frac{6{\rm i}k_{\rm s}^2}{\omega} \int_{V}\frac{j_3(k_{\rm s}r)}{(k_{\rm s}r)^3}[5({\bf r} \cdot {\bf j})xz-r^2(xj_z+j_xz)]dV. \end{aligned}$$

Magnetic quadrupole (MQ) Tensor of the magnetic quadrupole moment

(13)$${\hat{M}} = 5 \int_{V}\frac{j_2(k_{\rm s}r)}{(k_{\rm s}r)^2}[({\bf r} \times {\bf j}){\bf r}+{\bf r}({\bf r} \times {\bf j})]dV.$$

The scattering cross-section of the magnetic quadrupole radiation is

(14)$$\sigma_{\rm MQ} = \frac{k_0^6\varepsilon_{\rm s}^2\mu_0}{80\pi\varepsilon_0|{\bf E}_0|^2} \sum_{\alpha\beta}|M_{\alpha\beta}|^2.$$

The $yz$ component of ${\hat {M}}$ is

(15)$$M_{yz} = 5 \int_{V}\frac{j_2(k_{\rm s}r)}{(k_{\rm s}r)^2}[(zj_x-xj_z)z+y(xj_y-yj_x)]dV.$$

Electric octupole (EO) Tensor of the electric octupole moment

(16)$${\hat{O}^{({\rm e})}} = \frac{15{\rm i}}{\omega} \int_{V}\frac{j_2(k_{\rm s}r)}{(k_{\rm s}r)^2}({\bf j}{\bf r}{\bf r}+{\bf r}{\bf j}{\bf r}+{\bf r}{\bf r}{\bf j}-\hat{A}^{({\rm e})})dV,$$

where the components of the tensor $\hat {A}^{({\rm e})}$ are

(17)$$A_{\alpha\beta\gamma}^{({\rm e})} = \delta_{\alpha\beta}V_{\gamma}^{({\rm e})}+\delta_{\alpha\gamma}V_{\beta}^{({\rm e})}+\delta_{\beta\gamma}V_{\alpha}^{({\rm e})},\quad {\bf V}^{({\rm e})} = \frac{1}{5}[2({\bf r} \cdot {\bf j}){\bf r}+r^2{\bf j}].$$

The scattering cross-section of the electric octupole radiation is

(18)$$\sigma_{\rm EO} = \frac{k_0^8\varepsilon_{\rm s}^2}{1890\pi\varepsilon_0^2|{\bf E}_0|^2} \sum_{\alpha\beta\gamma}|{O}^{({\rm e})}_{\alpha\beta\gamma}|^2.$$

The $xzz$ component of ${\hat {O}^{({\rm e})}}$ is

(19)$$O^{({\rm e})}_{xzz} = \frac{15{\rm i}}{\omega} \int_{V}\frac{j_2(k_{\rm s}r)}{(k_{\rm s}r)^2}[(j_xz^2+2xzj_z)-\frac{1}{5}(2({\bf r} \cdot {\bf j})x+r^2j_x)]dV.$$

Magnetic octupole (MO) Tensor of the magnetic octupole moment

(20)$${\hat{O}^{({\rm m})}} = \frac{105}{4} \int_{V}\frac{j_3(k_{\rm s}r)}{(k_{\rm s}r)^3}[({\bf r} \times {\bf j}){\bf r}{\bf r}+{\bf r}({\bf r} \times {\bf j}){\bf r}+{\bf r}{\bf r}({\bf r} \times {\bf j})-\hat{A}^{({\rm m})})]dV,$$

where the components of the tensor $\hat {A}^{({\rm m})}$ are

(21)$$A_{\alpha\beta\gamma}^{({\rm m})} = \delta_{\alpha\beta}V_{\gamma}^{({\rm m})}+\delta_{\alpha\gamma}V_{\beta}^{({\rm m})}+\delta_{\beta\gamma}V_{\alpha}^{({\rm m})},\quad {\bf V}^{({\rm m})} = \frac{1}{5}r^2({\bf r} \times {\bf j}).$$

The scattering cross-section of the magnetic octupole radiation is

(22)$$\sigma_{\rm MO} = \frac{k_0^8\varepsilon_{\rm s}^3\mu_0}{1890\pi\varepsilon_0|{\bf E}_0|^2} \sum_{\alpha\beta\gamma}|{O}^{({\rm m})}_{\alpha\beta\gamma}|^2.$$

The $yzz$ component of ${\hat {O}^{({\rm m})}}$ is

(23)$$O^{({\rm m})}_{yzz} = \frac{105}{4} \int_{V}\frac{j_3(k_{\rm s}r)}{(k_{\rm s}r)^3}[(zj_x-xj_z)z^2+2(xj_y-yj_x)yz - \frac{1}{5}r^2(zj_x-xj_z)]dV.$$

Electric 16-pole (E16) Tensor of the electric 16-pole moment

(24)$${\hat{S}^{({\rm e})}} = \frac{105{\rm i}}{\omega} \int_{V}\frac{j_3(k_{\rm s}r)}{(k_{\rm s}r)^3}({\bf j}{\bf r}{\bf r}{\bf r}+{\bf r}{\bf j}{\bf r}{\bf r}+{\bf r}{\bf r}{\bf j}{\bf r}+{\bf r}{\bf r}{\bf r}{\bf j}-\hat{B}^{({\rm e})})dV,$$

where the components of the tensor $\hat {B}^{({\rm e})}$ are

(25)$$\begin{aligned} B_{\alpha\beta\gamma\tau}^{({\rm e})} =& \frac{1}{7}(\delta_{\alpha\beta}C_{\gamma\tau}^{({\rm e})}+\delta_{\alpha\gamma}C_{\beta\tau}^{({\rm e})}+\delta_{\alpha\tau}C_{\beta\gamma}^{({\rm e})}+\delta_{\beta\tau}C_{\alpha\gamma}^{({\rm e})}+\delta_{\beta\gamma}C_{\alpha\tau}^{({\rm e})}+\delta_{\gamma\tau}C_{\alpha\beta}^{({\rm e})})\\ +&\frac{4}{15}r^2({\bf r}\cdot{\bf j})(\delta_{\alpha\beta}\delta_{\gamma\tau}+\delta_{\alpha\gamma}\delta_{\beta\tau}+\delta_{\alpha\tau}\delta_{\beta\gamma}),\\ {\hat{C}^{({\rm e})}} =& 2({\bf j}\cdot{\bf r}){\bf r}{\bf r}+r^2({\bf j} {\bf r}+{\bf r}{\bf j})-\frac{4}{3}r^2({\bf r} \cdot {\bf j})\hat{U}. \end{aligned}$$

The scattering cross-section of the electric 16-pole radiation is

(26)$$\sigma_{\rm E16} = \frac{k_0^{10}\varepsilon_{\rm s}^3}{72576\pi\varepsilon_0^2|{\bf E}_0|^2} \sum_{\alpha\beta\gamma\tau}|{S}^{({\rm e})}_{\alpha\beta\gamma\tau}|^2.$$

The $xzzz$ component of ${\hat {S}^{({\rm e})}}$ is

(27)$$S^{({\rm e})}_{xzzz} = \frac{105{\rm i}}{\omega} \int_{V}\frac{j_3(k_{\rm s}r)}{(k_{\rm s}r)^3}[(3j_zxz^2+j_xz^3)-\frac{3}{7}(2({\bf r}\cdot{\bf j})xz+r^2(j_zx+j_xz))]dV.$$

Magnetic 16-pole (M16) Tensor of the magnetic 16-pole moment

(28)$${\hat{S}^{({\rm m})}} = 189 \int_{V}\frac{j_4(k_{\rm s}r)}{(k_{\rm s}r)^4}[({\bf r} \times {\bf j}){\bf r}{\bf r}{\bf r}+{\bf r}({\bf r} \times {\bf j}){\bf r}{\bf r}+{\bf r}{\bf r}({\bf r} \times {\bf j}){\bf r}+{\bf r}{\bf r}{\bf r}({\bf r} \times {\bf j})-\hat{B}^{({\rm m})}]dV,$$

where the components of the tensor $\hat {B}^{({\rm m})}$ are

(29)$$\begin{aligned} B_{\alpha\beta\gamma\tau}^{({\rm m})} &= \frac{1}{7}(\delta_{\alpha\beta}C_{\gamma\tau}^{({\rm m})}+\delta_{\alpha\gamma}C_{\beta\tau}^{({\rm m})}+\delta_{\alpha\tau}C_{\beta\gamma}^{({\rm m})}+\delta_{\beta\tau}C_{\alpha\gamma}^{({\rm m})}+\delta_{\beta\gamma}C_{\alpha\tau}^{({\rm m})}+\delta_{\gamma\tau}C_{\alpha\beta}^{({\rm m})}),\\ {\hat{C}^{({\rm m})}} &= r^2[({\bf r} \times {\bf j}){\bf r}+{\bf r}({\bf r} \times {\bf j})]. \end{aligned}$$

The scattering cross-section of the magnetic 16-pole radiation is

(30)$$\sigma_{\rm M16} = \frac{k_0^{10}\varepsilon_{\rm s}^4\mu_0}{72576\pi\varepsilon_0|{\bf E}_0|^2} \sum_{\alpha\beta\gamma\tau}|{S}^{({\rm m})}_{\alpha\beta\gamma\tau}|^2.$$

The $yzzz$ component of ${\hat {S}^{({\rm m})}}$ is

(31)$$S^{({\rm m})}_{yzzz} = 189 \int_{V}\frac{j_4(k_{\rm s}r)}{(k_{\rm s}r)^4}[3y(xj_y-yj_x)(z^2-\frac{1}{7}r^2)+z(zj_x-xj_z)(z^2-\frac{3}{7}r^2)]dV.$$

The sum of these individual scattering cross-sections approximate the numerical total cross-section $\sigma _{\rm Tot}$:

(32)$$\sigma_{\rm Tot}\approx\sigma_{\rm Sum} = \sigma_{\rm ED}+\sigma_{\rm MD}+\sigma_{\rm EQ}+\sigma_{\rm MQ}+\sigma_{\rm EO}+\sigma_{\rm MO}+\sigma_{\rm E16}+\sigma_{\rm M16}.$$

In the above equations: $\mu _0$ is the vacuum permeability, $\omega$ is the circular frequency, r is the vector of the coordinates of a point dipole located in the particle volume $V$, $r=|{\bf r}|$, ${\bf j}={\bf j}({\bf r})$ is the current density vector inside the particle volume $V$, ${\bf r}\cdot {\bf j}=xj_x+yj_y+zj_z$, $j_n(x)$ is the $n$-th order (in general $n \in \mathbb {C}$) spherical Bessel function (dot not mix it with components of the current density $j_x$, $j_y$ and $j_z$ which uses letter indices), $\delta _{\alpha \beta }$ is the Kronecker delta, $\hat {U}$ is the unit tensor. The Greek letter indices $\alpha,\beta$ stand for $x$, $y$ and $z$. Terms like ${\bf r}{\bf j}\equiv {\bf r}\otimes {\bf j}$ means outer product (also known as tensor, dyadic product) between vectors ${\bf r}$ and ${\bf j}$.

Appendix B: how to calculate the multipole moments numerically?

We calculated the multipole moments in ANSYS Lumerical’s FDTD or DGTD 3D electromagnetic solver. To calculate the multipole moments, we need to know the electric field ${\bf E}$ inside a particle. When ${\bf E}$ is known, the current density can be calculated as

(33)$${\bf j} = {-{\rm i}}\omega\varepsilon_0(\varepsilon_{\rm p}-\varepsilon_{\rm s}){\bf E}={-{\rm i}}\omega({\bf D}-\varepsilon_0\varepsilon_{\rm s}{\bf E}),$$

where ${\bf D}$ is the electric displacement field. We use "Displacement field" analysis object of ANSYS Lumerical (it is available within the FDTD solver, see "Advanced analysis" category in "Object Library") which calculates E, D and the 3D grid of the analysis region. The latter allows us to extract the fields only inside the particle.

One can note that Eqs. (2) and the scattering cross-sections contain the electric field ${\bf E}_0$ of the incident wave at the center position of the particle. To calculate ${\bf E}_0$ first, we run the simulation without the particle by disabling it (or by assigning to it the refractive index of the surrounding material) and save the ${\bf E}_0$ at the center position of the particle. This can be done by using a "point field monitor". After that, we can use the saved ${\bf E}_0$ in the second simulation with the particle. Furthermore, the "CW normalization (first source)" has to be selected in the settings of ANSYS Lumerical (Setting→Normalization state) in both simulations. Note that for Eqs. (2) the fields must be calculated with periodic boundary conditions along $x$ and $y$ directions, and perfectly matched layer (PML) boundary condition along $z$ direction, while for the scattering cross-section calculations, PML must be used in all $x$, $y$ and $z$ directions.

Appendix C: calculation of reflectance of a finite size metasurface

Reflectance of a finite size metasurface are simulated by using the “Transmission box” analysis group object of ANSYS Lumerical. We define the reflectance as a backscattered power normalized to the total-field scattered-field (TFSF) source power. Sketch of the setup is shown in Fig. 11. The backscattered power is calculated by placing the “Transmission box” outside of the TFSF source box and power going through the solid blue line is counted. To minimize the influence of the TFSF source size on simulation results, we keep half period $P/2$ gap between metasurface, TFSF box, and “Transmission box” as shown in Fig. 11. The size of the TFSF box in the propagation direction is 2 $\mathrm {\mu }$m.

 

Fig. 11. Finite size metasurface and TFSF source. Positions of the total-field scattered-field (TFSF) source box and “Transmission box” analysis object to calculate the reflectance of a finite size metasurface.

Download Full Size | PDF

Funding

Deutsche Forschungsgemeinschaft (EXC 2122, Project ID 390833453).

Acknowledgments

We acknowledge the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy within the Cluster of Excellence PhoenixD (EXC 2122, Project ID 390833453). Simulations were performed on the central computing cluster operated by Leibniz University IT Services (LUIS) , which is funded by the DFG (project number INST 187/742-1 FUGG). We also acknowledge the computing time granted by the Resource Allocation Board and provided on the supercomputer Lise and Emmy at NHR@ZIB and NHR@Göttingen as part of the NHR infrastructure (project nip00059).

Disclosures

The authors declare no conflicts of interest.

Data availability

All the data supporting this study are available from the authors upon reasonable request.

References

1. A. B. Evlyukhin, C. Reinhardt, A. Seidel, et al., “Optical response features of Si-nanoparticle arrays,” Phys. Rev. B 82(4), 045404 (2010). [CrossRef]  

2. P. Genevet, F. Capasso, F. Aieta, et al., “Recent advances in planar optics: From plasmonic to dielectric metasurfaces,” Optica 4(1), 139–152 (2017). [CrossRef]  

3. A. I. Kuznetsov, A. E. Miroshnichenko, M. L. Brongersma, et al., “Optically resonant dielectric nanostructures,” Science 354(6314), aag2472 (2016). [CrossRef]  

4. A. Arbabi, Y. Horie, M. Bagheri, et al., “Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission,” Nat. Nanotechnol. 10(11), 937–943 (2015). [CrossRef]  

5. S. Sun, Z. Zhou, C. Zhang, et al., “All-dielectric full-color printing with TiO₂ metasurfaces,” ACS Nano 11(5), 4445–4452 (2017). [CrossRef]  

6. B. H. Chen, P. C. Wu, V.-C. Su, et al., “GaN metalens for pixel-level full-color routing at visible light,” Nano Lett. 17(10), 6345–6352 (2017). [CrossRef]  

7. W. Yang, S. Xiao, Q. Song, et al., “All-dielectric metasurface for high-performance structural color,” Nat. Commun. 11(1), 1864 (2020). [CrossRef]  

8. C.-S. Park, I. Koirala, S. Gao, et al., “Structural color filters based on an all-dielectric metasurface exploiting silicon-rich silicon nitride nanodisks,” Opt. Express 27(2), 667–679 (2019). [CrossRef]  

9. J.-H. Yang, V. E. Babicheva, M.-W. Yu, et al., “Structural colors enabled by lattice resonance on silicon nitride metasurfaces,” ACS Nano 14(5), 5678–5685 (2020). [CrossRef]  

10. H. Li, Y. Xu, X. Zhang, et al., “All-dielectric high saturation structural colors enhanced by multipolar modulated metasurfaces,” Opt. Express 30(16), 28954–28965 (2022). [CrossRef]  

11. G. Li, S. Zhang, and T. Zentgraf, “Nonlinear photonic metasurfaces,” Nat. Rev. Mater. 2(5), 17010 (2017). [CrossRef]  

12. K. Koshelev, S. Kruk, E. Melik-Gaykazyan, et al., “Subwavelength dielectric resonators for nonlinear nanophotonics,” Science 367(6475), 288–292 (2020). [CrossRef]  

13. C. Gigli, G. Marino, A. Artioli, et al., “Tensorial phase control in nonlinear meta-optics,” Optica 8(2), 269–276 (2021). [CrossRef]  

14. R. C. Devlin, A. Ambrosio, N. A. Rubin, et al., “Arbitrary spin-to–orbital angular momentum conversion of light,” Science 358(6365), 896–901 (2017). [CrossRef]  

15. M. Khorasaninejad, A. Ambrosio, P. Kanhaiya, et al., “Broadband and chiral binary dielectric meta-holograms,” Sci. Adv. 2(5), e1501258 (2016). [CrossRef]  

16. S.-Q. Li, X. Xu, R. Maruthiyodan Veetil, et al., “Phase-only transmissive spatial light modulator based on tunable dielectric metasurface,” Science 364(6445), 1087–1090 (2019). [CrossRef]  

17. J. Xiong, E.-L. Hsiang, Z. He, et al., “Augmented reality and virtual reality displays: emerging technologies and future perspectives,” Light: Sci. Appl. 10(1), 216 (2021). [CrossRef]  

18. E. Højlund-Nielsen, J. Weirich, J. Nørregaard, et al., “Angle-independent structural colors of silicon,” J. Nanophotonics 8(1), 083988 (2014). [CrossRef]  

19. Z. Dong, J. Ho, Y. F. Yu, et al., “Printing beyond sRGB color gamut by mimicking silicon nanostructures in free-space,” Nano Lett. 17(12), 7620–7628 (2017). [CrossRef]  

20. I. Allayarov, A. B. Evlyukhin, D. J. Roth, et al., “Dynamic nonlocal dielectric metasurfaces: Tuning collective lattice resonances via substrate–superstrate permittivity contrast,” Adv. Photonics Res. 5(1), 2300268 (2024). [CrossRef]  

21. X. Zhu, W. Yan, U. Levy, et al., “Resonant laser printing of structural colors on high-index dielectric metasurfaces,” Sci. Adv. 3(5), e1602487 (2017). [CrossRef]  

22. Z.-X. Zhou, M.-J. Ye, M.-W. Yu, et al., “Germanium metasurfaces with lattice Kerker effect in near-infrared photodetectors,” ACS Nano 16(4), 5994–6001 (2022). [CrossRef]  

23. A. Fedotova, L. Carletti, A. Zilli, et al., “Lithium niobate meta-optics,” ACS Photonics 9(12), 3745–3763 (2022). [CrossRef]  

24. J. Hu, M. Lawrence, and J. A. Dionne, “High quality factor dielectric metasurfaces for ultraviolet circular dichroism spectroscopy,” ACS Photonics 7(1), 36–42 (2020). [CrossRef]  

25. J. Gu, Y. Liu, N. Meng, et al., “Structural colors based on diamond metasurface for information encryption,” Adv. Opt. Mater. 11(6), 2202826 (2023). [CrossRef]  

26. B. Yang, W. Liu, Z. Li, et al., “Ultrahighly saturated structural colors enhanced by multipolar-modulated metasurfaces,” Nano Lett. 19(7), 4221–4228 (2019). [CrossRef]  

27. D. Rodrigo, A. Tittl, N. Ait-Bouziad, et al., “Resolving molecule-specific information in dynamic lipid membrane processes with multi-resonant infrared metasurfaces,” Nat. Commun. 9(1), 2160 (2018). [CrossRef]  

28. E. F. Pecora, A. Cordaro, P. G. Kik, et al., “Broadband antireflection coatings employing multiresonant dielectric metasurfaces,” ACS Photonics 5(11), 4456–4462 (2018). [CrossRef]  

29. M. J. Huttunen, O. Reshef, T. Stolt, et al., “Efficient nonlinear metasurfaces by using multiresonant high-Q plasmonic arrays,” J. Opt. Soc. Am. B 36(7), E30–E35 (2019). [CrossRef]  

30. Y. Nagasaki, M. Suzuki, and J. Takahara, “All-dielectric dual-color pixel with subwavelength resolution,” Nano Lett. 17(12), 7500–7506 (2017). [CrossRef]  

31. F. Aieta, M. A. Kats, P. Genevet, et al., “Multiwavelength achromatic metasurfaces by dispersive phase compensation,” Science 347(6228), 1342–1345 (2015). [CrossRef]  

32. A. Calà Lesina, D. Goodwill, E. Bernier, et al., “Tunable plasmonic metasurfaces for optical phased arrays,” IEEE J. Sel. Top. Quantum Electron. 27(1), 1–16 (2021). [CrossRef]  

33. T. Santiago-Cruz, S. D. Gennaro, O. Mitrofanov, et al., “Resonant metasurfaces for generating complex quantum states,” Science 377(6609), 991–995 (2022). [CrossRef]  

34. S. Murai, D. R. Abujetas, L. Liu, et al., “Engineering bound states in the continuum at telecom wavelengths with non-bravais lattices,” Laser Photonics Rev. 16(11), 2100661 (2022). [CrossRef]  

35. A. B. Evlyukhin, C. Reinhardt, U. Zywietz, et al., “Collective resonances in metal nanoparticle arrays with dipole-quadrupole interactions,” Phys. Rev. B 85(24), 245411 (2012). [CrossRef]  

36. G. W. Castellanos, P. Bai, and J. Gómez Rivas, “Lattice resonances in dielectric metasurfaces,” J. Appl. Phys. 125(21), 213105 (2019). [CrossRef]  

37. O. Reshef, M. Saad-Bin-Alam, M. J. Huttunen, et al., “Multiresonant high-Q plasmonic metasurfaces,” Nano Lett. 19(9), 6429–6434 (2019). [CrossRef]  

38. T.-L. Lim, Y. Vaddi, M. S. Bin-Alam, et al., “Fourier-engineered plasmonic lattice resonances,” ACS Nano 16(4), 5696–5703 (2022). [CrossRef]  

39. P. D. Terekhov, V. E. Babicheva, K. V. Baryshnikova, et al., “Multipole analysis of dielectric metasurfaces composed of nonspherical nanoparticles and lattice invisibility effect,” Phys. Rev. B 99(4), 045424 (2019). [CrossRef]  

40. M. Poleva, K. Frizyuk, K. Baryshnikova, et al., “Multipolar theory of bianisotropic response of meta-atoms,” Phys. Rev. B 107(4), L041304 (2023). [CrossRef]  

41. W. Liu and Y. S. Kivshar, “Generalized Kerker effects in nanophotonics and meta-optics,” Opt. Express 26(10), 13085–13105 (2018). [CrossRef]  

42. D. A. Bobylev, D. A. Smirnova, and M. A. Gorlach, “Nonlocal response of Mie-resonant dielectric particles,” Phys. Rev. B 102(11), 115110 (2020). [CrossRef]  

43. I. H. Malitson, “Interspecimen comparison of the refractive index of fused silica,” J. Opt. Soc. Am. 55(10), 1205–1209 (1965). [CrossRef]  

44. A. Taylor, P. Ashcheulov, P. Hubík, et al., “Comparative determination of atomic boron and carrier concentration in highly boron doped nano-crystalline diamond,” Diamond Relat. Mater. 135, 109837 (2023). [CrossRef]  

45. A. Jolivet, C. Labbé, C. Frilay, et al., “Structural, optical, and electrical properties of TiO₂ thin films deposited by ald: Impact of the substrate, the deposited thickness and the deposition temperature,” Appl. Surf. Sci. 608, 155214 (2023). [CrossRef]  

46. D. R. Abujetas, J. Olmos-Trigo, J. J. Sáenz, et al., “Coupled electric and magnetic dipole formulation for planar arrays of particles: Resonances and bound states in the continuum for all-dielectric metasurfaces,” Phys. Rev. B 102(12), 125411 (2020). [CrossRef]  

47. V. E. Babicheva and A. B. Evlyukhin, “Resonant lattice Kerker effect in metasurfaces with electric and magnetic optical responses,” Laser Photonics Rev. 11(6), 1700132 (2017). [CrossRef]  

48. V. E. Babicheva and A. B. Evlyukhin, “Analytical model of resonant electromagnetic dipole-quadrupole coupling in nanoparticle arrays,” Phys. Rev. B 99(19), 195444 (2019). [CrossRef]  

49. A. E. Miroshnichenko, S. Flach, and Y. S. Kivshar, “Fano resonances in nanoscale structures,” Rev. Mod. Phys. 82(3), 2257–2298 (2010). [CrossRef]  

50. A. B. Evlyukhin and B. N. Chichkov, “Multipole decompositions for directional light scattering,” Phys. Rev. B 100(12), 125415 (2019). [CrossRef]