A simple transfer-matrix model for metasurface multilayer systems

2.1 Transfer matrix method for a standard homogeneous stack of layers

The optical properties of each layer are specified by its thickness d and refractive index n=ε. While our work is easily generalized to arbitrary incidence angle, here we deal with normal incidence only. At the heart of the transfer matrix method is the notion that one can relate the parallel field components E_x and H_y at the front side (z = 0) and the back side (z = d) of a slab by matrix multiplication

where M (k_d, d) is the transfer matrix. For a single slab the transfer matrix M is easily derived through the uses as auxiliary variables of forward and backward propagating fields in the slab, which are written as E_x(z) = E_fe^ikz + E_be^−ikz (suppressing time dependence e−iωt). The auxiliary field H_y is proportional to the magnetic field, with multiplicative factor iω suppressed throughout, since it is common to all H-fields independent of layer. The transfer matrix now follows by calculating H_y from E_x, and evaluating the fields at z = 0 and z = d. From the resulting equations one eliminates the auxiliary field E_f,b, to obtain

in which k_d is the wave vector nω/c in the slab with thickness d and refractive index n.

The strength of the transfer matrix method is that the transfer matrix of an arbitrary stack of layers is obtained simply as the matrix product of all the single layer transfer matrices. This fact hinges on the fact that the boundary equations for electric and magnetic fields require both E_x and H_y to be continuous across interfaces, so that the fields at the back of a stack can be simply obtained from the field at the front by multiplying with each of the single layer transfer matrices in turn. The field behind the stack of thickness d_st is thus constructed from the product of transfer matrices and the field in front of the stack using

From here on, we will denote with M_st the transfer matrix of a full stack, which we assume to be composed of layers m = 1, 2, … N, with transfer matrix M_m.

Once the transfer matrix of an arbitrary stack has been constructed, it is straightforward to calculate the complex reflection and transmission coefficients r and t by inserting as field at the front of the stack (subscript L for ‘left hand side of the slab’ or z < 0).

and at the back of the stack (R for right hand side)

This Ansatz specifies that the stack is illuminated only from the left (z < 0), and that at the other side, beyond the total thickness d_st of the stack, only a single outgoing wave exists. Evaluating these Ansatz fields at z_front = 0 and z_behind = d_st, and noting that these must be related through the stack transfer matrix, we obtain,

This equation results in explicit expressions for the complex reflection and transmission coefficients r and t of the entire stack in terms of the stack transfer matrix M_st.

2.2 Metasurface transfer matrix

Having revisited the standard transfer matrix for layered optical systems [49], [50], we now introduce the transfer matrix of a metasurface. Metasurfaces generally provide a very large, resonant, optical response that cannot be cast in a combination of a material thickness and refractive index. Instead, a metasurface is more naturally thought of as an infinitely thin sheet with a specified complex reflection and transmission coefficient r_a, t_a. We assume that these reflection and transmission coefficients for the metasurface in a homogeneous dielectric remain valid also when the metasurface is embedded in a thin dielectric slab of the same index, but inserted in a complex stack. The coefficients can thus be calculated by separate means, to serve as input for constructing a transfer matrix. In principle any full-wave method could be used. For our work we use a semi-analytical approach that constructs the response of a metasurface on the basis of scatterer polarizability, and using Ewald lattice summation to account for all the interactions in the lattice [28], [54], [55].

In order to construct a metasurface transfer matrix we specify the fields on both the left and right hand side of a metasurface of thickness 0 placed at z = 0 in terms of its complex-valued reflection and transmission coefficients. We suppose that the metasurface is simultaneously illuminated from the left by a field propagating towards the right (E_→e^ikz for z < 0) and by a field impinging on it from the right, hence propagating towards the left (E_←e^−ikz for z > 0). Accounting for the complex-valued reflection and transmission amplitudes yields the total electric field on the left E_L and right E_R.

Here we have employed the constraint that the metasurface is assumed to be embedded inside a layer with the same refractive index on either side of the metasurface, so that both r_a and t_a are independent of the incidence side. To obtain the transfer matrix, once calculates from E_L, E_R also the magnetic counterparts, and relates the fields at z = 0, evaluated on the left and right of the metasurface, by eliminating the auxiliary fields E_→,←. In matrix form this yields

and

so that the metasurface transfer matrix that satisfies (ER, HR)T=Mmeta(ra, ta)(EL, HL)T reads

For a metasurface in a homogeneous medium that is mirror symmetric in the z = 0 plane, it is necessary that t_a = 1 + r_a. It is important to note that this does not derive from energy conservation, but only from assumed mirror symmetry in the plane of the 2D lattice [56]. With this requirement, the metasurface transfer matrix becomes:

