© The Author(s), published by EDP Sciences, 2024

1 Introduction

Surface Electromagnetic Waves (SEW) represent an interesting option for controlling optical signals on miniaturized chips for integrated optics and sensing applications. Surface Plasmon Polaritons (SPP) are probably the most widely known SEWs, but they exhibit inherent issues related to the ohmic losses introduced by the metallic materials involved. As an alternative, SEWs sustained by dielectric multilayers (ML) have attracted a growing interest in the past decade. This kind of SEW [1] is also referred to as Bloch Surface Waves (BSWs) [2] to highlight the role of the underlying periodic multilayer structure required for their existence. BSWs offer several advantages as compared to SPPs, such as a wide spectral tunability and low losses thanks to the large choice of transparent dielectric materials available for multilayer manufacturing. In addition, BSWs can be either TE- or TM-polarized [3, 4], depending on the multilayer design. Their excitation by pulsed fields has also been recently studied numerically [5].

Being surface waves, BSWs are evanescent in the medium above the multilayer surface. The coupling with free-space radiation in a BSW-based device is therefore critical as it must provide momentum matching beyond the light-line. In most of the applications proposed so far, BSW coupling is performed by means of bulky prisms, either in Kretschmann or Otto configuration [6]. However, more sophisticated approaches have been recently implemented, involving, for example, the use of individual scatterers [7–9] or miniature prisms [10] placed onto the multilayer surface. Another promising option is represented by integrating diffraction gratings within the BSW-supporting structure [11]. This has been done mainly in two different ways: with the grating being fabricated on top of the multilayer [12–14] or buried beneath the multilayer [15]. In the first case, the multilayer is substantially planar, with the exception of the top layer, where the grating unavoidably perturbs the dispersion of the BSW mode (dielectric loading/unloading effect). In the second case, the grating is fabricated on the substrate surface prior to the multilayer deposition, which occurs on the same side. The multilayer itself results to not be perfectly planar because it (partially) conforms to the underlying corrugation. In both configurations, the BSW dispersion is altered by the presence of the grating, which may lead to some difficulties regarding the precise control of optical functions of complex, possibly resonating, BSW-based architectures.

We propose an alternative approach on diffractive coupling for BSWs, with gratings fabricated on the bottom surface of a transparent substrate having the multilayer deposited on the top surface. In particular, we explore the possibility of using two-dimensional gratings to simultaneously couple BSWs propagating in more than two directions by exploiting the momentum distribution of several diffraction orders. Once the mode dispersion of the multilayer is known, our approach facilitates BSW coupling in a controllable way, as far as wavelengths/numbers and propagation directions are concerned. The directional coupling of BSWs has been already tackled in a few previous articles [16–18], although never considered for multiple directions at once. When the optical path through the substrate is also taken into account, our approach allows a predictable control onto the coupling locations of BSWs launched in different directions.

The present paper is composed as follows. We begin, in Section 2, by introducing the grating-based BSW excitation principle and the assumed geometrical configuration. The theoretical framework for grating design, for which we use a rigorous technique known as the Fourier Modal Method (FMM) [19], is described in Section 3. The design process is analogous with the synthesis of grating-based multiple free-space beam splitters [20], but here we need to account for the BSW excitation conditions and the polarization state of the input wave. In Section 4, we first consider BSW stack design, providing a “benchmark” stack employed in the rest of the work, and then cover the design of linear gratings for simultaneous excitation of two counter-propagating BSWs. Such designs are extended in Section 5 to two-dimensional periodic gratings for excitation of either four or six BSWs propagating at 90° or 60° intervals along the stack, respectively. After a discussion presented in Section 6, conclusions are drawn in Section 7.

2 Excitation principle and geometry

Figure 1 illustrates the geometry for the simplest case of excitation of two counter-propagating BSWs. A flat fused silica substrate with refractive index n_sub = 1.462, such as a 0.5 mm or 3 mm-thick SiO₂ plate, is illuminated by a normally incident monochromatic beam (wavelength λ₀) from the medium underneath (air). A linear grating, with period d of the order of λ₀, provided on the air-substrate surface, splits the beam into three transmitted orders propagating within the substrate: the zeroth order m = 0 and the first diffracted orders m = ±1. The orders m = ±1 propagate in directions θ₊₁ and θ₋₁ given by sinθ_±1 = ±λ₀/n_subd towards the multilayer stack on the top surface of the substrate. If |θ_±1| matches the Kretschmann-incidence BSW excitation angle θ_BSW for the given wavelength and polarization state (TE or TM), two counter-propagating BSWs are generated simultaneously. The excitation is efficient as long as the angular spectrum of each diffracted order, which defines the beam divergence, falls essentially within the (stack-dependent) BSW momentum bandwidth. The polarization state of illumination affects the coupling significantly; we will consider only BSW excitation in TE polarization, which generally requires a smaller number of stack layers than TM-polarized BSW excitation.

