In this work we present a simple transfer-matrix based modeling tool for arbitrarily layered stacks of resonant plasmonic metasurfaces interspersed with dielectric and metallic multilayers. We present the application of this model by analyzing three seminal problems in nanophotonics. These are the scenario of perfect absorption in plasmonic Salisbury screens, strong coupling of microcavity resonances with the resonance of plasmon nano-antenna metasurfaces, and the hybridization of cavities, excitons and metasurface resonances.
Since almost a decade metasurfaces have emerged as a very fast growing field of research in nanophotonics , , . In contrast to the notion of metamaterials that were designed to realize 3D bulk materials with peculiar designer permittivity and permeability, the notion of metasurfaces is that they shape waves by scattering off of an abrupt designer boundary condition with a completely controlled and locally varying amplitude and phase. This boundary condition is realized by engineering subdiffractive arrays of building blocks that are arrayed in a 2D sheet. This notion revitalized a venerable history in grating science, particularly the notion of echelette gratings to shape diffraction . Recent demonstrations of metasurfaces include a large variety of flat optics components, such as lenses , , , , waveplates and polarization optics , , , as well as holograms , , , diffusers , and even computational metasurfaces that are designed to perform simple linear mathematical operations on incident wavefronts , , . In this work we focus on 2D arrays where the constituents are strong, resonant scatterers with a dipolar plasmonic resonance, as opposed to dielectric metasurfaces , .
Understanding the physics of layered systems in which metasurfaces are stacked in succession, or interspersed with normal metal and dielectric layers, is of large relevance from several viewpoints. First, if one considers metasurfaces as objects that transform incident wavefronts according to designed mathematical functions in the context of signal processing and computation , , , an important question is how you can concatenate optical functions stacking metasurfaces. Second, seminal early examples of metasurfaces use stacking to design function. An important family of examples is that of perfect absorbers like Salisbury screens, and patch-antenna reflective metasurfaces, in which arrays of resonant antennas are placed in front of a reflector , , , , , , , , , , . Such structures control reflected waves in amplitude, phase and polarization through constructive and destructive interference of scattering contributions from a ground plane and a scatterer array , . These structures essentially constitute a metasurface etalon, where both the round trip phase, controlled by separation, and the reflection phase of the metasurface reflector, control response , , , . Beyond applications as antireflection devices, absorption enhancers in photodetection scenarios, and as metasurface pixels from which to build metasurface devices, these metasurface etalons also have been proposed for plasmonic color printing and polarization multiplexing therein , , , , . Expanding on this idea of metasurface etalons in which one reflector is a metasurface with engineered reflectivity, one can also envision more complex resonant structures, such as Fabry–Pérot resonators in which metasurfaces are inserted inside the cavity mode. These were proposed by Ameling et al. , , , ,  as high-Q plasmonic structures for sensing, in which classical analogs of vacuum Rabi splitting occurs , a hallmark of strong coupling. These types of structures are currently under intense scrutiny in the field of plasmon strong coupling and polaritonic chemistry , , , , , , , , wherein researchers construct plasmonic–photonic resonators that are in their turn coupled to excitonic media. The relevance of that field lies in the pursuit of collective light–matter interaction phenomena at room temperature, such as exciton–polariton physics, condensation and lasing, and in the pursuit of controlling chemical reaction pathways and transport properties by shifting of energy levels through vacuum Rabi splitting , , .
Full-wave simulation of stratified systems that include metasurfaces quickly becomes infeasible with increasing complexity, owing to the large overall structural dimensions of stacks, yet the nanoscale detail needed to resolve metasurface building blocks. Also, it would be inefficient to do full-wave calculations for each envisioned stack if instead there would be a simple modeling strategy in which a library of simulations for a basic set of metasurfaces could be re-used when layers are rearranged and reordered, without having to start numerical analysis from scratch. Here we present a simple modeling approach based in the standard transfer matrix method for stratified systems , , yet capable to include also metasurfaces. To this end we rearrange the scattering S-matrices of metasurfaces, i.e., the input–output relations that can be extracted from any other suited and dedicated simulation tool, into transfer matrix method form. This approach is similar to the work of Menzel et al. and Sperrhake et al., who developed S-matrix based models to analyze the polarimetric behavior of metasurface stacks , .
This paper is structured as follows. We first present an overview of how to extend the simple transfer matrix method to include metasurfaces, and indicate how one should extract relevant observables from it. Next we show that this simple model gives surprising new insights in Salisbury screens for perfect absorption , , , , , , , , , , , as well as in the problem of strong coupling of microcavity modes with metasurface resonances , , , , . Finally, we discuss the application of our modeling approach to strong coupling of metasurface etalons with excitonic media , [ 53].
2 Transfer matrix model for multilayer stack with interspersed metasurfaces
Figure 1 shows a sketch of the family of systems that we model in this work, which consists of an arbitrary stack of layers with parallel interfaces. The layers can be either standard homogeneous dielectric or metallic materials, or instead metasurfaces embedded inside a dielectric layer. For completeness, we first briefly introduce the standard transfer matrix method for layered systems, as is described at length in the seminal paper by Yeh, Yariv and Hong , [ 50].
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 . 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 Ex and Hy at the front side (z = 0) and the back side (z = d) of a slab by matrix multiplication
where M (kd, 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 Ex(z) = Efeikz + Ebe−ikz (suppressing time dependence ). The auxiliary field Hy is proportional to the magnetic field, with multiplicative factor suppressed throughout, since it is common to all H-fields independent of layer. The transfer matrix now follows by calculating Hy from Ex, and evaluating the fields at z = 0 and z = d. From the resulting equations one eliminates the auxiliary field Ef,b, to obtain
in which kd is the wave vector 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 Ex and Hy 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 dst 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 Mst the transfer matrix of a full stack, which we assume to be composed of layers m = 1, 2, … N, with transfer matrix Mm.
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 dst of the stack, only a single outgoing wave exists. Evaluating these Ansatz fields at zfront = 0 and zbehind = dst, 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 Mst.
2.2 Metasurface transfer matrix
Having revisited the standard transfer matrix for layered optical systems , , 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 ra, ta. 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 , , .
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→eikz 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 EL and right ER.
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 ra and ta are independent of the incidence side. To obtain the transfer matrix, once calculates from EL, ER 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
so that the metasurface transfer matrix that satisfies reads
For a metasurface in a homogeneous medium that is mirror symmetric in the z = 0 plane, it is necessary that ta = 1 + ra. 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 . 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. , .
2.3 An analytical model for the response of a simple resonant particle array
In principle, the metasurface transfer matrix Mmeta (ra, ta) 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 , , . 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 , while implementation details for arbitary lattice symmetries and polarizability tensors are listed in , . In this work we assume that the polarizability is scalar, purely electric () and we use a unit system in which polarizability has units of volume. The electrostatic polarizability of a single plasmonic particle as a function of frequency can be parametrized as 
with V a measure of scattering strength with units of volume, 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 , 
The resulting polarizability is known as dynamic polarizability, or alternatively as t-matrix of the point scatterer , . 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 , 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 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 , , 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 , while the imaginary part adds radiative damping according to,
where is the unit cell area . 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 , the lattice reflectivity is,
while the lattice transmission is ta = 1 + ra. Note that the lattice reflectivity converges to a perfect reflector ra = −1 in the limit of very strong and dense scatterers (large V and small ).
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 rst and transmission tst of the stack using the full stack transfer matrix Mst. 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 with the partial stack matrix Mm−1⋅Mm−2 … M2⋅M1 that accounts for all layers up to m. Given the fields at the front of the mth slab, the electric field inside it simply reads (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 is the phase difference between the complex numbers Ef,m and Eb,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 |E|2 with n0 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 (with rst the full stack reflectivity already solved for) with the partial stack matrix Mm⋅Mm−1 … M2⋅M1 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.,. The induced dipole moment per metasurface scatterer is simply given by . 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 Salisbury screen as a metasurface etalon
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 , , , , , , , , , , . This structure can also be viewed as an etalon , ,  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−23 m3, , . This corresponds with the scattering properties of antennas in our previously reported experiment in which we used gold nanorod antennas , as verified by matching the resonance frequency, Ohmic damping, and extinction cross section 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 (nglass = 1.45), while we assume as back reflector 50 nm of gold (nAu = 0.25 + 4.5i, taken frequency independent). The frequency window in our plots corresponds to ca. 650–900 nm vacuum wavelength.
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 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 . 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 ⋅ 1015 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 , , , , , , , , , , , 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 -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 , 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 , 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 , , , , , , ,  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. , 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.  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 phase swing for 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.
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 , 
with ra the metasurface reflectivity, rm the back reflector reflectivity and n the etalon spacer index. This expression implies rFPI = 0 for the condition that
This essentially prescribes a matching condition on the metasurface reflectivity relative to the backreflector and etalon spacing. At the points of matching, reflection vanishes and moreover since rFPI is a meromorphic function of ω and d, a phase singularity in parameter space accompanies each zero by the Cauchy argument principle. The left hand side rme2inkd circumscribes a circle of radius |rm| 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 (|rm| = 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 −ra/1 + 2ra and rme2iknd to occur. One route to obtain crossings is to increase antenna density or oscillator strength, which increases the orange circle (−ra/1 + 2ra) 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 rme2iknd) to match the metasurface response. This is shown in Figure 3d) where rm 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 rx and ry 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 , , , , , , , , , .
We choose nanorod antennas to be resonant with x-polarization, such that rx 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 ry 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 for co-resp. cross-polarization.
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 to , where 0 encodes for horizontal orientation (along x). These parameters can be directly calculated from a field vector (Ex, Ey) 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 , , , , , , , , , . Finally we note that the structures at hand might be an interesting venue to realize so-called Voigt exceptional points . 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 , .
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 , ,  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 , , , , , , , , , , . Already Ameling et al. , , ,  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 .
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 achieved in absence of the mirrors and on antenna resonance, and the induced dipole moment normalized to the dipole moment 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.
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 , , , , , , , , , , , , , , 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 , . 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 , , , 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 , , , . 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 , , , . 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 . 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 . 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 , , .
Accordingly, one would also expect the splitting to increase with the square root of particle density, meaning an overall scaling (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/a2 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. , 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. . Indeed, the scaling 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  should be relevant to determine the anticrossing.
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 allows matching of the resonance condition more than once for the same mode order when 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 da, so that the round trip path length reads . This is an approximate result that uses the fact that can be approximated to for excitonic media where . The seminal work of Zhu  applies this with da = 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 da filled with polarizable medium and inserted into the cavity, and calculating the propagation delay over its thickness. This is the approach suggested by Ameling . Since for a metasurface the behavior is conceptually that of an infinitely thin polarizable sheet, one should then take a limit where da 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 . Since for a metasurface the polarizability is bunched in a small thickness with one object per 2D unit cell area (see sketch in Figure 7b) one has in the limit da → 0, since now the polarizable matter term dominates. This immediately leads to an inconsistent scaling of Rabi splitting with polarizability at fixed small da, and furthermore is demonstrably inconsistent since in the limit of vanishing da, the resultant splitting is not independent of da. This appears to have been overlooked by Ameling et al. .
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 . This model assumes a Fabry–Pérot cavity with outer mirror amplitude transmission and reflection constant r, t and membrane reflection and transmission constant ra, ta. Solving in the specific case of a metasurface placed in the middle leads to the transmission
where we have already used ta = 1 + ra. 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. . For near-perfect end-reflectors the denominator is of the form , showing resonances when . Therefore, one can analyze φ, which for the ‘empty’ etalon simply traces the roundtrip phase 2nkd, showing a linear behavior versus frequency crossing zero modulo at each resonance. If the cavity is filled with a traditional linear dispersive medium one sets , yet rd = 0 (no metasurface), reproducing exactly the model of Ref. . Instead one can treat also the case of a metasurface-in-the-middle (set , yet ra 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 with , and Vatom = 8 ⋅ 10−27 m3 inserted in (Eq. (11)), at density so that the maximum change in refractive Δn 0.12. For the metasurface we use ra given by Eq. (14) with identical resonance frequency and damping as for the atomic gas, but with V = 4.8 ⋅ 10−24 m3 and A = (200 nm)2 (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.
Following Zhu , 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 (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 , and increasingly moving away from when the oscillator strength is increased. For an excitonic medium this exactly reproduces the work in Ref. , 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 (). Linear increments in oscillator strength and density 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 . 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 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.
4.3 Splitting in the limit of weak polarizability
In the limit of a weakly reflective metasurface, for which , and as long as splittings remain small compared to the resonance frequency , and for electrostatic Lorentzian polarizabilities of the form 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.
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 , is the dynamic polarizability , and not the static polarizability . The polarizability that one would measure for an antenna in a lattice is in fact different from both, as the lattice polarizability 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.
5 Excitonic materials combined with resonant metasurfaces in etalons
There is currently a large interest in strong coupling of Fabry–Pérot resonators with excitons in 2D transition metal dichalcogenide (TMDC) materials, large oscillator strength organic systems like J-aggregates, as well as various types of 2D and 3D semiconductors , , , , , , , , , , . The rationale is that strong coupling gives access to, for instance, strongly nonlinear photonics via excitonic nonlinearities, which in turn gives rise to exciting opportunities for bistable optical devices, nonlinear sensing, and classical and quantum hardware optical Ising simulators , . Recent experimental reports include strong coupling of microcavities and plasmonic arrays with organic dye molecules , , , , extremely dense organic molecular ensembles that form J-aggregates , and for instance 2D transition metal dichalcogenide flakes , . These recent reports in fact build on experiments pioneered by Weisbuch  (III–V microcavities and quantum wells) and Lidzey  (microcavities with organic excitonic layers).
For decades, the transfer matrix has been used to predict classical strong coupling in microcavities filled with J-aggregates, quantum wells, or other excitonic materials , , , , . As model for the dielectric constant of such layers we use
The model accounts for a background refractive index nbg, and reparametrizes oscillator strength such that Δn is approximately the magnitude of the complex refractive index change on resonance. Figure 11 shows example etalon transmission for J-aggregates (, and background index nbg = 1.55, (values commensurate with Ref.  for TDBC in a PVA, with adjusted resonance frequency to match the plasmon antennas assumed throughout in this work) and for WS2 flakes (, , consistent with Ref. ). Strong coupling is predicted in both systems at sufficient layer thickness, where the stronger exciton line of WS2 means that a smaller thickness is required. The J-aggregate system with these parameters suffers from significant loss (broadening of bands obscures anticrossing). Peculiar to the WS2 system is the fact that the very large constant background index of the semiconductor significantly changes the mode frequencies even in absence of the excitonic line (due to nbg only). For small thickness the antisymmetric modes (node at the excitonic layer) hardly shift, but the resonance conditions for even ones are dramatically changed. This re-orders the modes at modest thickness (100 nm), and modifies the mode spacing for full filling.
Given that excitonic matter can strongly couple to resonances of metasurface arrays, there is interest in understanding what the optical signatures are when etalons are constructed that contain both excitonic media and plasmon particle arrays. Bisht et al.  recently reported measurements of Rabi splitting in cavities with simultaneously plasmonic and excitonic materials loaded in them, reporting observation of a joint Rabi splitting even exceeding the summed Rabi splitting from just particles, and just excitonic material alone. We present absorption calculations, as absorption partitioned over its plasmonic and excitonic contribution also report on the degree to which each type of polariton contributes to the hybrid plasmonic–excitonic polariton. As parameters we assume a material similar to the J-aggregate (200 nm slab), and a metasurface array with antennas also as before, at a pitch 350 nm. We take the mirrors to be non-absorbing (setting for gold purely negative) so that absorption calculations strictly trace plasmonic and excitonic loss channels. Figure 12a shows that if the etalon, which we take to have lossless mirrors, is filled with just a plasmon particle array, then absorption as function of frequency and etalon spacing simply traces the dispersion of the anticrossing even etalon modes (mode maximum at the etalon center). The odd modes are absent in absorption, even if they do appear as transmission resonances, as they have an antinode at the absorbing particles. If instead a slab of excitonic material is placed in the middle of the etalon, both the even and odd modes show in absorption (Figure 12b) since the slab thickness required for strong coupling means that also the anti-symmetric etalon modes have overlap with it. Finally, the absorption spectrum of the joint system is shown in Figure 12c. Evidently, the even modes now show a wider anticrossing then is the case for either excitonic, or plasmonic system alone. This is qualitatively consistent with the experimental report of Bisht .
To more quantitatively assess the anticrossing, we examine absorption spectra at the first-order bare-etalon resonant thickness in Figure 13, where we furthermore separate out absorption in the metasurface (panel a) and in the excitonic material (panel b). For reference also the absorption is shown for either of the layers alone, in presence and in absence of surrounding mirrors. The absorption spectrum of the individual plasmonic and exitonic layer is a single peak, though not Lorentzian in shape as at the assumed oscillator strength extinction is strong. Upon introduction of the mirrors, the peak in both cases splits, indicating that strong coupling occurs for the case of etalon and metasurface, and for the case of etalon and excitonic material (orange curves). Coincidental to the choice of parameters is that the splitting in both systems is similar in magnitude. In the joint plasmonic–excitonic system the apparent Rabi splitting is larger than in either system separately, but not larger than the sum of Rabi splittings. Moreover, it is noteworthy that in both subsystems the absorption maxima occur at the same frequencies, indicating that these are the eigenfrequencies of the tripartite hybrid modes. A first order approximation to the Rabi splitting Ωmeta+exciton predicted by the membrane-in-the-middle model is
in terms of the Rabi splittings Ωmeta and Ωexciton of the etalon with just the metasurface alone, and with just the excitonic material. Thus it appears that as reported by Bisht et al.  a large Rabi splitting is indeed easier to achieve by combining excitonic and plasmonic constituents in an etalon, even if in distinction to  the joint Rabi splitting cannot exceed the sum of Rabi splittings. Unfortunately, this larger splitting would not necessarily be of help for envisioned scenarios in cavity exciton–polariton physics, where one usually seeks to imbue photons with nonlinearities by casting a significant fraction of the excitation into the intrinsically nonlinear exciton. Figure 13b reveals that in the tripartite system the absorption in the upper and lower plasmon–exciton–polariton branch is not larger than that in case of the exciton-only etalon. An interesting question outside the scope of our transfer matrix model is in how far near-field enhancement, which is poorly represented in our model, contributes to further enlarged Rabi splitting and the exciton fraction, as is suggested by the work of Bisht .
6 Domain of validity and benchmark
While the metasurface multilayer transfer model as we sketched it already applies to many interesting problems, it is important to demarcate its validity and provide a benchmark. The model as we sketched it is specific to normal incidence, and, as it is scalar, does not include polarization effects. This is only a valid approach in case all the plasmon arrays in the stack have the same principal polarization axis, and in case the stack as a whole is interrogated with a polarization coincident with such a principal axis. As we showed for the particular case of the birefringent Salisbury screen, in case the antennas are, for instance, rectangular and aligned with the principal axis of rectangular arrays, one can trivially extend the formalism to deal with polarization-dependent scattering phenomena by decomposing the incident polarization in the two principal polarization components, for which separately the formalism applies. The model can also be easily extended to more complicated scenarios, such as off-normal incidence, metasurfaces with intrinsically magneto-electric or chiral unit cells , , , or when metasurfaces are stacked in a twisted manner. Since one must then account for two polarizations, this requires a 4×4 transfer matrix approach to allow for the polarization cross coupling effects. This would bring the proposed method on par with the S-matrix multilayer approach for metasurface stacks developed by Menzel and Sperrhake et al., who focused specifically on polarimetric applications , . A more fundamental limitation of our approach is that it cannot deal with propagating and evanescent diffraction orders. This implies, firstly, that grating diffraction phenomena, such as surface lattice resonances and waveguide lattice resonances that are reliant on grating diffraction orders that are propagative in at least one of the layers are beyond the scope of our model. Evanescent diffraction orders are required to describe near-field hybridization effects, which become relevant when particle arrays come to within a fraction of a wavelength () from interfaces or other metasurfaces . Here we discern two distinct cases. First, many metasurface realizations deposit scatterers on a dielectric interface, e.g., in air, but on a glass substrate. In this case a practical approach is to model the antennas as being positioned at the location of their center of mass (slightly away from an interface), and adjusting the polarizability V, and γ to match full wave calculations. This applies because for moderate index contrasts between either side of the interface, the presence of the interface only slightly renormalizes the polarizability but without changing the resonance character of the antennas. Instead, as a second extreme one could consider the applicability of our model to patch antenna arrays, i.e., antenna arrays at deep sub wavelength separation from metal interfaces. The strong near-field interaction then fundamentally changes the resonance character of the antennas, and the model does not apply. Even within the dipole approximation, our model ignores near-field interactions with interfaces that are encode in evanescent diffraction orders. In addition near-fields calculated in this model are only of qualitative use, as the model assumes dipolar scatterers instead of the real particle geometry.
To illustrate some of these limitations, Figure 14 provides finite element, i.e., full-wave, simulations for select examples from Figure 2 (Salisbury screen) and Figure 5 (strong coupling). These calculations use identical parameters for mirrors and glass layers as in the multilayer transfer matrix approach, and as particles assume rectangular nanorods (88 by 40 nm, with a height of 20 nm), and a Drude model for the antenna dielectric constant of the form with , and . The particle size was chosen to match the bare lattice response at pitch of 450 nm in COMSOL with that given by Eq. (14) (resonance frequency, peak reflectivity, transmission and absorption). The benchmark uses symmetry planes to simulate just one quarter of the unit cell, and measures input/output through ports located at either end of the stack (somewhat away from the mirror interfaces). Figure 14 shows that the full wave benchmark calculations confirm all our predictions regarding the peculiar amplitude and phase response of Salisbury screens, the appearance and magnitude of strong coupling in etalons, and the progression with antenna density. At the same time the benchmark calculations provide instructive insights in the limitations of the method, where there are three main effects at play. First, the resonance feature slightly blue shifts with increasing lattice density, as a consequence of dipole–dipole interactions between antennas (guides to the eye to indicate bare lattice resonance incl. particle interactions indicated by horizontal lines in graphs, labeled I). These are in fact contained in the real part of the lattice Green function in Eq. (14) which renormalizes the real part of the resonance frequency. While for clarity of presentation we set this real part to zero for our illustration examples so as to be able to compare the salient features with strictly aligned resonances, they in fact can be accounted for in lattice sum theory. Second, at distances well below 100 nm between particles and mirror near-field interactions set in that strongly red shift the antenna resonance frequency (features labeled with letter “II” in graphs). Description of this effect both requires evanescent grating orders that are not contained in our model, and would require to go beyond the dipole approximation. Since neither evanescent grating orders nor microscopic near-field detail is contained in the transfer matrix approach, this highlights a limitation on the domain of validity. Finally, the third effect is that for the lowest antenna density the large pitch allows launching of real grating diffraction into glass, as well as grating-vector assisted coupling into a variety of resonances that include the symmetric and antisymmetric surface plasmon polariton modes supported by the metal mirror, as well as grazing angle etalon resonances in the spacer between mirror and particles. This explains the sharp features that occur for the largest pitch only (labeled with “III” in the graphs). To include such grating diffraction in a dipole antenna scattering model in stratified systems is possible, but far from easy to implement . Finally we note that despite the fact that these are very small COMSOL models (just of order 5⋅104 elements) that take just between 3 and 10 s to evaluate (Intel I7-6950X 10 cores, running at 3.3 GHz), the total evaluation time for Figure 14 is of order 18 h for the Salisbury screen case (leftmost three panels, ca. 6000 parameter combinations per panel), and over three times that for the strong coupling. This should be contrasted to the multilayer approach, where it takes on order of minutes to simulate the bare lattice reflectivity in COMSOL that is input in our model, and which suffices to generate all the graphs.
To conclude, we have presented a simple yet powerful transfer matrix model for analyzing arbitrary stacks of dielectric and metallic layers combined with metasurfaces that can be inserted in any of the dielectric layers. The power of this approach is that once a library of metasurface reflection and transmission coefficients has been calculated, complex structures can be analyzed easily. We have demonstrated the power of this approach to gain insight into topical nanophotonic problems on basis of three examples. The first is the response of Salisbury screens, i.e., etalons in which one mirror is replaced by a resonant metasurface. The phenomenon of perfect absorption in such Salisbury screens actually has a deeper origin as a topological property in parameter space that is evident as pairs of phase singularities. These singularities in turn imply strong linear dichroism and linear birefringence properties, if one elicits the Salisbury screen response in just one resonant polarization axis. The second example that we analyzed concerns strong coupling between plasmonic metasurfaces and etalon resonances. While scaling laws for the magnitude of the splitting in function of oscillator strength and density are similar to those for classical models of vacuum Rabi splitting in a cavity filled with an excitonic medium, the underlying physics is different. Indeed, in metasurface etalons, the metasurface acts as a dispersive version of the ‘membrane in the middle’ problem with a strong phase response due to its interface reflectivity and transmission phase, and not due to accumulation of propagation phase. Finally, we have explored the combination of plasmonic and excitonic media in etalons, which provide Rabi splittings beyond the values possible with just either of the constituents.
The model that we proposed extends far beyond the three simple examples that we analyzed. As already foreseen by Menzel and Sperrhake ,  a rich variety of polarization conversion phenomena with twisted chiral and achiral metasurfaces can be readily analyzed by a simple expansion of our model. Also, one can envision the analysis of stacks of metasurfaces in the context of computational optics, where each metasurface may execute a mathematical function. Of large importance for this agenda is to deal with space-varying amplitude, phase and polarization responses in the plane of each metasurface. In current state of the art metasurface design, space varying metasurfaces are often built using fixed lattice sites, but individually varying scatterers, chosen on basis of library simulations for purely periodic systems. For structures discussed here we anticipate that similarly our results apply as a library to construct the space varying response of space-varying metasurfaces if they would be tessellated over tiles over ∼5 unit cells in size. In how far this approach would hold also for vertically much more extended structures (widely separated layers) is an open question. Our model can also deal with a variety of light–matter interaction scenarios in which stacks of metasurfaces interact with layers of excitonic media, of large interest in the field of collective strong coupling of organic molecules to metasurfaces and etalons. The approach can be easily extended to off-normal propagation and magneto-electric (chiral, non-reciprocal) scattering phenomena, as well as inclusion of gain alongside loss. A drawback of the model is the neglect of grating diffraction orders, both of propagating and evanescent nature. This implies that microscopic details on near-fields and near-field hybridization between very close layers and interfaces are neglected. Nonetheless our model can provide crucial guidance to researchers, helping to delineate which observed phenomena originate truly from hot spots and near-field hybridization effects, and which observed phenomena originate strictly from the zeroth-order phase and amplitude response of metasurfaces, in their multiple scattering interaction with dielectric stacks.
Funding source: Nederlandse Organisatie voor Wetenschappelijk Onderzoek
Award Identifier / Grant number: 680.47.621
This work is part of the research programme Hybrid nanophotonic architectures for ultrafast quantum optics [NWO-Vici] with project number 680.47.621, which is financed by the Dutch Research Council (NWO). The research was performed at the NWO research institute AMOLF. The authors are grateful to the Resonant Nanophotonics group at AMOLF, as well as to Randall Goldmith, for constructive feedback, and to Radek Kolkowski for help in 3D rendering.
Author contribution: All the authors have accepted responsibility for the entire content of this submitted manuscript and approved submission.
Research funding: This research was funded by Nederlandse Organisatie voor Wetenschappelijk Onderzoek (No. 680.47.621).
Conflict of interest statement: The authors declare no conflicts of interest regarding this article.
 P. Genevet, F. Capasso, F. Aieta, M. Khorasaninejad, and R. Devlin, “Recent advances in planar optics: from plasmonic to dielectric metasurfaces,” Optica, vol. 4, no. 1, pp. 139–152, 2017, https://doi.org/10.1364/optica.4.000139. Search in Google Scholar
 F. Ding, Y. Yang, R. A. Deshpande, and S. I. Bozhevolnyi, “A review of gap-surface plasmon metasurfaces: fundamentals and applications,” Nanophotonics, vol. 7, pp. 1129–1156, 2018, https://doi.org/10.1515/nanoph-2017-0125. Search in Google Scholar
 P. Lalanne, S. Astilean, P. Chavel, E. Cambril, and H. Launois, “Blazed binary subwavelength gratings with efficiencies larger than those of conventional échelette gratings,” Opt. Lett., vol. 23, no. 14, pp. 1081–1083, 1998, https://doi.org/10.1364/ol.23.001081. Search in Google Scholar
 F. Aieta, P. Genevet, M. A. Kats, N. Yu, R. Blanchard, Z. Gaburro, et al., “Aberration-free ultrathin flat lenses and axicons at telecom wavelengths based on plasmonic metasurfaces,” Nano. Lett., vol. 12, no. 9, pp. 4932–4936, 2012, https://doi.org/10.1021/nl302516v. Search in Google Scholar
 A. Arbabi, Y. Horie, A. J. Ball, M. Bagheri, and A. Faraon, “Subwavelength-thick lenses with high numerical apertures and large efficiency based on high-contrast transmitarrays,” Nat. Commun., vol. 6, no. 1, p. 7069, 2015, https://doi.org/10.1038/ncomms8069. Search in Google Scholar
 M. Khorasaninejad, W. T. Chen, R. C. Devlin, J. Oh, A. Y. Zhu, and F. Capasso, “Metalenses at visible wavelengths: Diffraction limited focusing and subwavelength resolution imaging,” Science, vol. 352, no. 6290, pp. 1190–1194, 2016, https://doi.org/10.1126/science.aaf6644. Search in Google Scholar
 T. Y. Huang, R. R. Grote, S. A. Mann, D. A. Hopper, A. L. Exarhos, G. G. Lopez, et al., “A monolithic immersion metalens for imaging solid-state quantum emitters,” Nat. Commun., vol. 10, 2019, Art no. 2392. https://doi.org/10.1038/s41467-019-10238-5. Search in Google Scholar
 N. Yu, F. Aieta, P. Genevet, M. A. Kats, Z. Gaburro, and F. Capasso, “A broadband, background-free quarter-wave plate based on plasmonic metasurfaces,” Nano. Lett., vol. 12, no. 12, pp. 6328–6333, 2012, https://doi.org/10.1021/nl303445u. Search in Google Scholar
 A. Arbabi, Y. Horie, M. Bagheri, and A. Faraon, “Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission,” Nat. Nanotechnol., vol. 10, no. 11, pp. 937–943, 2015, https://doi.org/10.1038/nnano.2015.186. Search in Google Scholar
 M. Khorasaninejad, W. Zhu, and K. B. Crozier, “Efficient polarization beam splitter pixels based on a dielectric metasurface,” Optica, vol. 2, no. 4, pp. 376–382, 2015, https://doi.org/10.1364/optica.2.000376. Search in Google Scholar
 H. Kwon, E. Arbabi, S. M. Kamali, M. Faraji-Dana, and A. Faraon, “Computational complex optical field imaging using a designed metasurface diffuser,” Optica, vol. 5, pp. 924–931, 2018, https://doi.org/10.1364/optica.5.000924. Search in Google Scholar
 A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” Science, vol. 343, no. 6167, pp. 160–163, 2014, https://doi.org/10.1126/science.1242818. Search in Google Scholar
 H. Kwon, D. Sounas, A. Cordaro, A. Polman, and A. Alù, “Nonlocal metasurfaces for optical signal processing,” Phys. Rev. Lett., vol. 121, no. 17, pp. 1–6, 2018, Art no. 173004, https://doi.org/10.1103/physrevlett.121.173004. Search in Google Scholar
 A. Cordaro, H. Kwong, D. Sounas, A. F. Koenderink, A. Alù, and A. Polman, “High-index dielectric metasurfaces performing mathematical operations,” Nano. Lett., vol. 19, no. 12, pp. 8418–8423, 2019, https://doi.org/10.1021/acs.nanolett.9b02477. Search in Google Scholar
 B. Sima, K. Chen, X. Luo, J. Zhao, and Y. Feng, “Combining frequency-selective scattering and specular reflection through phase-dispersion tailoring of a metasurface,” Phys. Rev. Appl., vol. 10, 2018, Art no. 064043, https://doi.org/10.1103/physrevapplied.10.064043. Search in Google Scholar
 W. W. Salisbury, Absorbent body for electromagnetic waves. US Patent No. 1952;p. 2,599,944. Search in Google Scholar
 R. L. Fante and M. T. McCormack, “Reflection properties of the Salisbury screen,” IEEE Trans. Antenn. Prop., vol. 36, pp. 1443–1454, 1988, https://doi.org/10.1109/8.8632. Search in Google Scholar
 A. Pors and S. I. Bozhevolnyi, “Plasmonic metasurfaces for efficient phase control in reflection,” Opt. Express, vol. 21, no. 22, pp. 27438–27451, 2013, https://doi.org/10.1364/oe.21.027438. Search in Google Scholar
 R. Alaee, M. Farhat, C. Rockstuhl, and F. Lederer, “A perfect absorber made of a graphene micro-ribbon metamaterial,” Opt. Express, vol. 20, pp. 28017–28024, 2012, https://doi.org/10.1364/oe.20.028017. Search in Google Scholar
 F. Huang, S. Drakeley, M. Millyard, A. Murphy, R. White, E. Spigone, et al., “Zero-reflectance metafilms for optimal plasmonic sensing,” Adv. Opt. Mater., vol. 4, pp. 328–335, 2015, https://doi.org/10.1002/adom.201500424. Search in Google Scholar
 R. Alaee, M. Albooyeh, and C. Rockstuhl, “Theory of metasurface based perfect absorbers,” J. Phys. D: Appl. Phys., vol. 50, no. 50, p. 503002, 2017, https://doi.org/10.1088/1361-6463/aa94a8. Search in Google Scholar
 P. T. Bowen, A. Baron, and D. R. Smith, “Theory of patch-antenna metamaterial perfect absorbers,” Phys. Rev. A, vol. 93, 2016, Art no. 063849, https://doi.org/10.1103/physreva.93.063849. Search in Google Scholar
 K. Aydin, V. E. Ferry, R. M. Briggs, and H. A. Atwater, “Broadband polarization-independent resonant light absorption using ultrathin plasmonic super absorbers,” Nat. Commun., vol. 2, p. 517, 2011, https://doi.org/10.1038/ncomms1528. Search in Google Scholar
 A. Berkhout and A. F. Koenderink, “Perfect absorption and phase singularities in plasmon antenna array etalons,” ACS Photonics, vol. 6, no. 11, pp. 2917–2925, 2019, https://doi.org/10.1021/acsphotonics.9b01019. Search in Google Scholar
 K. Kumar, H. Duan, R. S. Hegde, S. C. W. Koh, J. N. Wei, and J. K. W. Yang, “Printing colour at the optical diffraction limit,” Nat. Nanotechnol., vol. 7, no. 9, pp. 557–561, 2012, https://doi.org/10.1038/nnano.2012.128. Search in Google Scholar
 A. Kwadrin, C. I. Osorio, and A. F. Koenderink, “Backaction in metasurface etalons,” Phys. Rev. B, vol. 93, 2016, Art no. 104301, https://doi.org/10.1103/physrevb.93.104301. Search in Google Scholar
 S. J. Tan, L. Zhang, D. Zhu, X. M. Goh, Y. M. Wang, K. Kumar, et al., “Plasmonic color palettes for photorealistic printing with aluminum nanostructures,” Nano. Lett., vol. 14, no. 7, pp. 4023–4029, 2014, https://doi.org/10.1021/nl501460x. Search in Google Scholar
 Y. Gu, L. Zhang, J. K. W. Yang, S. P. Yeo, and C. W. Qiu, “Color generation via subwavelength plasmonic nanostructures,” Nanoscale, vol. 7, pp. 6409–6419, 2015, https://doi.org/10.1039/c5nr00578g. Search in Google Scholar
 A. Kristensen, J. Yang, S. Bozhevolnyi, S. Link, P. Nordlander, N. Halas, et al., “Plasmonic colour generation,” Nat. Rev. Mater., vol. 2, p. 16088, 2016, https://doi.org/10.1038/natrevmats.2016.88. Search in Google Scholar
 M. Song, X. Li, M. Pu, Y. Guo, K. Liu, H. Yu, et al., “Color display and encryption with a plasmonic polarizing metamirror,” Nanophotonics, vol. 7, pp. 323–331, 2017, https://doi.org/10.1515/nanoph-2017-0062. Search in Google Scholar
 M. Song, Z. A. Kudyshev, H. Yu, A. Boltasseva, V. M. Shalaev, and A. V. Kildishev, “Achieving full-color generation with polarization-tunable perfect light absorption,” Opt. Mater. Express, vol. 9, no. 2, pp. 779–787, 2019, https://doi.org/10.1364/ome.9.000779. Search in Google Scholar
 R. Ameling, L. Langguth, M. Hentschel, M. Mesch, P. Braun, and H. Giessen, “Cavity-enhanced localized plasmon resonance sensing,” Appl. Phys. Lett., vol. 97, p. 253116, 2010, https://doi.org/10.1063/1.3530795. Search in Google Scholar
 R. Ameling, D. Dregely, and H. Giessen, “Strong coupling of localized and surface plasmons to microcavity modes,” Opt. Lett., vol. 36, pp. 2218–20, 2011, https://doi.org/10.1364/ol.36.002218. Search in Google Scholar
 R. Ameling and H. Giessen, “Cavity plasmonics: Large normal mode splitting of electric and magnetic particle plasmons induced by a photonic microcavity,” Nano. Lett., vol. 10, pp. 4394–4398, 2010, https://doi.org/10.1021/nl1019408. Search in Google Scholar
 R. Ameling and H. Giessen, “Microcavity plasmonics: Strong coupling of photonic cavities and plasmons,” Laser Photon. Rev., vol. 7, no. 2, pp. 141–169, 2013, https://doi.org/10.1002/lpor.201100041. Search in Google Scholar
 S. Alrasheed and E. Di Fabrizio, “Effect of surface plasmon coupling to optical cavity modes on the field enhancement and spectral response of dimer-based sensors,” Sci. Rep., vol. 7, p. 10524, 2017, https://doi.org/10.1038/s41598-017-11140-0. Search in Google Scholar
 Y. Zhu, D. J. Gauthier, S. E. Morin, Q. Wu, H. J. Carmichael, and T. W. Mossberg, “Vacuum Rabi splitting as a feature of lineardispersion theory: Analysis and experimental observations,” Phys. Rev. Lett., vol. 64, no. 21, pp. 2499–2502, 1990, https://doi.org/10.1103/physrevlett.64.2499. Search in Google Scholar
 D. G. Baranov, M. Wersäll, J. Cuadra, T. J. Antosiewicz, and T. Shegai, “Novel nanostructures and materials for strong light–matter interactions,” ACS Photonics, vol. 5, no. 1, pp. 24–42, 2018, https://doi.org/10.1021/acsphotonics.7b00674. Search in Google Scholar
 P. Törmä and W. L. Barnes, “Strong coupling between surface plasmon polaritons and emitters: a review,” Rep. Prog. Phys., vol. 78, no. 1, 2014, Art no. 013901, https://doi.org/10.1088/0034-4885/78/1/013901. Search in Google Scholar
 B. Munkhbat, D. G. Baranov, M. Stührenberg, M. Wersäll, A. Bisht, and T. Shegai, “Self-hybridized exciton-polaritons in multilayers of transition metal dichalcogenides for efficient light absorption,” ACS Photonics, vol. 6, no. 1, pp. 139–147, 2019, https://doi.org/10.1021/acsphotonics.8b01194. Search in Google Scholar
 S. R. K. Rodriguez, Y. T. Chen, T. P. Steinbusch, M. A. Verschuuren, A. F. Koenderink, and J. G. Rivas, “From weak to strong coupling of localized surface plasmons to guided modes in a luminescent slab,” Phys. Rev. B, vol. 90, p. 235406, 2014, https://doi.org/10.1103/physrevb.90.235406. Search in Google Scholar
 S. R. K. Rodriguez and J. G. Rivas, “Surface lattice resonances strongly coupled to rhodamine 6G excitons: tuning the plasmonexciton- polariton mass and composition,” Opt. Express, vol. 21, pp. 27411–27421, 2013, https://doi.org/10.1364/oe.21.027411. Search in Google Scholar
 S. Wang, “Strong Light-molecule Coupling: Routes to New Hybrid Materials,” PhD thesis, Strasbourg, France: Université de Strasbourg, 2015. Search in Google Scholar
 S. Wang, S. Li, T. Chervy, A. Shalabney, S. Azzini, E. Orgiu, et al., “Coherent coupling of WS2 monolayers with metallic photonic nanostructures at room temperature,” Nano. Lett., vol. 16, no. 7, pp. 4368–4374, 2016, https://doi.org/10.1021/acs.nanolett.6b01475. Search in Google Scholar
 A. Bisht, J. Cuadra, M. Wersäll, A. Canales, T. J. Antosiewicz, and T. Shegai, “Collective strong light-matter coupling in hierarchical microcavity-plasmon-exciton systems,” Nano. Lett., vol. 19, no. 1, pp. 189–196, 2019, https://doi.org/10.1021/acs.nanolett.8b03639. Search in Google Scholar
 M. Ramezani, A. Halpin, A. I. Fernández-Domínguez, J. Feist, S. R. K. Rodriguez, F. J. Garcia-Vidal, et al., “Plasmon-exciton polariton lasing,” Optica, vol. 4, pp. 31–37, 2017. https://doi.org/10.1364/optica.4.000031. Search in Google Scholar
 P. Yeh, A. Yariv, and C. S. Hong, “Electromagnetic propagation in periodic stratified media. I. General theory*,” J. Opt. Soc. Am., vol. 67, no. 4, pp. 423–438, 1977, https://doi.org/10.1364/josa.67.000423. Search in Google Scholar
 A. Yariv and P. Yeh, “Electromagnetic propagation in periodic stratified media. II. Birefringence, phase matching, and x-ray lasers,” J. Opt. Soc. Am., vol. 67, pp. 438–447, 1977, https://doi.org/10.1364/josa.67.000438. Search in Google Scholar
 C. Menzel, J. Sperrhake, and T. Pertsch, “Efficient treatment of stacked metasurfaces for optimizing and enhancing the range of accessible optical functionalities,” Phys. Rev. A, vol. 93, 2016, Art no. 063832, https://doi.org/10.1103/physreva.93.063832. Search in Google Scholar
 J. Sperrhake, M. Decker, M. Falkner, S. Fasold, T. Kaiser, I. Staude, and T. Pertsch, “Analyzing the polarization response of a chiral metasurface stack by semi-analytic modeling,” Opt. Express, vol. 27, no. 2, pp. 1236–1248, 2019, https://doi.org/10.1364/oe.27.001236. Search in Google Scholar
 D. Baranov, B. Munkhbat, N. Länk, R. Verre, M. Käll, and T. Shegai, “Circular dichroism mode splitting and bounds to its enhancement with cavity-plasmon-polaritons,” Nanophotonics, vol. 9, pp. 283–293, 2019, https://doi.org/10.1515/nanoph-2019-0372. Search in Google Scholar
 F. J. García de Abajo, “Colloquium: Light scattering by particle and hole arrays,” Rev. Mod. Phys., vol. 79, pp. 1267–1290, 2007, https://doi.org/10.1103/revmodphys.79.1267. Search in Google Scholar
 P. Lunnemann, I. Sersic, and A. F. Koenderink, “Optical properties of two-dimensional magnetoelectric point scattering lattices,” Phys. Rev. B, vol. 88, no. 24, p. 245109, 2013, https://doi.org/10.1103/physrevb.88.245109. Search in Google Scholar
 S. Thongrattanasiri, F. H. L. Koppens, and F. J. García de Abajo, “Complete optical absorption in periodically patterned graphene,” Phys. Rev. Lett., vol. 108, 2012, Art no. 047401, https://doi.org/10.1103/physrevlett.108.047401. Search in Google Scholar
 P. de Vries, D. V. van Coevorden, and A. Lagendijk, “Point scatterers for classical waves,” Rev. Mod. Phys., vol. 70, pp. 447–466, 1998, https://doi.org/10.1103/revmodphys.70.447. Search in Google Scholar
 A. Krasnok, D. Baranov, H. Li, M. A. Miri, F. Monticone, and A. Alú, “Anomalies in light scattering,” Adv. Opt. Photon., vol. 11, no. 4, pp. 892–951, 2019, https://doi.org/10.1364/aop.11.000892. Search in Google Scholar
 S. Richter, H. G. Zirnstein, J. Zúñiga Pérez, E. Krüger, C. Deparis, L. Trefflich, et al., “Voigt exceptional points in an anisotropic ZnO-based planar microcavity: square-root topology, polarization vortices, and circularity,” Phys. Rev. Lett., vol. 123, 2019, Art no. 227401, https://doi.org/10.1103/physrevlett.123.227401. Search in Google Scholar
 H. M. Doeleman, F. Monticone, W. den Hollander, A. Alù, and A. F. Koenderink, “Experimental observation of a polarization vortex at an optical bound state in the continuum,” Nat. Photon., vol. 12, pp. 397–401, 2018, https://doi.org/10.1038/s41566-018-0177-5. Search in Google Scholar
 H. J. Carmichael, R. J. Brecha, M. G. Raizen, H. J. Kimble, and P. R. Rice, “Subnatural linewidth averaging for coupled atomic and cavity-mode oscillators,” Phys. Rev. A, vol. 40, pp. 5516–5519, 1989, https://doi.org/10.1103/physreva.40.5516. Search in Google Scholar
 J. M. Raimond and S. Haroche, “Atoms in cavities,” in Confined Electrons and Photons: New Physics and Applications, E. Burstein and C. Weisbuch, Eds., NATO Science Series B, Springer, 1995, pp. 383–426. Search in Google Scholar
 C. Weisbuch, M. Nishioka, A. Ishikawa, and Y. Arakawa, “Observation of the coupled exciton-photon mode splitting in a semiconductor quantum microcavity,” Phys. Rev. Lett., vol. 69, pp. 3314–3317, 1992, https://doi.org/10.1103/physrevlett.69.3314. Search in Google Scholar
 D. G. Lidzey, D. C. Bradley, M. S. Skolnick, T. Virgili, S. Walker, and D. M. Whittaker, “Strong exciton–photon coupling in an organic semiconductor microcavity,” Nature, vol. 395, pp. 53–55, 1998, https://doi.org/10.1038/25692. Search in Google Scholar
 H. S. Wei, C. C. Jaing, Y. T. Chen, C. C. Lin, C. W. Cheng, C. H. Chan, et al., “Adjustable exciton-photon coupling with giant Rabisplitting using layer-by-layer J-aggregate thin films in all-metal mirror microcavities,” Opt. Express, vol. 21, no. 18, pp. 21365–21373, 2013, https://doi.org/10.1364/oe.21.021365. Search in Google Scholar
 M. R. Foreman and F. Vollmer, “Theory of resonance shifts of whispering gallery modes by arbitrary plasmonic nanoparticles,” New. J. Phys., vol. 15, no. 8, 2013, Art no. 083006, https://doi.org/10.1088/1367-2630/15/8/083006. Search in Google Scholar
 J. Xavier, S. Vincent, F. Meder, and F. Vollmer, “Advances in optoplasmonic sensors – combining optical nano/microcavities and photonic crystals with plasmonic nanostructures and nanoparticles,” Nanophotonics, vol. 7, pp. 1–38, 2018. https://doi.org/10.1515/nanoph-2017-0064. Search in Google Scholar
 N. Thakkar, M. T. Rea, K. C. Smith, K. D. Heylman, S. C. Quillin, K. A. Knapper, et al., “Sculpting fano resonances to control photonic–plasmonic hybridization,” Nano. Lett., vol. 17, no. 11, pp. 6927–6934, 2017, https://doi.org/10.1021/acs.nanolett.7b03332. Search in Google Scholar
 S. Rodriguez, “Classical and quantum distinctions between weak and strong coupling,” Eur. J. Phys., vol. 37, 2016, Art no. 025802, https://doi.org/10.1088/0143-0807/37/2/025802. Search in Google Scholar
 A. M. Jayich, J. C. Sankey, B. M. Zwickl, C. Yang, J. D. Thompson, S. M. Girvin, et al., “Dispersive optomechanics: a membrane inside a cavity,” New. J. Phys., vol. 10, 2008, Art no. 095008, https://doi.org/10.1088/1367-2630/10/9/095008. Search in Google Scholar
 F. Bernal Arango, T. Coenen, and A. F. Koenderink, “Underpinning hybridization intuition for complex antennas by magnetoelectric quadrupolar polarizability retrieval,” ACS Photonics, vol. 1, no. 5, pp. 444–453, 2014, https://doi.org/10.1021/ph5000133. Search in Google Scholar
 S. Rodriguez, “Enhancing the Speed and sensitivity of a nonlinear optical sensor with noise,” Phys. Rev. Appl., vol. 13, no. 2, pp. 1–11, 2020, Art no. 024032, https://doi.org/10.1103/physrevapplied.13.024032. Search in Google Scholar
 D. Ballarini and S. D. Liberato, “Polaritonics: From microcavities to sub-wavelength confinement,” Nanophotonics, vol. 8, pp. 641–654, 2019, https://doi.org/10.1515/nanoph-2018-0188. Search in Google Scholar
 I. Sersic, C. Tuambilangana, T. Kampfrath, and A. F. Koenderink, “Magnetoelectric point scattering theory for metamaterial scatterers,” Phys. Rev. B, vol. 83, no. 24, Art no. 245102, 2011, https://doi.org/10.1103/physrevb.83.245102. Search in Google Scholar
 Y. T. Chen, Y. Zhang, and A. F. Koenderink, “General point dipole theory for periodic metasurfaces: magnetoelectric scattering lattices coupled to planar photonic structures,” Opt. Express, vol. 25, pp. 21358–21378, 2017, https://doi.org/10.1364/oe.25.021358. Search in Google Scholar
© 2020 Annemarie Berkhout et al., published by De Gruyter, Berlin/Boston
This work is licensed under the Creative Commons Attribution 4.0 International License.