In essence this construction can be viewed as casting the S-matrix of a metasurface (relating outgoing to incident fields) into a transfer matrix (relating fields on either side of the interface to each other. This approach is hence strongly related to the metasurface S-matrix formalism of Menzel and Sperrhake et al. [51], [52].

2.3 An analytical model for the response of a simple resonant particle array

In principle, the metasurface transfer matrix M_meta (r_a, t_a) allows to insert the reflection and transmission of a metasurface as calculated by any modeling approach. In this work we focus for demonstration purposes on behavior of metasurface etalons in which the metasurface is simply a non-diffractive sheet of resonant nanoscale polarizabilities placed in a lattice, as for instance obtained by making sufficiently dense arrays of plasmon antenna particles [28], [54], [55]. To model reflection and transmission of such a layer one must first specify the polarizability of a single scatterer, then include its radiative damping to obtain a self-consistent t-matrix, and finally account for all the near-field and far-field multiple scattering interactions with its neighbors in the lattice. A didactic introduction is provided by de Abajo [54], while implementation details for arbitary lattice symmetries and polarizability tensors are listed in [28], [55]. In this work we assume that the polarizability is scalar, purely electric (p=4πεαE) and we use a unit system in which polarizability αE has units of volume. The electrostatic polarizability α0 of a single plasmonic particle as a function of frequency can be parametrized as [57]

with V a measure of scattering strength with units of volume, ω0 a resonance frequency and γ an Ohmic damping rate. This expression is exact for metallic spheres with a dielectric constant specified by a Drude model, where γ is the Ohmic damping rate and V is proportional to the scatterer physical volume. Radiative losses can be taken self-consistently into account for a single scatterer by implementing a radiation damping term [54], [57]

The resulting polarizability is known as dynamic polarizability, or alternatively as t-matrix of the point scatterer [54], [57]. This t-matrix results in self-consistent extinction, scattering and absorption cross sections that for a Drude sphere match the dipolar contribution in its Mie expansion.

Following de Abajo [54], in order to find the effective polarizability of an array of particles one needs to incorporate an Ewald lattice summation technique to account for all the near- and far-field interaction between antenna particles. Accordingly, under normal incidence, the induced dipole moment in each scatterer reads,

The term G represents a lattice Green function, i.e., a summation over the free space Green function that is generally dependent on incidence angle, lattice pitch, lattice symmetry, and polarization. While this expression thus can encode the full band structure and diffraction physics of plasmon antenna arrays [54], [55], here we focus on metasurfaces under normal incidence, that are so dense that there are no propagating grating diffraction orders. In this limit, the real part of the lattice Green function only induces a shift in resonance frequency that we incorporate into ω0, while the imaginary part adds radiative damping according to,

where A is the unit cell area [54]. It should be noted that the single-particle radiation damping is exactly canceled by the lattice sum, and replaced by a term that takes into account the super radiant collective damping that increases with antenna density. Following de Abajo [54], the lattice reflectivity is,

while the lattice transmission is t_a = 1 + r_a. Note that the lattice reflectivity converges to a perfect reflector r_a = −1 in the limit of very strong and dense scatterers (large V and small A).

2.4 Extracting local fields, induced dipole moments and dissipated power in a layer

The framework as described so far allows to calculate the complex-valued reflection and transmission of arbitrary multilayer stacks in which metasurfaces can be interspersed with dielectric and metallic planar layers, with the understanding that the two layers directly adjacent to a metasurface are chosen to have identical refractive index. One can obtain not only far field reflection and transmission but also the field inside any layer in the stack. To this end, the formalism stipulates that one must first solve for the reflection r_st and transmission t_st of the stack using the full stack transfer matrix M_st. Next, one can calculate the field at the front side of a layer m (normalized to a unit strength incident field offered to the stack as illumination) by multiplying the field on the incident side (1+rst, ik (1−rst))T with the partial stack matrix M_m−1⋅M_m−2 … M₂⋅M₁ that accounts for all layers up to m. Given the fields (Ex(zm), Hy(zm)T) at the front of the mth slab, the electric field inside it simply reads E(z)=Ef,meik(z−zm)+Eb,me−ik(z−zm) (i.e., sum of forward and backward traveling waves, defined at Eqs. (1) and (2), with

Local absorption directly follows from the field inside a layer. The absorption in layer m [energy per unit of time and unit of area lost in layer m is given by

where δm is the phase difference between the complex numbers E_f,m and E_b,m and where k′ and k″ are resp. the real and imaginary part of the wavenumber k in layer m. The first two terms can be interpreted as the incoherent sum of contributions of the forward and backward propagating waves, which are exponentially damped due to absorption. The last term arises from interference, i.e., standing wave effects in the slab. For reference, the assumed input intensity is 12cϵ0n0|E|² with n₀ the z < 0 refractive index, and |E| = 1 V/m the incident field strength.

In a similar vein, one can also determine the induced dipole moments in any metasurface inside the complex stack. Supposing that the metasurface is situated immediately after layer m in the stack, one can simply calculate the field at the front side of the metasurface by multiplying the field on the incident side (1+rst, ik(1−rst))T (with r_st the full stack reflectivity already solved for) with the partial stack matrix M_m⋅M_m−1 … M₂⋅M₁ that accounts for all layers up to the metasurface. Immediately to the left of the metasurface the field reads

Inverting Eq. (8) immediately yields the impinging fields driving the metasurface, i.e.,(E→, E←)T. The induced dipole moment per metasurface scatterer is simply given by p=αlat(E→+E←). The power absorbed by the antenna array per unit of area follows by determining the deficit between the Poynting fluxes immediately to the left and right of the metasurface

The first terms in this sum reports the absorbed power if the waves driving the metasurface from the left and right would be added incoherently. The last term accounts for interference.

3.1 Perfect absorption in metasurface etalons

In the remainder of this paper, we illustrate the power of the simple transfer matrix method to predict and understand the physics of simple metasurface stacks that are of large relevance in current research. The first example that we consider is that of a Salisbury screen, consisting of a simple metal mirror, a dielectric spacer and a resonant metasurface mirror [3], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26]. This structure can also be viewed as an etalon [23], [26], [28] in which the back reflector (metal mirror) is thick, while the front reflector is the metasurface which depending on density and frequency can range from semi-transparent to strongly reflective, and furthermore has a peculiar phase response. According to literature over the past decade, Salisbury screens can show perfect absorption in the metasurface antennas at particular resonance conditions, even if by themselves the particle arrays are only weakly absorptive. A schematic view is presented in Figure 2 as sketches on the far left. This figure presents both the response of just particle arrays, and of envisioned Salisbury screens, as function of frequency ω (vertical axis on color plots) and etalon spacing d, for different antenna array densities. We choose each scatterer to have the properties V = 6.9 ⋅ 10⁻²³ m³, ω0=2.4⋅1015 rad/s, γ=9.3⋅1013 Hz. This corresponds with the scattering properties of antennas in our previously reported experiment in which we used gold nanorod antennas [26], as verified by matching the resonance frequency, Ohmic damping, and extinction cross section (σext=4πkImαdyn≈0.088 μm2) to finite element modeling of Au nanorods in glass (circa 100 × 50 × 40 nm length, width and height). We assume the scatterers to be in a square lattice of pitch a, with a increasing from 100, 200, 350, 433 to 500 nm. The medium on the incoming side and in between mirror and metasurface is glass (n_glass = 1.45), while we assume as back reflector 50 nm of gold (n_Au = 0.25 + 4.5i, taken frequency independent). The frequency window in our plots corresponds to ca. 650–900 nm vacuum wavelength.

Figure 2:

Calculated response of Salisbury screens in reflection. From left to right the columns display metasurfaces consisting of scatterers in a square geometry with a pitch increasing from 100 to 550 nm and with fixed antenna parameters V = 6.9⋅10⁻²³ m³, ω0=2.4⋅1015 rad/s and γ=9.3⋅1013 Hz. Top row (a–e): amplitude response of the particle arrays as a function of frequency in absence of the backing mirror, showing a Lorentzian response around the resonance frequency ω0. Center (f–j) and bottom (k–o) rows show the full etalon response in amplitude and phase respectively, as a function of frequency and metasurface–mirror spacing d, employing a 50 nm Au back reflector. Reflection phase is referenced to that in absence of the array.

Figure 2a–e shows that for just a single particle array with no back reflector, the reflectivity not only strengthens, but also significantly broadens with antenna density. This is a signature of collective effects in antenna arrays, whereby super radiant line width broadening causes the lattice polarizability αlat to be significantly different from the particle response for dense lattices. For the bare particle arrays, the maximum achievable absorption is fundamentally limited by symmetry to 50%, which is reached on particle resonance near 300 nm pitch [56]. Absorption drops for larger pitch as the metasurface becomes less strongly populated and hence more transparent, and conversely also drops for smaller pitch as the metasurface becomes more reflective. For the etalon-geometry the reflectivity maps as function of frequency ω and mirror–metasurface spacing d are significantly different (Figure 2f–j). At the highest antenna density, the metasurface is practically a mirror across the entire frequency range, giving rise to etalon resonances, evident as pockets of absorption close to hyperbolas (contours of constant kd) in Figure 2f. The main signature indicating that the metasurface is not a regular mirror but a dispersive reflector is that the etalon resonances show as asymmetric spectra, and with signature vanishing right at particle resonance (2.4 ⋅ 10¹⁵ rad/s). This can be traced to the metasurface reflection phase. With decreasing antenna density the response evolves from a Fabry–Pérot-like response at large antenna density, to that of an isolated mirror with just isolated pockets of reduced reflection near antenna resonance (panel 2g). These appear only when the particle array is not in a node of the field reflected from the back mirror, i.e., in between the etalon resonances of panel Figure 2f.

Looking more closely at the low reflectivity areas in Figure 2f–i we find that there are points of identically zero reflection and hence perfect absorption. This is the well-known phenomenon of perfect absorption [3], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26], which in literature has been described as a condition of critical coupling between radiative and nonradiative loss channels of a single-port resonance. More surprising is that corresponding phase plots (Panels k–n, phase referenced to the back-reflector in absence of the particles) show singular points in the ωd-parameter space, around which the phase circulates and in which the phase is not defined. These points evidence the truly vanishing reflection observed in the corresponding amplitude plots. Moreover, it is evident that these singular points do not come alone but in fact appear as a countably finite numbers of pairs of opposite topological charge, with each pair contained between two perfect-etalon conditions (hyperbolas kd=mπ, m integer). For dense arrays, these singularities are located very close to the resonance condition for a standard etalon, though at frequency detuning that is far from the bare particle resonance. For larger pitch and hence diluted lattices, (Figure 2m and n), the singular points move away from the etalon condition, and towards the points in parameter space that are centered on the bare lattice resonance frequency ω0, and on spacer thicknesses exactly in between the etalon condition. Here the lattice is placed not at half-integer wavelengths from the back reflector, but is a quarter-wave offset, so as to be in an antinode of the standing wave generated by the back reflector. At pitches above circa 500 nm, the phase singularities annihilate, and absorption in the particles is no longer perfect (Figure 2j and o). These findings thus reproduce the well-known appearance of conditions of perfect absorption in Salisbury screens in a model that is a simple transfer matrix model. At the same time it provides insights that go beyond the simple coupled mode analysis in literature [17], [18], [19], [20], [21], [22], [23], [24] that neither foresees perfect absorption points to come in pairs, nor suggests singular phase response in parameter space. The topological nature of the phase singularities in parameter space might lead one to suspect that there is an underlying anomaly in the etalon eigen mode structure, as it is well known that, for instance, exceptional points and bound states in the continuum have topological properties. Following the general classification by Krasnok et al. [58], we assert that such physics is not at play here. According to Krasnok et al., eigen modes correspond to complex-frequency poles of the S-matrix, while the zero reflection in the system that we describe is a zero of the S-matrix on the real axis. If one analyzes the S-matrix for this system (see Eq. (19) for a simplified model) at the etalon thickness where the antenna is right in the mirrors standing wave antinode, the S-matrix is found to carry a single pole in the complex frequency plane that is accompanied by a single nearby zero. For low antenna density (no singular response), the zero is in the same half-plane as the S-matrix pole. As the antenna density is tuned, the zero reaches the real frequency axis, at which point perfect absorption appears. For larger antenna strength the zero has crossed the real frequency axis. At this point the structure does not provide its zero in reflection for the exact thickness where the antennas sit in the standing wave antinode, with reflection zeroes instead happening at adjacent spacing, i.e., in the pairs of points evident in Figure 2. It should be noted that throughout, the S-matrix pole has multiplicity one. We refer the reader to Ref. [26] for experimental evidence for the proposed appearance and annihilation of pairs of phase singularities in these structures, as function of density and in the parameter spaced spanned by frequency and etalon spacing.

3.2 Metasurface etalons with near transparent mirrors

A particular advantage of the matrix method in which one can intersperse arbitrary homogeneous layers and metasurfaces, is that one can easily explore a variety of scenarios with realistic experimental conditions that are not accessible in even more simplified models, e.g., taking into account realistic refractive indices or layer thicknesses, e.g., ones that correspond to highly ‘imperfect’ Salisbury screens. Figure 2 showed that for a near-perfect back reflector, a Salisbury screen shows phase singularities and points of perfect absorption on the proviso that the metasurface in itself is sufficiently strongly scattering. Figure 2, however, also showed that a poor meta-mirror (antenna array oscillator strength below a certain threshold) causes the singularities to disappear. This raises the question whether the use of a similarly poor back mirror could retrieve these points of singular reflectivity. Qualitatively one might argue that if perfect cancellation of reflection is due to matching of the metasurface reflection constant and that of the back reflection, then the singular response should re-appear for dilute metasurfaces, if one concomitantly reduces also the back reflector reflection constant. Figure 3a exemplifies exactly this point, for the least dense array considered in Figure 2j. While for the near perfect back reflector there are no singular points, once the mirror is replaced by one of just 5 nm thickness (partial reflectance), the pairs of points of zero reflection return, as shown in Figure 3a). These are again accompanied by phase singularities in parameter space (shown in Figure 3b), and a concomitant 2π phase swing for ω,d combinations boxed in by the singularities. A subtle point is that in this case zero reflection does not imply perfect absorption, as there is also a transmitted channel.