Figure 1 Principle of MBSW generation: the two-beam case. A binary linear surface-relief grating defined by period d, ridge width c, and ridge height h on the bottom of a substrate of thickness H splits the input beam into two diffracted orders m = −1 and m = +1, which excite BSWs on the top surface of the substrate by interaction with the multilayer stack (ML). We assumed nsub = 1.462 and nsup is air.

The parameters of the system are chosen such that the two BSWs shown in Figure 1 are spatially separated under finite-beam illumination. This feature can be useful in BSW-based platforms such as interferometers [21] and integrated components [22]. First-order diffracted beams are partially reflected at the top interface, thus propagating back into the substrate. The reflected beams continue to propagate according to multiple-reflection paths inside the substrate unless they are extracted by means of diffusers or gratings. At each reflection with the ML interface, coupling to BSW occurs. Stated differently, BSWs are launched at different locations on the ML surface each time the beam is incident on the bottom interface of the dielectric stack, thus leading to the appearance of BSW interference effects unless the substrate thickness H is sufficiently large to minimize spatial overlaps.

3 Theoretical framework

Let us consider a rectangularly periodic grating of period d_x × d_y in the cartesian xy coordinate system and assume a plane-wave illumination (at frequency ω) normally incident onto the substrate from air. In view of the grating equations, the wave vectors of the propagating diffraction orders (m, n) in the substrate are$${\mathit{k}}_{\mathrm{mn}}=k_{\mathrm{xm}}\widehat{\mathit{x}}+k_{\mathrm{yn}}\widehat{\mathit{y}}+k_{\mathrm{zmn}}\widehat{\mathit{z}}$$(1)where$$k_{\mathrm{xm}}=m{K}_x=2\pi m/d_x$$(2a) $$k_{\mathrm{yn}}=n{K}_y=2\pi n/d_y$$(2b) $$k_{\mathrm{zmn}}=\sqrt{k_0^2{n}_{\mathrm{sub}}^2-k_{\mathrm{xm}}^2-k_{\mathrm{yn}}^2}$$(2c)and k₀ = ω/c₀ = 2π/λ₀ is the wave number in vacuum. After defining the radial spatial frequency of the generic order (m, n) as$$k_{\rho \mathrm{mn}}=\sqrt{k_{\mathrm{xm}}^2+k_{\mathrm{yn}}^2}=\sqrt{(m{K}_x{)}^2+(n{K}_y{)}^2}$$(3)the condition k_ρmn < k₀n_sub identifies those diffraction orders propagating within the substrate, the others being evanescent. If we denote the refractive index of the superstrate by n_sup and assume n_sup < n_sub, order (m, n) is evanescent in the superstrate when k_ρmn > k₀n_sup. Considering BSW excitation, we are therefore interested in orders with radial spatial frequencies in the range k₀n_sup < k_ρmn < k₀n_sub. We are primarily interested in the nearest neighbors of the zeroth transmitted order, while higher orders are made evanescent by appropriate choices of d_x and d_y. In the illustrative example presented in Figure 2a, orders (m, n) = (−1, 0) and (m, n) = (+1, 0) fall on the yellow line of radius k_ρBSW, which defines the BSW excitation condition dictated by the ML design.

Figure 2 (a) Diffraction orders of a rectangular lattice in spatial-frequency representation at normal incidence. Diffraction orders are represented by dots at positions kxm = mKx, ky = nKy. Blue and red circles represent the cut-off radial spatial frequencies kρ = k0nsub and kρ = k0nsup, respectively, between which BSW excitation is possible. The yellow circle indicates the radial spatial frequency of BSW for a given ML. (b) Definition of the propagation angles (θmn, ϕmn) of a single transmitted diffracted order (m, n) in its exit plane (the grey rectangle) and the π − σ basis of the diffracted electric field.

Following reference [23], we define the “exit plane” of diffraction order (m, n) as the plane containing the wave vector k _mn and the unit vector $\widehat{\mathit{z}}$. Further, propagation angles θ_mn and ϕ_mn of the transmitted orders, are defined as$$k_{\mathrm{xm}}=k_0{n}_{\mathrm{sub}}\sin {\theta }_{\mathrm{mn}}\cos {\varphi }_{\mathrm{mn}}$$(4a) $$k_{\mathrm{yn}}=k_0{n}_{\mathrm{sub}}\sin {\theta }_{\mathrm{mn}}\sin {\varphi }_{\mathrm{mn}}$$(4b) $$k_{\mathrm{zmn}}=k_0{n}_{\mathrm{sub}}\cos {\theta }_{\mathrm{mn}}$$(4c)as illustrated in Figure 2b. Here ϕ_mn is the azimuthal angle in the range [0, 2π), measured counter-clockwise from the k_x axis, and θ_mn in the range [0, π/2) is the propagation angle measured from the k_z axis. It will prove convenient to use the so-called π–σ basis (or local TM/TE basis) to define the polarization states of the transmitted orders. As described in reference [23], this basis allows us to treat incident fields with any polarization state, including partial polarization. Here, however, we are mainly interested in either fully polarized or unpolarized illumination.