Figure 3:

Calculated amplitude (a) and phase (b) response for the Salisbury screen in Figure 2j and o), in which the reflectivity r_m of the back reflector has been decreased to match the reflectivity of the weakly scattering array, such that singular points in reflection are retrieved. (c and d) Graphical construction of Eq. (20) in the complex plane, showing the trajectories that the right-hand and left-hand side (gray and orange respectively) of the equation trace out as a function of frequency. (c) When a weakly scattering metasurface (pitch 550 nm) is combined with a highly reflective back reflector, the circles traced out by r_me^inkd (gray) and −ra(ω)/(1+2ra(ω)) (orange) do not intersect. (d) Reducing the back mirror reflectivity to that of a 5 nm Au layer results in crossing of the circles and retrieval of the singular points shown in (a) and (b).

An a posteriori understanding of this re-appearance of points of zero reflection and phase singular appearance can be constructed from a simple two-layer Fabry–Pérot model wherein the reflectivity of an etalon reads [26], [56]

with r_a the metasurface reflectivity, r_m the back reflector reflectivity and n the etalon spacer index. This expression implies r_FPI = 0 for the condition that

This essentially prescribes a matching condition on the metasurface reflectivity ra(ω) relative to the backreflector and etalon spacing. At the points of matching, reflection vanishes and moreover since r_FPI is a meromorphic function of ω and d, a phase singularity in ω−d parameter space accompanies each zero by the Cauchy argument principle. The left hand side r_me^2inkd circumscribes a circle of radius |r_m| centered on the origin in the complex plane as one sweeps frequency. Instead, the right hand side for a metasurface with a Lorentzian resonance circumscribes a circle exactly touching the origin in the complex plane, displaced along the real axis, and with a diameter that increases with increasing antenna oscillator strength and density. On this circle, frequencies far from resonance are close to the origin, while the resonance frequency corresponds to the point on the circle furthest from the origin. Figure 3c illustrates this behavior for a strong back reflector (|r_m| = 0.95, as in the case of the 50 nm thick mirror in Figure 2) yet a dilute metasurface (550 nm pitch). For these parameters, the metasurface is not sufficiently strongly scattering (i.e., orange circle is small) for crossings between −r_a/1 + 2r_a and r_me^2iknd to occur. One route to obtain crossings is to increase antenna density or oscillator strength, which increases the orange circle (−r_a/1 + 2r_a) in size. This explains the emergence of pairs of singularities and perfect absorption above a threshold antenna density in Figure 2. Another route instead, is to reduce the back reflector reflectivity (radius of the circle r_me^2iknd) to match the metasurface response. This is shown in Figure 3d) where r_m has been decreased to match the reflectivity of a 5 nm Au mirror. Crossings denoted by stars show that the condition in Eq. (20) is met and that, for the dilute metasurface that does not show any singular response when combined with a strong back reflector, zero reflection and phase singularities in reflection now must occur.