If the incident plane wave is fully polarized, we can use any suitable rigorous grating analysis method (in our case FMM) to determine the transverse Cartesian components e_xmn and e_ymn of the polarization vector for any transmitted order, as discussed shortly below. The longitudinal component of e _mn is fixed by Maxwell’s divergence equation, which gives k _mn ∙  e _mn = 0 and$$e_{\mathrm{zmn}}=-\frac{1}{k_{\mathrm{zmn}}}\left(k_{\mathrm{xm}}{e}_{\mathrm{xmn}}+k_{\mathrm{yn}}{e}_{\mathrm{ymn}}\right).$$(5)

In the π–σ basis the polarization state of any order is described by a two-dimensional vector e _πσmn = [e_πmn, e_σmn]^T, where the π and σ components are explicitly given by$$e_{\pi \mathrm{mn}}=e_{\mathrm{xmn}}\cos {\theta }_{\mathrm{mn}}\cos {\varphi }_{\mathrm{mn}}+e_{\mathrm{ymn}}\cos {\theta }_{\mathrm{mn}}\sin {\varphi }_{\mathrm{mn}}-e_{\mathrm{zmn}}\sin {\theta }_{\mathrm{mn}}$$(6a) $$e_{\sigma \mathrm{mn}}=-e_{\mathrm{xmn}}\sin {\varphi }_{\mathrm{mn}}+e_{\mathrm{ymn}}\cos {\varphi }_{\mathrm{mn}}.$$(6b)

As shown in Figure 2b, the component e_πmn lies in the exit plane, whereas e_σmn is perpendicular to it. Hence, they represent the TM and TE components of the electric field in the exit plane, respectively.

In diffraction by two-dimensionally periodic gratings, the polarization states of the transmitted (and reflected) diffracted orders generally depend on the state of input polarization. We represent the polarization vector of a (generally, elliptically polarized) unit-amplitude input plane wave as$$\mathit{e}=e_x\widehat{\mathit{x}}+e_y\widehat{\mathit{y}}=\widehat{\mathit{x}}\cos \alpha +\widehat{\mathit{y}}\sin \alpha \exp \left(\mathrm{i\delta }\right)$$(7)normalized such that e  = 1. The effect of the grating on transmitted radiation can be analyzed by calculating (by FMM) the transmission coefficients$$T_{\mathrm{xmn}}^{\left(x\right)}, T_{\mathrm{ymn}}^{\left(x\right)},T_{\mathrm{xmn}}^{\left(y\right)},T_{\mathrm{ymn}}^{\left(y\right)}$$(8)

for all diffraction orders, where the superscripts (x) and (y) refer to illumination by a purely x-polarized ($\mathit{e}=\widehat{\mathit{x}}$) or y-polarized ($\mathit{e}=\widehat{\mathit{y}}$) incident wave. The coefficients in equation (8) are precisely the complex vector amplitudes that appear in the Rayleigh plane-wave expansion of the field at the output plane of the grating; see, e.g., equation (5) in reference [20]. For an arbitrarily (fully) polarized incident wave the transverse electric-field components of the transmitted orders are [23]$$e_{\mathrm{xmn}}=T_{\mathrm{xmn}}^{\left(x\right)}{e}_x+T_{\mathrm{xmn}}^{\left(y\right)}{e}_y$$(9a) $$e_{\mathrm{ymn}}=T_{\mathrm{ymn}}^{\left(x\right)}{e}_x+T_{\mathrm{ymn}}^{\left(y\right)}{e}_y.$$(9b)

The longitudinal components e_zmn are obtained from equation (5), and the π–σ representation of each order is given by equations (6). Since the input polarization state affects both the π and σ components, it can be used as a design degree of freedom in multiple-BSW excitation, in addition to the geometrical grating parameters.

It is customary to describe the state of polarization of a fully polarized field by a 2 × 2 polarization matrix J = e ^* e ^T ([24], Sec. 6.3.2). Explicitly, for the incident field,$$\mathbf{J}=\left[\begin{array}{ll}{J}_{\mathrm{xx}}& J_{\mathrm{xy}}\\ J_{\mathrm{yx}}& J_{\mathrm{yy}}\end{array}\right]=\left[\begin{array}{ll}{\left|e_x\right|}^2& e_x^{\mathrm{*}}{e}_y\\ e_y^{\mathrm{*}}{e}_x& {\left|e_y\right|}^2\end{array}\right]$$(10)where the asterisk denotes complex conjugation. Correspondingly, the polarization state of any transmitted order in the π–σ basis is described by ${\mathbf{J}}_{\pi \sigma \mathrm{mn}}={\mathbf{e}}_{\pi \mathrm{mn}}^{*}{\mathbf{e}}_{\sigma \mathrm{mn}}^{\mathrm{T}}$ [23], explicitly$${\mathbf{J}}_{\pi \sigma \mathrm{mn}}=\left[\begin{array}{ll}{J}_{\pi \pi \mathrm{mn}}& J_{\pi \sigma \mathrm{mn}}\\ J_{\sigma \pi \mathrm{mn}}& J_{\sigma \sigma \mathrm{mn}}\end{array}\right]=\left[\begin{array}{ll}{\left|e_{\pi \mathrm{mn}}\right|}^2& e_{\pi \mathrm{mn}}^{*}{e}_{\sigma \mathrm{mn}}\\ e_{\sigma \mathrm{mn}}^{*}{e}_{\pi \mathrm{mn}}& {\left|e_{\sigma \mathrm{mn}}\right|}^2\end{array}\right].$$(11)

The polarization states of the diffracted orders can also be characterized by the Stokes parameters [23]$$S_{0\mathrm{mn}}=J_{\pi \pi \mathrm{mn}}+J_{\sigma \sigma \mathrm{mn}}$$(12a) $$S_{1\mathrm{mn}}=J_{\pi \pi \mathrm{mn}}-J_{\sigma \sigma \mathrm{mn}}$$(12b) $$S_{2\mathrm{mn}}=2\mathfrak{R}\left(J_{\pi \sigma \mathrm{mn}}\right)$$(12c) $$S_{3\mathrm{mn}}=2\mathfrak{I}\left(J_{\pi \sigma \mathrm{mn}}\right)$$(12d)where $\mathfrak{R}$ and $\mathfrak{I}$ denote the real and imaginary parts. The normalized forms of the Stokes parameters are defined as s_jmn = S_jmn/S_0mn (j = 1, 2, 3), and the degree of polarization associated with order (m, n) is given by$$P_{\mathrm{mn}}=\sqrt{s_{1\mathrm{mn}}^2+s_{2\mathrm{mn}}^2+s_{3\mathrm{mn}}^2}.$$(13)

For a fully polarized incident wave, P_mn = 1 for all orders, even though the values of the individual Stokes parameters generally depend on order indices.

In addition to fully polarized illumination, we consider the opposite extreme case of unpolarized illumination. The matrix J for partially polarized light is defined as J = 〈e ^* e ^T〉, where the brackets denote ensemble averaging over all polarization realizations. For unpolarized illumination it has a diagonal form ([24], Sec. 6.3.3)$$\mathbf{J}=\frac{1}{2}\left[\begin{array}{ll}1& 0\\ 0& 1\end{array}\right]$$(14)and the degree of input polarization is P = 0. The polarization matrix associated with order (m, n) can be represented as an average$${\mathbf{J}}_{\pi \sigma \mathrm{mn}}=\frac{1}{2}\left[{\mathbf{J}}_{\pi \sigma \mathrm{mn}}^{\left(x\right)}+{\mathbf{J}}_{\pi \sigma \mathrm{mn}}^{\left(y\right)}\right]$$(15)which remains diagonal because e_x and e_y are uncorrelated. However, since the grating treats these components differently, in general J_ππmn ≠ J_σσmn, implying that the individual orders become partially polarized with P_mn > 0.

In standard beam splitting problems in resonance-domain diffractive optics [20] one is interested in the distribution of the diffraction efficiencies of the propagating orders. Since we have normalized the intensity of the incident field such that S₀ = 1, the diffraction efficiencies are defined as [23]$${\eta }_{\mathrm{mn}}=n_{\mathrm{sub}}\cos {\theta }_{\mathrm{mn}}{S}_{0\mathrm{mn}}.$$(16)

In BSW excitation problem, the design goal is to maximize and equalize the coupling of the incident field to a set of BSW modes with the angle θ_mn equal to θ_BSW. If θ_BSW is the excitation angle for TE polarization, the component e_σmn excites a BSW while e_πmn is non-resonant, and vice versa.

We choose the geometry such that several diffraction orders have the same radial spatial frequency k_ρmn = k₀n_subsinθ_BSW, and therefore lie on the yellow circle depicted in Figure 2a. The relative amplitudes of the excited BSWs are determined by the σ-polarized components J_σσmn in TE polarization and J_ππmn in TM polarization. The fraction$${\kappa }_{\mathrm{mn}}=J_{\sigma \sigma \mathrm{mn}}/J_{\pi \pi \mathrm{mn}}$$(17)provides the ratio of the coupled and uncoupled parts of the incident wave in BSW excitation.

4 Plane-wave design with linear gratings

As evident from the preceding discussion, the π–σ representation returns the BSW excitation problem to the basic TE or TM polarized problem. In addition, since we assume a substrate thickness H ≫ λ₀, the evanescent parts of the diffracted fields above the grating and the BSW field below the stack are spatially well separated. Hence, we may treat the BSW stack design and the grating design as two separate problems. In order to obtain an illustrative stack design useful for our purposes, we fix λ₀ = 514 nm and provide a stack geometry sustaining BSWs at angles between the blue and red lines in Figure 2a. The resulting stack can then be used to design gratings for excitation of BSWs that lie on the yellow circle in Figure 2a.

4.1 Multilayer stack design