3.3 Birefringent Salisbury screens

In the analysis of the Salisbury screen we have so far considered polarization exactly along the main axis of the nanorod antennas. For rectangular array metasurfaces with plasmonic resonators that have their principal axis aligned with the array symmetry axis, it is straightforward to also analyze polarimetric responses, as one can separately calculate the polarization dependent responses r_x and r_y for the two principal axes from two scalar calculations. The total response to any incoming polarization can then be constructed by decomposing the input polarization in two orthogonal linear components, and coherent addition of the corresponding reflected fields. By way of example we examine the response of nanorod-array Salisbury screens, for which the strong phase and amplitude response for the polarization channel of the nanorod resonance should induce a strong linear dichroism (differential absorption between linear polarization components) and a strong linear birefringence (differential retardance). The following example highlights the potential of our transfer matrix method for the rapid analysis of metasurface-etalon based polarimetric components [9], [10], [11], [29], [30], [31], [32], [33], [51], [52].

We choose nanorod antennas to be resonant with x-polarization, such that r_x is the response previously shown (Figure 2). We assume the orthogonal antenna resonance to lie far outside the frequency window of interest so that we can model r_y as the response in absence of the antenna array (i.e., a plain mirror). We report the amplitude and phase response in the co-polarized reflection channel in panels Figure 4a, c, e, g (amplitude) and Figure 4b, d, f, h (phase) for horizontal, vertical, diagonal linear, and circular polarization. When illuminating with polarization along one of the principal axes of the system, one simply obtains a flat response for the vertical y-polarization for which the metasurface is assumed transparent, and instead the resonance with strong phase features for the horizontal x-polarization already reported in Figure 2g and l. With rotation of the polarization from horizontal to diagonal (see Figure 4c and e), the points of zero reflection disappear. This is intuitively expected since the y-polarization component is still fully reflected, irrespective of whether the x-component experiences a strong phase and amplitude effect. At the same time it is at first sight remarkable that for a band in parameter space in between the two x-polarization singular points the reflection is well below that in both the x and y channel. This is due to the fact that the x and y reflections are not in phase, leading to destructive interference. It is easy to show that the co- and cross-polarized reflectivity signature is strictly identical for diagonal and circular polarization, equaling 12[rx±ry] for co-resp. cross-polarization.