Figure 3a shows the assumed stack structure, which consists of N high/low (H/L) refractive index bilayers and a terminating top (T) layer with refractive indices n_H, n_L, n_T and thicknesses h_H, h_L, h_T, respectively. To reduce the number of variable parameters, we consider TE polarization, fix the number of bilayers to N = 6, use refractive indices n_H = 2.520 (TiO₂), n_L = 1.476 (SiO₂), n_T = n_H = 2.520. The thicknesses h_H, h_L, h_T are used to design the stack such that the BSW resonance occurs at an angle θ_mn in the exit plane.

Figure 3 (a) Definition of the multilayer structure and notation. (b) Stack parameters as a function of the ratio kx/k0 = neff for TE-mode BSW excitation with N = 6 bilayers: hH/λ0 (blue), hL/λ0 (red), and hT/λ0 (black). The dots mark the position kx/k0 = 1.1209 for BSW excitation at 50° angle of incidence.

Figure 3b shows the design results. The horizontal axis is k_x/k₀ = n_sub sinθ_mn = n_eff, where n_eff can be interpreted as the effective index of the stack. The plotting range starts from the critical angle of BSW generation and extends to k_x/k₀n_sub, i.e., it spans the region between the blue and red circles in Figure 2a. As n_eff increases, the BSW becomes increasingly buried within the multilayer and acts less like a surface mode. At the same time, all layer thicknesses show a monotonically increasing trend.

4.2 Two-way splitting

As illustrated in Figure 1, the coupling of two counter-propagating BSWs is possible with linear gratings (d_x = d, d_y = ∞). The exit plane of both orders, (m, n) = (−1, 0) and (+1, 0), is the xz plane and the π–σ representation reduces to the standard TM/TE decomposition. Since, by symmetry, η_−1,0 = η_+1,0 for binary profiles defined in the inset of Figure 1, we need to maximize η_+1,0. This also leads to the optimum value of J_σσmn, while J_ππmn = 0. Now we only need to find the values of the fill factor f = c/λ₀ and grating height h/λ₀ that maximize η_+1,0 (= η_−1,0) to also maximize J_σσ.

The grating-design results are summarized in Figure 4. The optimum fill factor remains fairly constant over the entire angular range considered here, whereas the optimum grating height decreases with increasing angle. The efficiencies of all propagating orders are plotted in Figure 4b. At around k_x/k₀ = 1.1209 (corresponding to an excitation angle 50°) we get η_±1,0 ≈ 0.4973. Some light is “lost” in zeroth reflected and transmitted orders when we move close to the cut-off at k_x/k₀ = 1 or towards larger values of k_x/k₀, but the designs remain acceptable over a relatively wide range of excitation angles.

Figure 4 Design of two-way beam splitters. (a) Optimum values of the fill factor f (red) and the relief depth h/λ0 (blue) as a function of the exit angle of the first diffracted order in TE polarization. (b) The corresponding first-order diffraction efficiency η+1,0 = η−1,0 (black), the efficiency η0,0 of the zeroth transmitted order (red), and that of the zeroth reflected order (blue), which is 1 − 2η+1,0 − η0,0 due to energy conservation. The inset shows the grating structure and direction of illumination.

The results in Figures 3b and 4 allow us to design two-way beam splitters for any BSW resonance angle of interest. The stack design for the desired angle is obtained from Figure 3b and the corresponding grating design from Figure 4a. The performance of the design can be evaluated from Figure 4b. To limit the number of variables further, we set θ_BSW = 50°, corresponding to k_x/k₀ = 1.1209. The stack design is marked by the dots in Figure 3b, the optimum parameters for TE excitation with N = 6 bilayers being h_H = 60 nm, h_L = 85 nm, and h_T = 20 nm. Correspondingly, the vertical lines in Figure 4a give a grating design f = 0.2536, h/λ₀ = 0.2752, with η_±1,0 = 0.4973.

Considering the optimized case represented by the dots on Figure 3 and the black dashed line on Figure 4, simulation of the reflected and transmitted coefficients has been performed for the full structure. It implies a grating of period d ≃ 459 nm and fill factor is f = 0.2536 and h/λ₀ = 0.2752 on the lower side of a fused silica wafer on top of which the multilayer is deposited. The multilayer design leads to a Bloch surface wave excited when the first diffracted order emerge from the grating at an angle of 50° (k_x/k₀ = 1.1209) at a wavelength of 514 nm. This is observed in Figure 5a, where a strong dip in reflection arises at this value of k_x/k₀. In Figure 5b the response in wavelength is presented.

Figure 5 Response of the full structure (grating, substrate, multilayer and superstrate). (a) and (b) Reflected (black curves) and transmitted (red curves) first diffracted orders as a function of the normalized wavevector (a) and wavelength (b).

5 Plane-wave design with biperiodic gratings

We proceed to design of two-dimensionally periodic gratings that allow simultaneous excitation of more than two BSWs. Two lattice geometries are considered: square lattices for four-way excitation and hexagonal lattices for six-way excitation.

5.1 Four-way splitting

Let us first consider biperiodic gratings with primitive direct-lattice vectors ${\mathit{a}}_1=d\widehat{\mathit{x}}$, ${\mathit{a}}_2=d\widehat{\mathit{y}}$. The (Wigner–Seitz) primitive cell is square-shaped, covering the area −d/2 < x <d/2, −d/2 < y < d/2. The spatial frequencies of the diffraction orders are then k_xm = mK, k_yn = nK, the coordination number is 4, and the nearest neighbors of the zeroth order (0, 0), namely (m, n) = (+1, 0), (0, +1), (−1, 0), (0, −1), propagate in directions ϕ_+1,0 = 0, ϕ_0,+1 = π/2, ϕ_−1,0 = π, ϕ_0,−1 = 3π/2, respectively. By an appropriate choice of d, all of these four orders can be placed simultaneously on the yellow ring in Figure 2a, thus enabling four-way BSW excitation.

In the design, we found it sufficient to consider binary (z-invariant) relative-permittivity profiles of the particular form$${ϵ}_{\mathrm{r}}\left(x,y,z\right)=\{\begin{array}{cc}{n}_1^2& \mathrm{when} x^2+y^2<r^2\\ n_2^2& \mathrm{otherwise}\end{array}$$(18)in 0 < z < h within the primitive cell. The circular feature defined by the radius r can be either a pillar (n₁ = n_sub, n₂ = 1 or a hole (n₁ = 1, n₂ = n_sub) etched in the substrate. This type of pillar/hole structures can be patterned at a nanometer-scale addressing resolution using electron beam lithography system available to us, and require only a single etching step.

The radius r and the relief depth h can be used as the structural design parameters. Some symmetry rules exist, which are helpful in the design. Since the unit cell and the structure are centered at the origin, the transmission coefficients in equation (8) satisfy the inversion symmetry rules$$T_{x,-m,-n}=T_{\mathrm{xmn}}, T_{y,-m,-n}=T_{\mathrm{ymn}}$$(19)for both (x) and (y) input polarizations. These rules hold regardless of the input polarization state, which however has an effect on the actual values of T_xmn and T_ymn. They reduce the number of orders that we need to (or can) control from four to two: we see from equation (19) that η_−m,−n = η_mn. Similar symmetry rules hold also for J_σσmn and J_ππmn.

We begin the design of four-way couplers by optimizing the structural parameters r and h to minimize the sum of the efficiencies of the reflected and transmitted zeroth orders. This maximizes the combined efficiency of the four nearest-neighbor diffraction orders, and leaves the polarization state of the incident field free for design. Choosing θ_BSW = 50° (d ≈ 0.892λ₀), for either 45° or circularly polarized illumination and considering pillars, we get a design r ≈ 0.201λ₀, h ≈ 0.53176λ₀, which gives reflected and transmitted zero-order efficiencies of ~3.5% and ~5.2%, respectively, leaving the rest of the incident energy to be distributed among the nearest-neighbor orders.

Our remaining target is to equalize (and maximize) the coupling strengths J_σσmn of the four signal orders by designing the input polarization state defined in equation (7). The symmetry in the 4-way splitting implies that there is no structurally induced polarization conversion: for (x)-polarized input we get J_σσmn = 0 for orders (m, n) = (±1, 0), while (y)-polarized input gives J_σσmn = 0 for orders (m, n) = (0, ±1). Considering linearly polarized light, the values of J_σσmn (and J_ππmn) vary rapidly with the angle α. Choosing α ≈ π/4 gives values J_σσmn ≈ 0.063 and κ_mn ≈ 0.359 for all four orders. The same result is obtained also for circularly polarized illumination with α ≈ π/4, δ = ±π/2. Both the optimized diffracted efficiencies and the maximized coupling strengths occur at the same illumination polarization.

Considering unpolarized illumination, the matrix J_πσmn becomes diagonal and the degree of polarization takes the form$$P_{\mathrm{mn}}=\left|s_{1\mathrm{mn}}\right|=\frac{|J_{\pi \pi \mathrm{mn}}-J_{\sigma \sigma \mathrm{mn}}|}{J_{\pi \pi \mathrm{mn}}+J_{\sigma \sigma \mathrm{mn}}}$$(20)

With the present numerical values we obtain P_mn ≈ 0.473 for all nearest-neighbor orders. Even though the excitation wave is partially polarized, we obtain the same values of J_σσmn as above; both of the two mutually uncorrelated components of the incident field contribute to TE-mode BSW excitation.

5.2 Six-way splitting

Let us consider a grating with hexagonal symmetry, which allows simultaneous excitation of six BSWs. The primitive vectors are now ${\mathit{a}}_1=d\widehat{\mathit{x}}$ and ${\mathit{a}}_2=(d/2)\widehat{\mathit{x}}+(\sqrt{3}d/2)\widehat{\mathit{y}}$, and the Wigner–Seitz primitive cell is a hexagon as shown in Figure 6a. It will, however, be convenient for our purposes to define a rectangular direct-lattice cell as in reference [20], which covers the spatial region −d/2 < x <d/2, $-\sqrt{3}d/2<y<\sqrt{3}d/2$ in Figure 6a. This alternative lattice representation simplifies the visualization of the geometry. It also allows the use of FMM in a cartesian instead of a non-orthogonal basis, as in reference [20], though in the present work we actually used the latter basis.

Figure 6 (a) The spatial structure of a hexagonal grating. The hexagon shows the spatial Wigner–Seitz primitive cell, the green features illustrate the grating structure, and the blue rectangle shows the non-primitive cartesian cell. (b) The spatial-frequency grid. The filled and empty dots represent the allowed and forbidden orders of the hexagonal lattice. The blue, yellow, and red circles have the same meaning as in Figure 2a.

The spatial-frequency structure defined by the reciprocal-lattice primitive vectors ${\mathit{b}}_1=K\widehat{\mathit{x}}-(K/\sqrt{3})\widehat{\mathit{y}}$, ${\mathit{b}}_2=(2K/\sqrt{3})\widehat{\mathit{y}}$ with K = 2π/d is illustrated in Figure 6b, where the solid green circles show the locations of the allowed orders in the cartesian (k_x, k_y) system. The empty circles represent the orders of the rectangular spatial lattice, which are forbidden by the hexagonal symmetry. The yellow circle connects the six nearest neighbors of the zeroth order that satisfy the condition for BSW excitation simultaneously: orders (m, n) = (+1, +1), (0, +2), (−1, +1), (−1, −1), (0, −2), (+1, −1) of the rectangular lattice, with exit planes at angles π/6 + qπ/3, q = 0,…,5. The excited BSWs propagate along the surface of the stack in these directions.

In hexagonal lattice geometry, the symmetry rules in equation (19) ensure J_{σσ,−1,−1} = J_σσ,1,1, J_σσ,0,−2 = J_σσ,0,2, J_σσ,−1,1 = J_σσ,1,−1. These symmetries leave us three pairs of orders to control, and we expect to need additional structural freedom compared to the 4-wave case. Let us nevertheless see what designs are possible with circular pillars by following the same strategy as above. An important difference is that in the hexagonal geometry we do have structurally induced polarization conversion.

By optimizing r and h for pillars, we get h ≈ 0.422λ₀ and r ≈ 0.254λ₀, which leaves a combined efficiency of ~0.884 available for the 6 orders of interest. The distribution of J_σσmn again depends on input polarization. We found that it is not possible to equalize the coupling exactly for all six orders, but using circularly polarized light with α = π/4 and δ = 0.486π we have J_σσ,1,1 = 0.097, J_σσ,0,2 = 0.102, and J_σσ,1,−1 = 0.113, respectively. Similarly, for κ_mn, we have κ_1,1 = 1.691, κ_0,2 = 1.903, and κ_1,−1 = 2.310. Though it is not of concern for the present purposes, it is worth noting that the diffraction efficiencies are: η_1,1 ≈ 0.152, η_0,2 ≈ 0.145, and η_1,−1 ≈ 0.144. As with the four-wave case, the design with circular pillars works also for circularly polarized or unpolarized illumination but the exact values of J_σσmn depend on polarization, but are within the same range as above. For unpolarized illumination, the degrees of polarization of the individual orders are nearly equal, P_mn ≈ 0.3226.

In Figure 7, we show the field amplitude distribution associated with the six-way coupling geometry in the xz-plane, i.e., crossing the multilayer, when illuminated with a 45° polarized light wave. The field is evaluated over 3-unit cells, i.e., 3 grating periods, in the x-direction. It shows, as expected, a strong field on the upper medium. Such a structure is ideal for sensing applications, especially when providing multiple sensing areas thanks to the splitting of the BSW excitation.

Figure 7 Field amplitude distribution across the multilayer (xz-plane). The illumination polarization was set to 45°. The dashed lines superimposed on the field represent the multilayer interfaces.

6 Discussion

Throughout the paper we have considered normally incident illumination. The use of non-normal incidence could potentially allow us to consider other combinations of diffracted orders being simultaneously resonant. Changing the angle of incidence moves the grid of diffracted orders transversely in Figure 2a with respect to the circles centered at the origin. For instance, if the propagation direction of the incident field is chosen as (k_xi, k_yi) = (0, k_yi), increasing k_yi moves the grid downwards in k_y direction, giving the orders at positions k_xm = mK_x, k_yn = k_yi +nK_y. Hence the three orders (m, n) = (0, 0) and (m, n) = (±1, 0) would have a common radial spatial frequency if $k_{y\mathrm{i}}=-(K_x^2+K_y^2)/2K_y$, being therefore available for 3-way BSW excitation. To avoid order (0, 2) from occupying the same ring as the zeroth order, we would need to choose K_x ≠ K_y. However, placing the (yellow) BSW resonance ring outside the blue ring in Figure 2a requires k_ρ00 > k₀n_sub, which is not possible with incidence from air. Hence a Kretschmann excitation geometry would be needed, thus sacrificing the compact footprint of the setup.

As an alternative to the geometry considered in this paper, we could consider having the splitter grating and the BSW stack as an integrated structure. This would still allow a compact platform at normal incidence, but the grating design and BSW stack design would not be independent anymore. As a first drawback, the splitter grating would most likely have to be rather thick (~λ₀) to suppress the zeroth transmitted order, preventing the possibility of etching it in the top ML layer. As a consequence, a strong degrading effect in the excited BSWs would be expected. Alternatively, one could use a highly index-modulated splitter grating with a flat top surface immediate below the stack. This would partially alleviate the dependency in the ML and the grating design, but presumably the BSWs would be less affected.

7 Conclusions

In summary, we considered grating design for 2-way, 4-way, and 6-way BSW coupling at normally incident but arbitrarily polarized illumination of gratings with linear, square, and hexagonal symmetries. The plane-wave designs feature ideal TE-mode BSW coupling in the two-wave case. In the other cases the non-resonant parts of the excitation orders cannot be eliminated simultaneously, and they actually dominate the resonant (σ polarized) parts by a factor of ~2.8 in the 4-wave case. The 6-wave case reveals the opposite observation with the resonant part dominating by a factor of ~1.968 making them ideal candidate for TE-mode BSW excitation. Nevertheless, in the 4-wave case all four coupling ratios can be made equal, and in the 6-wave case practically equal, for several input polarization states of practical significance.

Acknowledgments

We want to pay our respects to our mentor, teacher, and friend, Dr Jari Turunen, who passed away in May 2023. He was a dedicated professor in the Center for Photonics Sciences of the Department of Physics and Mathematics at the University of Eastern Finland, Joensuu, who was passionate about research.

Funding

The work was partially funded by the Academy of Finland through project 333938 and the Flagship Programme PREIN (346518). E. Descrovi acknowledges the funding received by Italian “Ministero dell’Universitá e della Ricerca” under the “Dipartimento di Eccellenza 2018-2022” program.

Conflicts of Interest

The authors declare no conflicts of interest.

Data availability statement

All data generated or analyzed during this study are included in this published article.

Author contribution statement

ALA developed the codes and performed the numerical simulations. HP contributed in developing the codes and to the discussion. ED and MR supervised the work. MR and JT proposed the concept. JT developed the theory. All authors contributed to the manuscript preparation, have read and approved the final manuscript.

References

All Figures

Figure 1 Principle of MBSW generation: the two-beam case. A binary linear surface-relief grating defined by period d, ridge width c, and ridge height h on the bottom of a substrate of thickness H splits the input beam into two diffracted orders m = −1 and m = +1, which excite BSWs on the top surface of the substrate by interaction with the multilayer stack (ML). We assumed nsub = 1.462 and nsup is air.

In the text

Figure 2 (a) Diffraction orders of a rectangular lattice in spatial-frequency representation at normal incidence. Diffraction orders are represented by dots at positions kxm = mKx, ky = nKy. Blue and red circles represent the cut-off radial spatial frequencies kρ = k0nsub and kρ = k0nsup, respectively, between which BSW excitation is possible. The yellow circle indicates the radial spatial frequency of BSW for a given ML. (b) Definition of the propagation angles (θmn, ϕmn) of a single transmitted diffracted order (m, n) in its exit plane (the grey rectangle) and the π − σ basis of the diffracted electric field.

In the text

	Figure 3 (a) Definition of the multilayer structure and notation. (b) Stack parameters as a function of the ratio kx/k0 = neff for TE-mode BSW excitation with N = 6 bilayers: hH/λ0 (blue), hL/λ0 (red), and hT/λ0 (black). The dots mark the position kx/k0 = 1.1209 for BSW excitation at 50° angle of incidence.
In the text

Figure 4 Design of two-way beam splitters. (a) Optimum values of the fill factor f (red) and the relief depth h/λ0 (blue) as a function of the exit angle of the first diffracted order in TE polarization. (b) The corresponding first-order diffraction efficiency η+1,0 = η−1,0 (black), the efficiency η0,0 of the zeroth transmitted order (red), and that of the zeroth reflected order (blue), which is 1 − 2η+1,0 − η0,0 due to energy conservation. The inset shows the grating structure and direction of illumination.

In the text

	Figure 5 Response of the full structure (grating, substrate, multilayer and superstrate). (a) and (b) Reflected (black curves) and transmitted (red curves) first diffracted orders as a function of the normalized wavevector (a) and wavelength (b).
In the text

Figure 6 (a) The spatial structure of a hexagonal grating. The hexagon shows the spatial Wigner–Seitz primitive cell, the green features illustrate the grating structure, and the blue rectangle shows the non-primitive cartesian cell. (b) The spatial-frequency grid. The filled and empty dots represent the allowed and forbidden orders of the hexagonal lattice. The blue, yellow, and red circles have the same meaning as in Figure 2a.

In the text

	Figure 7 Field amplitude distribution across the multilayer (xz-plane). The illumination polarization was set to 45°. The dashed lines superimposed on the field represent the multilayer interfaces.
In the text