Figure 4:

(a–h) Polarization dependent amplitude and phase response of Salisbury screens in reflection, for a lattice of pitch 200 nm (other antenna parameters taken the same as in Figure 2) and a back reflector of thickness 50 nm. The top right arrows denote the input polarization, and results shown are for the co-polarized detection channel. (a and b) polarization orthogonal to the antenna resonance axis, (c and d) polarization aligned with antenna resonance axis (i.e., plots identical to Figure 2g and l), (e and f) diagonal and (g and h) right handed circular (RHC) polarization. (i–l) Ellipticity parameters ϵ and α, a measure of circular polarization and orientation respectively, for the cases of diagonal and RHC input polarization. Diagonal and RHC plots are equivalent in amplitude and phase (e–h) for the co-polarized reflected channel, but entirely distinct in ellipticity (i–l).

The linear dichroism and birefringence of nanorod-array Salisbury screens become evident when analyzing polarization conversion. We quantify the polarization state of reflected light upon diagonal and RHC polarized illumination through the polarization ellipticity ϵ and polarization ellipse major axis orientation α that fully characterize the ellipse that the electric field vector traces as a function of time. The ellipticity ϵ takes values between −1 and 1, such that ϵ is 0 for a linearly polarized field and +1 or −1 for right resp. left handed (RHC, LHC) circular polarisation. The parameter α represents the orientation angle of the major axis of the polarization ellipse, and takes values from −π/2 to π/2, where 0 encodes for horizontal orientation (along x). These parameters can be directly calculated from a field vector (E_x, E_y) as

To conclude, nanorod–antenna based mirror–metamirror etalons not only show pairs of perfect absorption and phase singularity points in their principal orientation axis, but also show a rich polarization behavior. In particular, they show very strong polarization rotation behavior due to a combination of linear dichroism and linear birefringence. This may have uses for realizing (lossy) quarter-wave and half-wave plate reflectors at specific operation points. Also since the polarimetric signature is directly related to the singular phase behavior for x-polarization, polarimetry could be used instead of interferometry to understand the phase response in x-polarization that accompanies perfect absorption. Our analysis method contributes to diversifying plasmonic color printing strategies to encode information in polarization degrees of freedom [9], [10], [11], [29], [30], [31], [32], [33], [51], [52]. Finally we note that the structures at hand might be an interesting venue to realize so-called Voigt exceptional points [59]. These correspond to singular polarimetric response and chiral eigen modes of etalons with broken cylindrical symmetry, for instance containing a birefringent medium. While in this work we discuss driven responses only, and not eigen modes, combined frequency and angle-resolved polarimetry would give access to eigen modes [59], [60].

4 Strong coupling in mirror–metamirror–mirror sandwiches

In the remainder of this paper we analyze the physics of a second highly pertinent example of a stratified metasurface stack problem. This is the physic of strong coupling between the resonances of a microcavity spanned by two mirrors, and a resonant object placed inside its field. This is a seminal problem in quantum optics [39], [61], [62] as strong coupling of a single two-level system and a high Q cavity underlies cavity QED. Also classical, and collective strong coupling is a topic of large past and current interest [39], [40], [41], [42], [43], [44], [45], [53], [63], [64], [65]. Already Ameling et al. [34], [35], [36], [37] suggested that Fabry–Pérot resonators with plasmon antenna arrays inside them should display strong coupling, and should have hybrid photonic–plasmonic resonances that have high Q, yet may have locally enhanced fields. Currently groups are exploring the combination of metasurfaces and excitonic layers inside microcavity resonators [47].

To explore the physics of such systems we examine the response of symmetric stacks where we take the scattering particle arrays from the previous section, yet now sandwiched between thicknesses d/2 of glass and 20 nm Au reflectors, effectively forming a Fabry–Pérot cavity with a metasurface in the middle. We chose the mirror thickness to tune the etalon finesse. Figure 5 plots the calculated complex reflectivity r and transmission t in amplitude, and plots the induced dipole moment at the particles. We provide both the induced dipole moment normalized to the peak dipole moment p(ω0) achieved in absence of the mirrors and on antenna resonance, and the induced dipole moment normalized to the dipole moment p(ω) achieved at the same frequency, in absence of the mirrors. Finally, the left column (panels Figure 5a and f) plots the response of the stack in absence of the scatterers, simply displaying the well-known Fabry–Pérot modes in reflection and transmission, of which the width is determined by the quality factor of the etalon through the reflectivity of the mirrors.

Figure 5:

Calculated response of scattering nanoparticle arrays sandwiched between 20 nm Au mirrors, as a function of cavity length and frequency. Nanoparticle parameters are taken the same as in Figure 2. Columns from left to right: no array, pitch 550, 350, 200 and 100 nm. The first two rows plot reflection (panels a–e) and transmission amplitude (panels f–j). Panel (g) is overplot with Fabry–Pérot resonance lines (dark blue dashed) and a line indicating the scatterer resonance frequency ω0 (yellow dashed). We take cross cuts (bright blue solid) through their points of intersection for the first and third order mode to estimate the mode splitting due to the presence of the array (see also Figure 6). The third row (panels k–n) plots the dipole moment |p| of the array, referenced by the maximum dipole moment in the bare array. The bottom row (panels o–r) plots |p| referenced by the frequency dependent |p(ω)| of the bare array. The bottom two rows are both plot between 0 and 4 for comparison.

Figure 6:

(a) Spectrum at the first order bare-etalon resonant thickness indicated as a vertical line in Figure 5g) (antenna array properties: pitch 550 nm, and V = 6.9⋅10⁻²³ m³, ω0=2.4⋅1015 rad/s, γ=9.3⋅1013 Hz as in Figure 2). Maxima denoted by blue dots are used to estimate the Rabi splitting Δω that measures the coupling strength between lattice and cavity mode. (b) Calculated splitting Δω as a function of scatterer volume V, for pitches 200, 350 and 550 nm. The blue dot corresponds to Δω from (a). (c) By plotting Δω as function of V, the expected scaling with the square root of oscillator strength is evident. (d) Δω as function of V/a2 shows that a given mode splitting is proportional to the square root of density and oscillator strength, with a proportionality constant dependent on the mode order (shown for first and third mode). In (b, c, and d), the lines serve as a guide to the eye.

Once one introduces a resonant metasurface, where we vary the pitch to vary the strength of the inserted perturbation, anticrossings appear around the resonance frequency of the scattering array (Figure 5b and g). As also observed in full-wave simulations by Ameling [34], [35], [36], [37], [38], [39], [40], [41], [42], [43], [44], [45], [46], [47], these anticrossings only appear for the symmetric modes, which have a maximum at d/2. The anti-symmetric modes have no field overlap with the particles and hence show no crossing. Increasing the density from 550 nm pitch further to 350, 200 and 100 nm pitch (Figure 5c–e and h–j), the crossings become increasingly pronounced and it becomes apparent that branches on either side of an unperturbed mode bend towards one another and merge. Qualitatively, these results match the simulations by Ameling [36], [37]. The magnitude of the induced dipole moment in the scattering array, normalised to the dipole moment that would be induced in the same scattering layer for the same incident field, but without the cavity is significantly enhanced only for a large array pitch (very low density) and for large detuning from particle resonance (intrinsically small response per antenna). The enhancement arises as a consequence of the fact that the bare etalon enhances the circulating power in the cavity, but is counteracted by collective scattering effects for dense arrays and small detunings. Indeed for dense arrays, and near resonance, the net dipole moment is comparable to, or even reduced, compared to that achievable with just the array alone. This is a direct indication that collective effects put a limit on the achievable polarization in a sheet on basis of flux arguments: the total flux that the polarized sheet radiates and/or absorbs can never exceed the flux of the input field. This should be contrasted to reports for cavity–antenna hybrids with single antennas [66], [67], [68], and to the notion of Ameling et al. that plasmon array etalons allow to boost the sensitivity of, e.g., refractive index sensors, by combining high Q with enhanced fields [34], [35], [36], [37]. It is only true for weakly scattering antennas that plasmon antenna enhancement and microcavity field enhancement effects add. For strongly scattering antennas collective effects, i.e., radiation damping limit the enhancement.

4.1 Rabi splitting

The analysis presented so far in essence reproduces the full-wave simulation results of Ameling for similar structures [34], [35], [36], [37]. The advantage of our matrix method is that parameters can be explored with ease, for instance to map the strong coupling as function of metasurface scattering strength quantitatively. A first order approach to quantify the strong coupling strength is to establish the (classical vacuum) Rabi splitting from response function spectra taken at etalon thicknesses d at which the bare etalon has its resonance frequency coincident with the particle resonance ω0 [69]. An example spectrum is shown in Figure 6a. While one could attempt to fit a full model to the transmission line shape, we take the splitting simply as the frequency difference Δω between the two maxima. This difference is plotted in Figure 6b as a function of scatterer volume V for pitches 200, 350 and 550 nm. A square root dependence on the scattering V is qualitatively apparent, and more clearly exemplified by Figure 6c that shows a linear behavior of Δω as a function of V. Qualitatively this square-root behavior is consistent with common knowledge for cavities filled with polarizable media that state that the Rabi-splitting should scale with the square of the oscillator strength of atoms in the medium [39], [40], [41].

Accordingly, one would also expect the splitting to increase with the square root of particle density, meaning an overall scaling Δω∝V/a2 (where a is lattice pitch). This expectation is verified in Figure 6d which shows that for each given mode order, the splittings as function of V/a² all collapse on a straight line. For higher order modes, the crossing reduces, as expected due to the more extended spatial mode (wider etalon) yet fixed amount of polarizable matter placed in it. The scaling of Rabi-splitting with the square root of plasmon antenna density was in fact observed in Ref. [47], where Rabi-splitting was observed in etalons with square lattices of disk-like plasmon particles at three different densities.

4.2 A metasurface is different from a dispersive medium in a cavity

At first sight, the physics of strong coupling between Fabry–Pérot modes and a particle array resonance may appear to be almost trivially identical to the signature of classical vacuum Rabi splitting obtained with linear dispersive media inserted in a cavity, as was first discussed in the seminal paper by Zhu et al. [39]. Indeed, the scaling Δω∝ρf where f is a measure for oscillator strength per scatterer and ρ for number density, is identical to that expected for an atomic gas filling an etalon (Figure 7a). That this analogy is nonetheless not straightforward is obvious firstly from the fact that the physics of a thin metasurface cannot be cast in terms of an assigned dispersive refractive index, and second from the fact that for a metasurface it would not be clear which polarizability [α0,αlat,αdyn] should be relevant to determine the anticrossing.

Figure 7:

(a) Sketch of an etalon entirely filled with a dispersive medium with refractive index n(ω). (b) Etalon in which all polarizable material from (a) is bunched into a slab of thickness d_a. (c) Etalon in which the perturbation is viewed as a strongly dispersive reflector, essentially transforming the system into two coupled cavities coupled through a partially transparent membrane in the middle. A metasurface is a strongly dispersive membrane. t_d (r_d) and t (r) are transmission (reflection) coefficients of the center mirror and the outer mirrors respectively. Importantly, both case (b) and (c) show strong coupling, but (c) cannot be described as the limit d_a → 0 of panel (b), i.e., a metasurface cannot be consistently viewed as the limit of an increasingly thin but polarizable slab.

When inserting a linear dispersive medium into a cavity (mirror spacing d, see Figure 7a), a Rabi splitting comes about because the optical round trip path length 2n(ω)d allows matching of the resonance condition more than once for the same mode order when n(ω) is dispersive. Usually this occurs in media in which only a small variation in refractive index occurs due to an atomic or excitonic resonance. A small index variation nonetheless translates to a significant phase increment due to propagation over the full thickness of the excitonic medium. For a regular dispersive medium, one can imagine a number density ρ of atoms with polarizability α filling a slab of width d_a, so that the round trip path length reads 2n(ω)d≈2(d−da)+2daRe[(1+12ρα)]. This is an approximate result that uses the fact that n2=ϵ=1+ρα can be approximated to n=1+12ρα for excitonic media where ρα≪1. The seminal work of Zhu [37] applies this with d_a = d. The linear contribution of α to the propagation delay ultimately guarantees the square root scaling of splitting with density and oscillator strength, under the assumption of a Lorentzian polarizability of the same form as chosen in Eq. (11). One could attempt to cast the physics of a thin metasurface into that of a linear dispersive medium by assuming the metasurface to also be a slab of thickness d_a filled with polarizable medium and inserted into the cavity, and calculating the propagation delay over its thickness. This is the approach suggested by Ameling [37]. Since for a metasurface the behavior is conceptually that of an infinitely thin polarizable sheet, one should then take a limit where d_a is taken to zero yet the integrated polarizability is kept constant (Figure 7b and c). In fact this limiting procedure cannot be consistently done. The round trip path length over just the hypothetically inserted slab is 2n(ω)da=2daRe[1+ρα]. Since for a metasurface the polarizability is bunched in a small thickness with one object per 2D unit cell area A (see sketch in Figure 7b) one has 2n(ω)da=2daRe[1+αAda]≈2Re[dα/A] in the limit d_a → 0, since now the polarizable matter term dominates. This immediately leads to an inconsistent scaling of Rabi splitting with polarizability at fixed small d_a, and furthermore is demonstrably inconsistent since in the limit of vanishing d_a, the resultant splitting is not independent of d_a. This appears to have been overlooked by Ameling et al. [37].

4.2.1 Meta-mirror in the middle model

Instead of treating the particle array layer as a dispersive slab that modifies phase through propagation delay, we argue that it should be viewed as a strongly dispersive reflective boundary condition in the middle of a Fabry–Pérot resonator (Figure 7c) which induces no propagation delay but a strong phase pickup upon reflection/transmission. As approximate model, we evaluate the “membrane in the middle” toy model proposed by Jayich et al. for cavity optomechanics [70]. This model assumes a Fabry–Pérot cavity with outer mirror amplitude transmission and reflection constant r, t and membrane reflection and transmission constant r_a, t_a. Solving in the specific case of a metasurface placed in the middle leads to the transmission

(22)tstack=−(ra+1)−t2eikdr2(2ra+1)e2ikd+2rraeikd−1

where we have already used t_a = 1 + r_a. In the limit of zero metasurface reflectivity this expression reverts to that for an etalon of thickness d, while for near unity reflection resonances characteristic for etalons of length d/2 appear. To analyze the appearance of strong coupling, one can use the method proposed by Zhu et al. [39]. For near-perfect end-reflectors the denominator is of the form −1+Aeiφ, showing resonances when φ=0mod2π. Therefore, one can analyze φ, which for the ‘empty’ etalon simply traces the roundtrip phase 2nkd, showing a linear behavior versus frequency crossing zero modulo 2π at each resonance. If the cavity is filled with a traditional linear dispersive medium one sets k=ωn(ω)/c, yet r_d = 0 (no metasurface), reproducing exactly the model of Ref. [39]. Instead one can treat also the case of a metasurface-in-the-middle (set k=ω/c, yet r_a as in Eq. (14)).

Figure 8a–f plot reflection, transmission and the phase φ obtained from Eq. (22) for both a traditional dispersive medium filling the cavity, and for a metasurface. As parameters we use (r, t) = (−0.92, 0.34) (equivalent to 11.5% transmission and 4% absorption) for the reflectors. For panels (a–c) we assume a dispersive atomic gas of refractive index n(ω)=1+12ρα0(ω) with ω0=2.4⋅1015 rad/s, γ=4⋅1013 s−1 and V_atom = 8 ⋅ 10⁻²⁷ m³ inserted in α0 (Eq. (11)), at density ρ=(10nm)−3 so that the maximum change in refractive Δn 0.12. For the metasurface we use r_a given by Eq. (14) with identical resonance frequency and damping as for the atomic gas, but with V = 4.8 ⋅ 10⁻²⁴ m³ and A = (200 nm)² (extinction circa 50%, lattice and gas chosen to result in the same Rabi splitting). In Figure 8a–c, the dispersive medium causes strong coupling for etalon modes of all orders. For the metasurface case such splitting is only observed for modes with a central anti-node in the cavity, commensurate with Figure 5. The phase increment φ for a bare etalon simply shows hyperbolic equiphase lines (constant kd). For a dispersive atomic gas filling the cavity, a dispersive feature around the material resonance becomes apparent, which essentially reports the dispersive, real part of refractive index n. For the assumed metasurface, however, the dispersive feature has a markedly different shape, with very sharp variations, and phase singularities.

Figure 8:

Reflection, transmission, and phase analysis of the meta-membrane in the middle model. Top row: the cavity contains a dispersive atomic medium, but no metasurface. The reflection and transmission show strong coupling for all modes. Analyzing the argument φ when identifying the denominator of Eq. (22) with −1+Aeiφ shows the dispersive medium determining the propagation phase. According to Ref. [39], etalon resonances occur as zero crossing of φ, when taking vertical cross cuts at bare-etalon resonant thicknesses. Bottom row: same for a metasurface, and no atomic medium. Only every other (symmetric) mode splits. The φ-landscape is markedly different from the top row. We assume a dispersive atomic medium with ω0=2.4⋅1015 rad/s, γ=4⋅1013 s−1 and V_atom = 8⋅10⁻²⁷ m³ inserted in α0, at density ρ=(10nm)−3. The metasurface has identical resonance frequency and damping as the atomic gas, but with V = 4.8⋅10⁻²⁴ m³ and A = (200 nm)².

Following Zhu [39], one finds the Rabi splitting by identifying the zero-crossing of φ in spectra taken at constant thickness taken right at the resonant thickness for zero detuning, i.e., taking vertical cross cuts through Figure 8c and f at d=mπcn/ω0 (m the mode order). Figure 9a and b shows as example φ for the first symmetric cavity mode for both the case of a linear dispersive homogeneous medium and a metamirror in the cavity. The phase evolves from a straight line in absence of polarizing material (V = 0), to a curve with additional zero crossings symmetrically placed around the resonance frequency ω0, and increasingly moving away from ω0 when the oscillator strength is increased. For an excitonic medium this exactly reproduces the work in Ref. [39], wherein the dispersive feature causes a modest phase advance and delay, in proportion to the small dispersive real part of the Lorentz line contribution to the refractive index (12Reρα). Linear increments in oscillator strength and density ρVatom cause the phase variation away from the V = 0-line at any given frequency to vary proportionally. The concomitant mode splitting then traces the well-known square root behavior with increasing ρVatom. For the metasurface (Figure 9b) the phase behavior is shown for densities and linearly increasing particle oscillator strengths V that correspond to identical Rabi splittings as for the atomic gas example. The phase profiles are nonetheless quite different since the phase variation near zero detuning is very large already for modest metasurface characteristics, and crosses through π. Thereby the dispersive correction to the bare etalon line is replaced by a monotonic behavior, i.e., an S-curve with a full 2π phase swing. Even if for frequencies in between the zero crossings the phase swing is large and no longer linearly proportional to V, the location of the zero crossings nonetheless still scales with the square root of oscillator strength V.

Figure 9:

The argument φ when viewing the resonant denominator of Eq. (22) as −1+Aeiφ as function of frequency, and at the etalon opening where the bare etalon has its resonance frequency at the resonance ω0=2.4⋅1015 s−1 of the inserted species. According to Ref. [39] zero-crossings correspond to the Rabi-split normal mode frequencies. Panel (a) is for a homogeneous atomic gas, and panel (b) for a dispersive metasurface. Different line colors vary oscillator strength and density in linear increments. Note how the zero crossings move away from ω0 in proportion to the square root of the oscillator strength. Finally, we have taken oscillator strengths that correspond to identical location of the zero crossing in (a) and (b). Nonetheless, the phase variation for frequencies between the crossings is strongly different. Parameters as in Figure 8, with oscillator strength V varied as indicated in the legend.

4.3 Splitting in the limit of weak polarizability

In the limit of a weakly reflective metasurface, for which rd=2πik/Aα≪1, and as long as splittings remain small compared to the resonance frequency ω0, and for electrostatic Lorentzian polarizabilities of the form α0≈ω02V/(ω02−ω2−iωγ) one can analytically derive the location of the zero crossings of φ. In the limit of zero damping, they read for a metasurface

and for a dispersive atomic gas filling the complete cavity

Although the underlying variation of φ is quite different, the splitting follows a very similar scaling at small polarizability, namely a proportionality to the square root of oscillator strength and the number of oscillators placed in the mode. An apparent difference is that for a metasurface the splitting will reduce with etalon length L (i.e., for higher order modes), while for a cavity completely filled with an excitonic medium, the splitting is instead constant. This is due to the fact that in the excitonic case a longer cavity also contains more oscillators, whereas for a metasurface the number of oscillators is fixed with varying L. Figure 10 shows color plots of the (full, unapproximated) |φ| as function of frequency and oscillator strength. The zero crossings can be traced as the contour of zero |φ| (darkblue), and correspond well with the approximate expressions Eqs. (23) and (24) over a large range of oscillator strengths.

Figure 10:

Color plots of |φ| with φ defined from equating the resonant denominator of Eq. (22) to −1+Aeiφ, as a function of oscillator strength and frequency, for fixed etalon opening chosen such that the bare etalon has its resonance frequency at the resonance ω0=2.4⋅1015 s−1 of the inserted dispersive medium (panel a) resp. metasurface (panel b). The white parabola in (a) over plots Eq. (24). Also for the resonant metasurface case, the splitting traces a parabola (white continuous line, Eq. (23). This is surprising, as this prediction takes the static polarizability as input, whereas the dynamic and lattice polarizabilities that determine the single-particle resp. collective antenna scattering strength are strongly different. Dashed and dotted lines show the Rabi splitting if one would have to correct V in Eq. (23) with the on-resonance dynamical resp. lattice radiative correction factor (see Eq. (25)). Parameters as in Figure 8.

4.4 Splitting at strong polarizability

A remarkable finding in Figure 10 is that the approximate expressions for the Rabi splitting versus oscillator strength even hold for very large oscillator strength. For the metasurface this is especially surprising, since one can ask which polarizability is actually a physical attribute of a scatterer. The polarizability that one would derive from a scattering or extinction measurement on a single scatterer, or instead from numerical full wave simulation [71], is the dynamic polarizability αdyn, and not the static polarizability α0. The polarizability that one would measure for an antenna in a lattice is in fact different from both, as the lattice polarizability αlatt is modified by super radiant, i.e., collective, radiative damping. The on-resonance lattice polarizability is smaller than the static polarizability by a factor:

As the dashed and dotted curves in Figure 10b show, if one would assume the on-resonance dynamic or lattice polarizability to determine the Rabi splitting instead of just V, Rabi splittings would be up to a factor 2 smaller. Thus Figure 10 demonstrates that even though the evaluation of the metasurface-in-the-middle model actually uses the lattice polarizability, the Rabi splitting finally traces out Eq. (23) which only contains the electrostatic oscillator strength V without the dynamic, i.e., k-dependent correction factor for the lattice. This is a remarkable finding, since the electrostatic polarizability is not actually an observable in any optical scattering experiment or full wave calculation. Our remarkable finding can be rationalized by noting that the usual statement that the on-resonance polarizability determines the magnitude of the Rabi splitting does not hold. The radiative correction is very large on resonance, but actually it is small far away from particle resonance. This makes the dynamic and lattice polarizability strongly non-Lorentzian. At the frequencies of the zero-crossings of φ the lattice polarizability is actually very close to the electrostatic polarizability at the same frequency, even though on resonance the difference is large.