*Article* **Meta-Atoms with Toroidal Topology for Strongly Resonant Responses**

**Odysseas Tsilipakos 1,\* , Zacharias Viskadourakis <sup>2</sup> , Anna C. Tasolamprou 2,3 , Dimitrios C. Zografopoulos <sup>4</sup> , Maria Kafesaki 2,5 , George Kenanakis 2,\* and Eleftherios N. Economou 2,6**

	- **\*** Correspondence: otsilipakos@eie.gr (O.T.); gkenanak@iesl.forth.gr (G.K.)

**Abstract:** A conductive meta-atom of toroidal topology is studied both theoretically and experimentally, demonstrating a sharp and highly controllable resonant response. Simulations are performed both for a free-space periodic metasurface and a pair of meta-atoms inserted within a rectangular metallic waveguide. A quasi-dark state with controllable radiative coupling is supported, allowing to tune the linewidth (quality factor) and lineshape of the supported resonance via the appropriate geometric parameters. By conducting a rigorous multipole analysis, we find that despite the strong toroidal dipole moment, it is the residual electric dipole moment that dictates the electromagnetic response. Subsequently, the structure is fabricated with 3D printing and coated with silver paste. Importantly, the structure is planar, consists of a single metallization layer and does not require a substrate when neighboring meta-atoms are touching, resulting in a practical, thin and potentially low-loss system. Measurements are performed in the 5 GHz regime with a vector network analyzer and a good agreement with simulations is demonstrated.

**Keywords:** toroidal dipole; metasurfaces; multipole expansion; broken symmetry; 3D printing; microwaves

#### **1. Introduction**

Toroidal multipoles are a class of fundamental electromagnetic excitations that complement the more familiar electric and magnetic multipole families [1–5]. The toroidal dipole, the lowest order member of the toroidal family, first considered by Zel'dovich [6], originates from conduction or displacement currents circulating on a torus along the meridians, producing a closed loop of magnetic field circulation. Although it is distinctly different in its construction from the electric dipole (**p**, a simple separation of electric charges), it emits radiation with the same angular momentum and parity properties as the electric dipole. As a result, the two cannot be easily discerned by observing the far field radiation.

Recently, metamaterials have provided fertile ground for the observation of the toroidal dipole and higher order multipoles through the ability to judiciously shape the meta-atom/meta-molecule geometry in the unit cell. The first demonstration was published in 2010 by Kaelberer et al. [7]. Since then, a broad range of structures have been studied, based on both dielectric materials (displacement currents) [8–12] and metals (conduction currents) [13–15]. In particular, planar and, ideally, single-metallization-layer structures are favorable for easier fabrication [16–20]. An important point of attention is the possibility

**Citation:** Tsilipakos, O.;

Viskadourakis, Z.; Tasolamprou, A.C.; Zografopoulos, D.C.; Kafesaki, M.; Kenanakis, G.; Economou, E.N. Meta-Atoms with Toroidal Topology for Strongly Resonant Responses. *Micromachines* **2023**, *14*, 468. https:// doi.org/10.3390/mi14020468

Academic Editors: Luigi Sirleto and Giancarlo C. Righini

Received: 24 January 2023 Revised: 9 February 2023 Accepted: 15 February 2023 Published: 17 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

to attain a quasi-dark (almost non-radiating) state, which is termed a "dynamic anapole", and results from the near-destructive interference of the toroidal dipole, **T**, with the electric dipole, **p** [21]. As perfect destructive interference is approached through **p** + *ik***T**→0, one can obtain arbitrarily high radiation quality factors, *Q*rad (provided that contribution from other multipoles is suppressed). The total quality factor, *Q*tot, can be very high as well, provided that a material system with low absorption is used [22]. This concept is related to the topics of trapped [23], broken-symmetry [24,25] and BIC (bound states in the continuum) [26] resonances. In such cases, sharply resonant responses and narrow spectral linewidths can be achieved, which are particularly useful for functionalities requiring enhanced local fields. Examples of such applications with metasurfaces supporting anapole states are, for instance, nonlinear effects [27–29], lasers [30] and (bio)sensing [31,32].

In recent years, a broad range of metasurfaces (single-layer metamaterial structures) have been proposed for demonstrating toroidal dipole- and anapole-based phenomena [33]. In most cases, the designs are based on a toroidal topology. Here, we adopt such a design aimed at supporting an anapole state [13]; it consists of a unit cell geometry meant to enhance the toroidal character, it is planar and it requires a single metallization layer that can be free-standing (without a substrate), which are important advantageous traits regarding fabrication and practical applications. We show that, indeed, a quasi-dark resonant state can be supported, which leads to an arbitrarily narrow linewidth, limited only by the resistive quality factor. Varying the appropriate geometric parameters, both the linewidth *and* lineshape of the resonance can be controlled, offering valuable degrees of freedom in shaping the spectral response. By performing a rigorous multipole expansion, we find that the toroidal dipole moment in the structure is strong. However, the corresponding far-field scattering is cancelled exactly by the magnetic quadrupole, **Q***m*. Thus, a quasi-dark resonance is achieved by controlling the residual electric dipole moment through asymmetry in the unit cell geometry (central vs. outer gaps) and not by satisfying the anapole condition. Subsequently, we experimentally verify our calculations through measurements within a rectangular metallic waveguide setup. More specifically, toroidal meta-atom samples, fabricated by a very low-cost technique based on 3D printing and subsequent metallization with a conductive paste, are loaded in the waveguide cross-section, and the S parameters (reflection and transmission coefficients) are measured with a vector network analyzer.

One goal of our work is to demonstrate a practical, easily realizable meta-atom geometry that provides freedom in shaping the electromagnetic response in terms of both linewidth and lineshape. Importantly, this study also clearly demonstrates that an interpretation suggested by the unit cell geometry (its topology) does not necessarily lead directly to the physical interpretation of the phenomenon taking place. Instead, in order to reveal the actual physical behavior, a careful analysis should be performed in all cases. Thus, our work can act as a reference for the mindful discussion of the physical mechanisms behind the resonant response of toroidal metamaterials (but not limited exclusively to this special class of structures). In this paper, the physical evidence is provided by the rigorous multipole expansion and is supported by the experimental verification of the spectral response.

#### **2. Materials and Methods**

#### *2.1. Metasurface Full-Wave Simulations and Multipole Expansion*

The metasurfaces studied in this work are composed of a unit cell made of conductive material. In order to study the effect of finite conductivity and increasing resistive loss, conductivities in the range of 103–10<sup>7</sup> S/m have been assumed. Full-wave simulations are conducted with the finite element method (FEM) using the commercial software COMSOL Multiphysics. They concern scattering and guided-wave simulations with a CW excitation and eigenvalue simulations without excitation.

Two types of metasurfaces are considered. The first is a free-space metasurface, where the excitation is a normally incident plane wave with *E<sup>y</sup>* polarization (Figure 1a). In this case, a single unit cell is simulated with periodic boundary conditions in the *x*- and *y*-axes.

In the second case, two periods of the metasurface along the *x*-axis are inserted within a rectangular metallic waveguide (Figure 1b), since the aspect ratio of the waveguide cross-section (*a* × *b*) is approximately 2:1. The structure is excited with the TE<sup>10</sup> mode of the waveguide, exhibiting the characteristic field distribution *Ey*(*x*) ∼ sin(*πx*/*a*). The waveguide walls are modeled as perfect electric conducting (PEC) walls.

The multipole expansion is performed for the free-space geometry by using the expressions found in ref. [3]. The field scattered by the terms of the multipole expansion can be used for reconstructing the reflected and transmitted field [3,34,35]. The correctness of the calculated multipole moments is checked through the agreement between the direct and reconstructed reflection/transmission. This way, we can also specify where to truncate the multipole expansion without sacrificing reconstruction accuracy.

**Figure 1.** (**a**) Periodic metasurface illuminated with a normally incident plane wave of *Ey* polarization. (**b**) Two unit cells fitted in a rectangular metallic waveguide. (**c**) 3D-printed unit cells using a PLA filament. (**d**) Structure after coating with a conductive silver paste. (**e**) Test fitting of the structure under study on the flange of the waveguide-to-coax adapter. The unit cells are lightweight and can be simply glued on a piece of paper and inserted at any junction between waveguide segments. (**f**) Measurement setup including a vector network analyzer allowing to measure reflection (S11) and transmission (S21) coefficients.

#### *2.2. Meta-Atom Fabrication with 3D Printing*

For a fast and cost-effective fabrication of the metasurface under study, 3D printing technology was employed. A commercial fused filament fabrication (FFF) 3D printer was utilized (MakerBot Replicator 2x, New York, NY, USA). A common polylactic acid (PLA) filament was used as a spool material for building the unit cell geometry (Figure 1c). Details regarding the printing process and conditions can be found in refs. [36,37]. Freestanding meta-atoms, such as those depicted in Figure 1c, were successfully fabricated. The fabricated metasurfaces exhibit negligible electrical conductivity, since PLA is an insulator. Therefore, the meta-atoms were subsequently coated with a thin (∼100 µm) layer of conductive silver paste, as shown in Figure 1d. The silver epoxy exhibits electrical conductivity as high as 104–10<sup>5</sup> S/m. We have deduced these values in earlier works [36,37] by comparing experimental results with simulations. Such a meta-atom configuration, in which an insulating core is coated with an adequately thick metallic paint, has been successfully employed at microwave frequencies [36,38].

#### *2.3. Electromagnetic Characterization with Rectangular Waveguide Setup*

The electromagnetic response of the investigated structure was experimentally confirmed through microwave measurements. To obtain a well-defined and highly reproducible measurement environment, which is not susceptible to external perturbations, we employed a rectangular metallic waveguide setup. More specifically, we fitted 2 × 1 unit cells inside the waveguide to fill the cross-section (aspect ratio 2:1). A test fitting of the two-unit-cell structure within the cross-section of a waveguide-to-coax adapter is shown in Figure 1e.

Measurements were conducted in the vicinity of 5 GHz by using an HP 8722ES vector network analyzer (Agilent Technologies Inc., Santa Clara, CA, USA). WR-187 rectangular waveguides (cross-section: 47.55 mm × 22.2 mm) were used, able to cover the frequency range of 3.95–5.85 GHz. The measurement setup is depicted in Figure 1f and allows for measuring both the reflection (*S*11) and transmission (*S*21) coefficients.

#### **3. Results**

#### *3.1. Free-Space Metasurface with Controllable Strongly Resonant Response*

The meta-atom considered in this work is depicted in Figure 2a. Its geometry is selected in such a way so as to lead to a circulation of induced magnetization giving rise to a toroidal dipole moment [13]; we will come back to this while discussing Figure 2c. The gap in the inner branch is denoted by *g*<sup>1</sup> and the gaps in the outer branches are denoted by *g*2. For the inner and outer radii, *R*inn and *R*out, it holds *R*inn = *R*out − *w*, where *w* is the width of the central and outer branches. The lattice constant (pitch) is denoted by *a* for both *x* and *y* axes (square periodicity). When *R*out > *a*/2, the neighboring meta-atoms are touching; in this case, the conductive meta-atoms can form an interconnected metasurface and no substrate is strictly required. This can be important for avoiding additional resistive loss from the substrate and avoiding excess thickness. Preserving the vertical symmetry of the structure has been also shown to lead to higher modulation depth in transmission [39].

First, we examined a periodic metasurface. The reflection/transmission/absorption power coefficients (R/T/A) for plane wave scattering under normal incidence (*E<sup>y</sup>* polarization) are depicted in Figure 2b. The dimensions of the specific example are *a* = 15 mm, *R*out = 7 mm, *R*inn = 5 mm (*w* = 2 mm), *g*<sup>1</sup> = 0.5 mm and *g*<sup>2</sup> = 0.5 mm. The height of the conductive meta-atom is *h* = 1 mm. In this case, the meta-atoms are not touching. We observed a Fano lineshape associated with the excitation of a resonance at ∼8.71 GHz and its interference with a non-resonant electric polarizability background. Even if the gaps are closed and the meta-atom is non-resonant, the conductive unit cell would still produce significant reflection. The conductivity adopted in these simulations is high (*σ* = 10<sup>7</sup> S/m), leading to low absorption.

We subsequently performed an eigenvalue analysis to identify the nature of the supported resonance. The mode profile is depicted in Figure 2c; the color corresponds to the *E<sup>y</sup>* component (real part) and the arrows to the magnetic field. The magnetic field circles around the central branch and this should produce a toroidal dipole moment along the *y* axis (*Ty*). Looking at the fields in the gaps, one realizes that the central electric dipole moment is opposed by the two contributions of the outer gaps; when the opposing *p<sup>y</sup>* contributions compensate each other exactly, no coupling via the electric dipole moment is possible. This balance can be used to tune the radiative strength of an electrically coupled meta-atom [40].

To further study the physics of the scattering process, we performed a multipole expansion based on the induced conduction currents on the meta-atom using the expressions found in ref. [3], and we focused on the multipole moments that produce scattered fields toward the *z*-axis with *E<sup>y</sup>* (*Hx*) polarization. The eight leading terms of the expansion that satisfy the above criterion are *py*, *mx*, *Ty*, *Q<sup>m</sup> xz*, *Q<sup>e</sup> yz*, *Q<sup>T</sup> yz*, *O<sup>m</sup> xzz* and *O<sup>e</sup> yzz*. The results are depicted in Figure 2d. It can be seen that the rationale behind the meta-atom geometry has indeed resulted in a strong toroidal dipole moment contribution in the spectral neighborhood of the resonance frequency. However, the magnetic quadrupole moment is equally

dominant. In fact, their scattered fields cancel out, since on resonance it holds exactly *E Ty* sca = −*E Q<sup>m</sup> xz* sca . As a result, the response is dictated by the residual electric dipole moment, *py*, which is the third most dominant contribution. This can be further verified in Figure 2e. We first confirm the validity of the multipole expansion by reconstructing the reflection coefficient. The agreement is excellent and shows that using the eight leading terms up to the electric octupole is more than adequate in this case. Importantly, if we use only the field scattered by the electric dipole moment, we obtain a fairly accurate description of the scattering process, especially near the resonance, where there is mutual cancellation of the dominant *T<sup>y</sup>* and *Q<sup>m</sup> xz* terms.

**Figure 2.** (**a**) Meta-atom geometry. When *R*out > *a*/2, the neighboring meta-atoms are touching; in this case, no substrate is necessary. (**b**) R/T/A power coefficients for plane-wave scattering by the periodic metasurface under normal incidence (*E<sup>y</sup>* polarization). The dimensions are *a* = 15 mm, *R*out = 7 mm, *R*inn = 5 mm (*w* = 2 mm), *g*<sup>1</sup> = 0.5 mm and *g*<sup>2</sup> = 0.5 mm. The thickness of the conductive meta-atom is *h* = 1 mm. (**c**) Supported resonance associated with the spectral feature in panel (**b**). The color corresponds to the *E<sup>y</sup>* component (real part) and the arrows correspond to the magnetic field distribution. The eigenmode is characterized by a residual electric dipole moment due to counteracting contributions from the inner gap vs. the outer gaps. The magnetic field circulation gives rise to a strong toroidal dipole moment, *Ty*. (**d**) Power scattered from each multipole. The toroidal dipole and magnetic quadrupole moments are dominant at the resonant frequency. However, their scattered fields cancel out, since *E Ty* sca = −*E Q<sup>m</sup> xz* sca exactly. Thus, the response is dictated by the residual electric dipole moment, *py*. (**e**) Comparison of the reflection coefficient calculated from the full-wave simulation with those reconstructed from the multipole moments. The response can be described fairly accurately when only the electric dipole moment is considered.

Next, we investigated different geometric degrees of freedom for tuning the resonance linewidth and lineshape. In Figure 3a, we varied the central gap (*g*1) while keeping the outer gaps (*g*2) constant at 0.5 mm. This controls the residual electric dipole moment and consequently the radiative strength of the resonance. For *g*<sup>1</sup> = 0.4 mm, the two opposing contributions (see Figure 2c) cancel each other out almost perfectly. In this case, very little coupling via the electric dipole moment is possible and the resonance becomes quasi-dark, resulting in a very narrow linewidth (high quality factor). The non-negligible reflection that remains is due to the non-resonant electric polarizability background. For *g*<sup>1</sup> > 0.4 mm, the central *p<sup>y</sup>* contribution dominates over those originating from the outer gaps, while the opposite occurs for *g*<sup>1</sup> < 0.4 mm. In both cases, the linewidth (and radiative strength) increases. The corresponding total quality factors are compiled in Table 1. They have been obtained from the solution of an eigenvalue problem using the complex eigenfrequency *ω*˜ = *ω*0 + *iω*00 (the eigenvalue) through the expression *Q* = *ω*0/(2*ω*00). For additional details regarding approaches to calculating the quality factor we refer the interested reader to ref. [22]. We document the total quality factor, *Q*tot, which is directly associated with the linewidth of the spectral feature, along with the corresponding resonant frequency. Note that *Q*rad ≈ *Q*tot holds, since *Q*res → ∞ ("res" stands for resistive) for such a high value of the conductivity (*σ* = 10<sup>7</sup> S/m).

**Figure 3.** (**a**) Tuning the central gap (*g*<sup>1</sup> ), while keeping other gaps (*g*2) constant at 0.5 mm. This controls the residual electric dipole moment and consequently the radiative strength of the resonance. (**b**) Varying the outer radius for *R*inn = *R*out − *w* with *w* = 2 mm. As the radius increases, the non-resonant (background) electric dipole moment changes and with it the Fano lineshape. For *R*out values of 7.6, 7.7, and 7.8 mm, the meta-atoms are touching; this modifies the residual electric dipole moment, as the outer gaps of each meta-atom begin to merge with those of the neighboring one.

**Table 1.** Quality factors for the cases depicted in Figure 3a. *g*<sup>2</sup> is held constant at 0.5 mm. They have been calculated using the complex eigenfrequency through *Q* = *ω*0/(2*ω*00).


A different option is explored in Figure 3b, where we tune the outer radius. In all cases, the inner radius is also appropriately modified, so that the branch width is kept constant, i.e., *R*inn = *R*out − *w*, with *w* = 2 mm. This primarily modifies the non-resonant (background) electric dipole moment, leading to changes in the Fano lineshape from highly asymmetric to quasi-Lorentzian. For *R*out values of 7.6, 7.7, and 7.8 mm, the meta-atoms are touching, this additionally modifies the residual electric dipole moment, as the outer gaps of each meta-atom begin to merge with those of the neighboring atom, producing high values of background reflectance. Importantly, by tuning *R*out we are able to achieve very narrow linewidths and, at the same time, a large resonance depth (compare Figure 3b with Figure 3a). The total quality factors for *R*out = 7.7 and 7.8 mm are ∼5000 and ∼30,000, respectively. Note that this quality factor corresponds to radiation losses and will drop if resistive losses increase.

#### *3.2. Meta-Atoms in a Rectangular Waveguide Setup—Experimental Verification*

In this Section, the unit cells under study are inserted inside a rectangular waveguide environment (see Figure 1). This is exercised in order to end up with a well-defined and highly reproducible measurement environment, which is not susceptible to external perturbations. Furthermore, the electromagnetic response is anticipated to be similar to the free-space periodic structure, since the top/bottom PEC walls emulate a periodic repetition along the *y* axis and the main difference is the *sin*(*πx*/*a*) profile of the incident guided wave along the *x* axis [35,37,41].

We first performed simulations. We tuned the dimensions in order to bring the resonant frequency near to the center of the band of the WR-187 rectangular waveguide (3.95–5.85 GHz) and fit 2×1 meta-atoms in the cross-section (*a* × *b* = 47.55 mm × 22.2 mm). The dimensions are *R*out = 12 mm, *R*inn = 8.6 mm, *w* = 3.4 mm, *g*<sup>1</sup> = 0.7 mm and *g*<sup>2</sup> = 1 mm. The thickness of the conductive meta-atom is *h* = 3.4 mm. The response is depicted in Figure 4, where the reflection, *R* = |*S*11| 2 , transmission, *T* = |*S*21| 2 , and absorption, *A* = 1 − |*S*11| <sup>2</sup> − |*S*21<sup>|</sup> 2 , are plotted. The meta-atoms are touching and the resonance lineshape resembles the respective cases in Figure 3b. The meta-atoms will be fabricated via first 3D printing a dielectric (PLA) scaffold and then by coating with a silver paste, which is characterized by a limited conductivity of approximately *σ* = 104–10<sup>5</sup> S/m. Thus, in Figure <sup>4</sup> we examine the conductivities *<sup>σ</sup>* <sup>=</sup> <sup>10</sup><sup>5</sup> S/m, *<sup>σ</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>4</sup> S/m, *<sup>σ</sup>* <sup>=</sup> <sup>10</sup><sup>4</sup> S/m and *<sup>σ</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>3</sup> S/m. As the conductivity decreases, the spectral features become broader and the absorption increases up to a maximum of 0.5 (for an electrically polarizable structure) before starting to decrease again (under-coupling regime). The mode profile is depicted in the inset of Figure 4d. Notice that the gaps near the edges of the waveguide do not accommodate strong fields, since the incident TE<sup>01</sup> waveguide mode is zeroed out at the side walls (cf. Figure 2).

**Figure 4.** Simulations of the 3D-printed conductive unit cells within the WR-187 rectangular waveguide. The reflection (*R* = |*S*11| 2 ), transmission (*T* = |*S*21| 2 ) and absorption (*A* = 1 − |*S*11| <sup>2</sup> − |*S*21<sup>|</sup> 2 ) power coefficients are plotted for different conductivities: (**a**) *<sup>σ</sup>* <sup>=</sup> <sup>10</sup><sup>5</sup> S/m, (**b**) *<sup>σ</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>4</sup> S/m, (**c**) *<sup>σ</sup>* <sup>=</sup> <sup>10</sup><sup>4</sup> S/m and (**d**) *<sup>σ</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>3</sup> S/m. As the conductivity is decreased, the spectral feature becomes broader and the absorption on resonance increases up to a maximum of 0.5 before starting to decrease again.

In Figure 5, we depict the measurements of the 3D-printed conductive unit cells within the WR-187 rectangular waveguide. The reflection (*R* = |*S*11| 2 ) and transmission (*T* = |*S*21| 2 ) power coefficients are plotted in the frequency range of 4.5–5.5 GHz. The lineshape of the spectral feature and the linewidth (full-width half maximum of approximately 100 MHz, as extracted from the spectral feature in reflection) are in good qualitative agreement with those predicted by the simulations and specifically Figure 4c, which corresponds to a conductivity value we anticipate for the silver epoxy. The resonant frequency is slightly different; we attribute this discrepancy to the limited accuracy of the geometric dimensions achieved in the actual fabricated sample.

**Figure 5.** Measurement of the 3D-printed conductive unit cells within the WR-187 rectangular waveguide. The reflection (*R* = |*S*11| 2 ) and transmission (*T* = |*S*21| 2 ) power coefficients are plotted in the frequency range of 4.5–5.5 GHz.

#### **4. Discussion and Conclusions**

We have studied a conductive meta-atom of toroidal topology, meant to enhance the toroidal dipole moment. We have performed simulations of both a periodic freespace metasurface comprised of this meta-atom and, subsequently, we have inserted two such meta-atoms into a metallic, rectangular waveguide setup. The latter configuration has allowed for experimentally verifying the resonance characteristics (lineshape and linewidth) and thus supports the physical claims arising from the theoretical analysis and multipole expansion.

One goal of our work was to demonstrate a practical, easily realizable meta-atom geometry that provides freedom in shaping the electromagnetic response in terms of both linewidth and lineshape. Another important goal was to reveal and highlight the physical mechanisms dictating the resonant response. Importantly, despite the physical interpretation suggested by the unit cell geometry (its toroidal topology), a rigorous multipole expansion has revealed that, although the toroidal dipole moment is indeed strong, the corresponding scattered field is cancelled out exactly by that of the magnetic quadrupole. As a result, the scattering response can be described accurately by simply considering the electric dipole moment alone. Thus, we conclude that a multipole expansion must be performed in all cases, in order to reach safe conclusions regarding the multipole components that mediate the scattering process.

The studied meta-atom allows to control the residual electric dipole moment through the geometric asymmetry between the central and outer gaps. This is a very similar concept to "accidental BIC" or "trapped/broken-symmetry" resonances. It becomes possible to achieve a controllably sharp resonant response, reaching very high radiation quality factors (*Q*rad). This is also possible in finite meta-atom arrays through spatially extended dark (sub-radiant) eigenmodes, as has been discussed in ref. [24]. Note, however, that ultimately the spectral linewidth is dictated by the total quality factor (*Q*tot), which will be limited by resistive losses as well. As a result, pushing for very high *Q*rad values which require slight asymmetry and precise fabrication is not necessary when it cannot be supported by low material losses. Interconnected meta-atoms that do not require a substrate can help to avoid excess losses stemming from the substrate material and excess radiation losses stemming from the breaking of the vertical symmetry.

**Author Contributions:** Conceptualization: E.N.E., O.T. and A.C.T.; methodology, O.T., A.C.T., D.C.Z. and E.N.E.; software: O.T. (implementation of multipole expansion and models); validation: O.T., Z.V. and G.K.; formal analysis, O.T.; investigation: O.T. and Z.V.; resources: M.K., D.C.Z., Z.V. and G.K.; data curation: O.T.; writing—original draft preparation, O.T.; writing—review and editing: Z.V., D.C.Z., M.K., A.C.T., G.K., E.N.E. and O.T.; visualization: O.T.; supervision: E.N.E., M.K. and G.K.; project administration: G.K.; funding acquisition: G.K. and D.C.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the project "METAmaterial-based ENERGY autonomous systems (META-ENERGY) (Project ID 2936)" which is implemented under the 2nd Call for H.F.R.I. "Research Projects to Support Faculty Members & Researchers" funded by the Operational Programme "Competitiveness, Entrepreneurship and Innovation" (NSRF 2014–2020). D.C.Z. acknowledges the support of the project ECS00000024 "Ecosistemi dell'Innovazione"- Rome Technopole of the Italian Ministry of University and Research, public call n. 3277, PNRR-Mission 4, Component 2, Investment 1.5, financed by the European Union, Next Generation EU and of the CNR-FAPESP biennial (2022–2023) bilateral project StReAM "Strongly Resonant All-dielectric Metasurfaces based on quasi-dark and toroidal modes".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The authors thank Varvara Mouzi for their help with silver paste coating of the 3D-printed meta-atoms.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

### *Article* **Bloch Surface Waves in Open Fabry–Perot Microcavities**

**Niccolò Marcucci <sup>1</sup> , Tian-Long Guo <sup>2</sup> , Ségolène Pélisset <sup>2</sup> , Matthieu Roussey <sup>2</sup> , Thierry Grosjean <sup>3</sup> and Emiliano Descrovi 1,\***


**\*** Correspondence: emiliano.descrovi@polito.it

**Abstract:** Thanks to the increasing availability of technologies for thin film deposition, all-dielectric structures are becoming more and more attractive for integrated photonics. As light–matter interactions are involved, Bloch Surface Waves (BSWs) may represent a viable alternative to plasmonic platforms, allowing easy wavelength and polarization manipulation and reduced absorption losses. However, plasmon-based devices operating at an optical and near-infrared frequency have been demonstrated to reach extraordinary field confinement capabilities, with localized mode volumes of down to a few nanometers. Although such levels of energy localization are substantially unattainable with dielectrics, it is possible to operate subwavelength field confinement by employing high-refractive index materials with proper patterning such as, e.g., photonic crystals and metasurfaces. Here, we propose a computational study on the transverse localization of BSWs by means of quasi-flat Fabry–Perot microcavities, which have the advantage of being fully exposed toward the outer environment. These structures are constituted by defected periodic corrugations of a dielectric multilayer top surface. The dispersion and spatial distribution of BSWs' cavity mode are presented. In addition, the hybridization of BSWs with an A exciton in a 2D flake of tungsten disulfide (WS<sup>2</sup> ) is also addressed. We show evidence of strong coupling involving not only propagating BSWs but also localized BSWs, namely, band-edge and cavity modes.

**Keywords:** Bloch Surface Waves; strong coupling; TMD; 2D materials

#### **1. Introduction**

The existence of optical modes strongly confined at the truncation interface of planar dielectric multilayers (one-dimensional photonic crystals, -1DPC) has been known for decades as a result of the pioneering studies by P. Yeh and A. Yariv [1,2] and several later works [3–5]. With the advent of the Plasmonics era, interest toward optical surface modes on planar dielectric stacks has seen a significant boost, mainly due to the possibility of overcoming some the inherent limitations of surface plasmons, especially in sensing applications [6–9]. In the last 20 years, optical surface modes (hereafter called Bloch Surface Waves, or BSW) on flat and patterned dielectric multilayers have been investigated in many frameworks, such as strong light–matter coupling [10–13], integrated [14–18], guided [19–22] and fiber [23,24] optics, label-free [25–28] and fluorescence-based [29–31] sensing, metrology, [32] photon management [33–36], light-driven particle manipulation, [37,38] emitting devices [39,40], microscopy imaging, and spectroscopy [41–43].

Besides the advantages of tuning the BSWs' spectral position, polarization [44], as well as propagation and penetration lengths [45], one of the most intriguing features offered by 1DPC-based platforms is the possibility to exploit the 1DPC surface to set up an accessible framework for controlling light–matter interaction. Furthermore, shallow patterning of the 1DPC surface allows the BSW to be manipulated and possibly confined along transverse directions [40,46], thus providing new degrees of freedom for exploring

**Citation:** Marcucci, N.; Guo, T.-L.; Pélisset, S.; Roussey, M.; Grosjean, T.; Descrovi, E. Bloch Surface Waves in Open Fabry–Perot Microcavities. *Micromachines* **2023**, *14*, 509. https://doi.org/10.3390/ mi14030509

Academic Editors: Luigi Sirleto and Giancarlo C. Righini

Received: 11 January 2023 Revised: 17 February 2023 Accepted: 20 February 2023 Published: 22 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

complex phenomena. Having a photonic mode confined both out-of-plane and in-plane on a surface is advantageous when alleviating some difficulties in technological tasks such as the integration of emitters/absorbers, which are typically buried/grown within the photonic structure itself, where the radiation energy of photonic modes is generally localized. complex phenomena. Having a photonic mode confined both out-of-plane and in-plane on a surface is advantageous when alleviating some difficulties in technological tasks such as the integration of emitters/absorbers, which are typically buried/grown within the photonic structure itself, where the radiation energy of photonic modes is generally localized. In this article, we present a computational work describing some peculiar features of

accessible framework for controlling light–matter interaction. Furthermore, shallow patterning of the 1DPC surface allows the BSW to be manipulated and possibly confined along transverse directions [40,46], thus providing new degrees of freedom for exploring

*Micromachines* **2023**, *14*, x FOR PEER REVIEW 2 of 13

In this article, we present a computational work describing some peculiar features of transverse electric (TE)-polarized BSWs confined within linear resonant structures recalling Fabry–Perot microcavities. In addition, the interaction of cavity BSWs with the A exciton in a WS<sup>2</sup> monolayer laying on top of the multilayered structure is investigated and discussed. WS<sup>2</sup> monolayers display a narrow and intense excitonic resonance at 2.03 eV that is well suited to promote mode hybridization with BSWs. Extending previous results reporting on strong coupling between the WS<sup>2</sup> exciton and BSWs on flat dielectric multilayers [47], we anticipate here that strong coupling can be also observed with BSWs confined within planar cavities covered by WS<sup>2</sup> flakes of down to about 1 µm in transverse size. transverse electric (TE)-polarized BSWs confined within linear resonant structures recalling Fabry–Perot microcavities. In addition, the interaction of cavity BSWs with the A exciton in a WS<sup>2</sup> monolayer laying on top of the multilayered structure is investigated and discussed. WS2 monolayers display a narrow and intense excitonic resonance at 2.03 eV that is well suited to promote mode hybridization with BSWs. Extending previous results reporting on strong coupling between the WS<sup>2</sup> exciton and BSWs on flat dielectric multilayers [47], we anticipate here that strong coupling can be also observed with BSWs confined within planar cavities covered by WS<sup>2</sup> flakes of down to about 1 m in transverse size.

#### **2. Computational Model 2. Computational Model**

The exemplary model for the 1D photonic crystal considered in this work is constituted by a stack of five alternating pairs of TiO<sup>2</sup> and SiO<sup>2</sup> layers on a semi-infinite glass half-space (refractive index *nglass* = 1.5). The TiO<sup>2</sup> and SiO<sup>2</sup> layers have thickness 85 nm and 110 nm, respectively, so that a main forbidden band is crossing the light line at energies of about 2 eV. On top of the stack, an additional pair of SiO<sup>2</sup> (90 nm) and TiO<sup>2</sup> (40 nm) layers is introduced, leading to a total of N<sup>L</sup> = 12 layers. These two layers are aimed at tailoring the position of the BSW across the forbidden band, in the energy–momentum space. In addition, the upmost TiO<sup>2</sup> layer is intended as a functional medium carrying a topographic modulation. In fact, a linear grating acting as a Distributed Bragg Reflector (DBR) is engraved therein, with a modulation height *h* = 25 nm, period ΛDBR = 275 nm, and fill factor 0.5. A linear defect (hereafter called a "spacer") with a width *W* is eventually introduced within the DBR, thus resulting in a linear Fabry–Perot cavity. A schematic of the structure is presented in Figure 1a. The exemplary model for the 1D photonic crystal considered in this work is constituted by a stack of five alternating pairs of TiO<sup>2</sup> and SiO<sup>2</sup> layers on a semi-infinite glass half-space (refractive index = 1.5). The TiO<sup>2</sup> and SiO<sup>2</sup> layers have thickness 85 nm and 110 nm, respectively, so that a main forbidden band is crossing the light line at energies of about 2 eV. On top of the stack, an additional pair of SiO<sup>2</sup> (90 nm) and TiO<sup>2</sup> (40 nm) layers is introduced, leading to a total of NL = 12 layers. These two layers are aimed at tailoring the position of the BSW across the forbidden band, in the energy–momentum space. In addition, the upmost TiO<sup>2</sup> layer is intended as a functional medium carrying a topographic modulation. In fact, a linear grating acting as a Distributed Bragg Reflector (DBR) is engraved therein, with a modulation height ℎ = 25 nm, period ΛDBR = 275 nm, and fill factor 0.5. A linear defect (hereafter called a "spacer") with a width is eventually introduced within the DBR, thus resulting in a linear Fabry–Perot cavity. A schematic of the structure is presented in Figure 1a.

**Figure 1.** (**a**) Sketch of the patterned 1DPC in the conical diffraction mounting. The topographic modulation is oriented along the x-direction. Light is incident from the glass substrate at an angle of incidence ; (**b**) refractive index dispersion of the materials considered in this work. **Figure 1.** (**a**) Sketch of the patterned 1DPC in the conical diffraction mounting. The topographic modulation is oriented along the x-direction. Light is incident from the glass substrate at an angle of incidence *θ*; (**b**) refractive index dispersion of the materials considered in this work.

As will be apparent from the following, the DBR is aimed at perturbing the BSW dispersion by opening a bandgap at the crossing point of the BSW dispersion with the boundary of the first Brillouin zone defined by the DBR. The bandgap width is determined by the BSW (effective) refractive index contrast associated to the corrugation. Generally As will be apparent from the following, the DBR is aimed at perturbing the BSW dispersion by opening a bandgap at the crossing point of the BSW dispersion with the boundary of the first Brillouin zone defined by the DBR. The bandgap width is determined by the BSW (effective) refractive index contrast associated to the corrugation. Generally speaking, the larger the topographic height modulation, the wider the BSW bandgap. However, practical limitations to the maximum attainable width of the BSW bandgap exist. For example, if the modulation height or the top layer thickness are too large, the dielectric

loading/unloading effect may cause BSWs to disappear in either the trenches (low effective refractive index) or the ridges (high effective refractive index), thus preventing a BSW bandgap to open. The geometry proposed here represents a (not unique) solution that avoids the occurrence of such an issue, while allowing the formation of a reasonably wide BSW bandgap for observing cavity BSWs and related hybridization with the WS<sup>2</sup> exciton.

This kind of dielectric multilayer can be fabricated by means of standard techniques, including RF sputtering [16], chemical vapor deposition (CVD) and plasma-enhanced CVD [48] plasma ion-assisted vacuum evaporation (PIAD) [49], and atomic layer deposition (ALD) [50]. In particular, the refractive indexes used to design the present stack are taken from ellipsometric data collected from ALD-deposited materials and plotted in Figure 1b. In addition, the refractive index of a WS<sup>2</sup> monolayer is plotted as well [51].

Calculations were performed by means of a freely available implementation of the Rigorous Coupled Wave Analysis (RCWA) RETICOLO [52,53]. The RCWA supercell includes the spacer surrounded by 180 DBR periods on each side. The periodic corrugation has an associated Bragg grating vector oriented along the x-axis, i.e., KDBR = 2*π*·Λ −1 DBR, 0, 0 . The conical mounting configuration is considered so that the incoming radiation can be made incident onto the bottom interface of the multilayer at an incident angle *θ* with respect to the normal (z-axis) and at an angle *ϕ* with respect to the x-axis. At *ϕ* = 0 o , the incident light is parallel to the KDBR and diffracted accordingly. Illumination is always provided from the glass substrate to reach momentum values beyond the light line in air, as required for BSW coupling. RCWA is a Fourier-based method involving a Rayleigh expansion of the diffracted field as well as a Fourier decomposition of the structure harmonics. For this reason, it is crucial to define the proper number of Fourier terms to be retained in the calculation. The case of the spacer surrounded by DBRs is addressed using N = 1200 Fourier terms. In the convergence plot in Figure 2, the spectral position of the BSW cavity mode and the corresponding reflectivity are shown as a function of the number of Fourier terms considered. Such a supercell contains a number of DBR periods large enough to avoid significant cross-talk effects between adjacent supercells occurring in the momentum– energy region wherein the cavity mode is located. Instead, the simpler case of a purely periodic corrugation (i.e., no spacer) is modeled as an elementary cell containing a single DBR period and N = 60 Fourier terms for assuring reliable results. *Micromachines* **2023**, *14*, x FOR PEER REVIEW 4 of 13

**Figure 2.** Reflectivity map (, ℏ) and **s**pectral position of the BSW resonance (cyan circles) in an exemplary cavity with spacer width = 500 nm and 360 DBR periods overall at = 0 o incidence. N indicates the number of Fourier terms retained in the RCWA code. The reflectivity values at the resonance center (red squares) are also plotted. The dashed pink line indicates the number of harmonics used in this work. **Figure 2.** Reflectivity map *R*(*N*, }*ω*) and spectral position of the BSW resonance (cyan circles) in an exemplary cavity with spacer width *W* = 500 nm and 360 DBR periods overall at *ϕ* = 0 o incidence. N indicates the number of Fourier terms retained in the RCWA code. The reflectivity values at the resonance center (red squares) are also plotted. The dashed pink line indicates the number of harmonics used in this work.

The multilayer supports TE-polarized Bloch Surface Waves whose dispersion changes from the near-infrared to the visible spectrum. The dispersion of the BSWs on the flat multilayer can be inferred upon calculation of the reflectivity *R*(*β*, }*ω*), wherein the effective

**Figure 3.** (**a**) Calculated TE−polarized log(1 − (, ℏ)) of the flat 1DPC, with illumination from the bottom glass substrate at an incidence angle such that = glass ∙ sin; (**b**) intensity profile of

The case of the planar stack corrugated with a purely periodic surface modulation (no defect spacer) is considered first. In Figure 4a, the energy and angularly resolved full

<sup>2</sup> = (ℏ). Reflectivity values refer to the diffraction efficiency of the 0

The periodic modulation along the x-direction results in a folding of the BSW dispersion and the formation of a forbidden band that is dispersed in energy. To facilitate the interpretation, the BSW folding and the corresponding energy gap at = 0 (corresponding to = 0, i.e., a direction parallel to the grating vector K) is presented in Figure 4b.

opening at about 2 eV. As pointed out elsewhere [40,54], the energy gap width depends on the effective refractive index contrast introduced by the periodic modulation. More interestingly, along directions different than K, the gap shifts to higher energies. Band

Brillouin zone boundary, K/2) with the full BSW dispersion, as shown in Figure 4c. The forbidden band can be better appreciated by projecting said BSW dispersion onto the

, <sup>y</sup> , )

st order results in a gap

) −1

(first

) = ℏ ∙ (Λ ∙

th

dispersion of the TE-polarized BSWs is shown as a 3D surface whose points (<sup>x</sup>

are extracted from the calculated reflectivity (,, ℏ). There, = ∙ cos and = ∙ sin, and they must match the BSW effective refractive index at any BSW energy ℏ,

the BSW at ℏ = 2 eV, normalized to the amplitude of the incident wave.

The folding of the BSW dispersion due to diffraction at the −1

edges are identified at the intersection of the surface (

**3. Results**

i.e., √

<sup>2</sup> +

order at different illumination conditions.

refractive index is *β* = *n*glass· sin *θ* and }*ω* is the photon energy of the TE-polarized incident radiation. Without loss of generality, the BSW dispersion can be identified by the pairs (*β*BSW, *ω*) corresponding to the bright line in Figure 3a. The TiO<sup>2</sup> and SiO<sup>2</sup> layers have been given a small imaginary refractive index *k*SiO<sup>2</sup> = *k*TiO<sup>2</sup> = 10−<sup>3</sup> to introduce enough losses to make the BSW dip visible on the calculated reflectivity maps. In Figure 3a,b, the log(1 − *R*(*β*, }*ω*)) map and the normalized intensity profile of the BSW at }*ω* = 2 eV are respectively shown. **Figure 2.** Reflectivity map (, ℏ) and **s**pectral position of the BSW resonance (cyan circles) in an exemplary cavity with spacer width = 500 nm and 360 DBR periods overall at = 0 o incidence. N indicates the number of Fourier terms retained in the RCWA code. The reflectivity values at the resonance center (red squares) are also plotted. The dashed pink line indicates the number of harmonics used in this work.

*Micromachines* **2023**, *14*, x FOR PEER REVIEW 4 of 13

**Figure 3.** (**a**) Calculated TE−polarized log(1 − (, ℏ)) of the flat 1DPC, with illumination from the bottom glass substrate at an incidence angle such that = glass ∙ sin; (**b**) intensity profile of the BSW at ℏ = 2 eV, normalized to the amplitude of the incident wave. **Figure 3.** (**a**) Calculated TE−polarized log(1 − *R*(*β*, }*ω*)) of the flat 1DPC, with illumination from the bottom glass substrate at an incidence angle *θ* such that *β* = *n*glass· sin *θ*; (**b**) intensity profile of the BSW at }*ω* = 2 eV, normalized to the amplitude of the incident wave.

#### **3. Results 3. Results**

The case of the planar stack corrugated with a purely periodic surface modulation (no defect spacer) is considered first. In Figure 4a, the energy and angularly resolved full dispersion of the TE-polarized BSWs is shown as a 3D surface whose points (<sup>x</sup> , <sup>y</sup> , ) are extracted from the calculated reflectivity (,, ℏ). There, = ∙ cos and = ∙ sin, and they must match the BSW effective refractive index at any BSW energy ℏ, i.e., √ <sup>2</sup> + <sup>2</sup> = (ℏ). Reflectivity values refer to the diffraction efficiency of the 0 th order at different illumination conditions. The case of the planar stack corrugated with a purely periodic surface modulation (no defect spacer) is considered first. In Figure 4a, the energy and angularly resolved full dispersion of the TE-polarized BSWs is shown as a 3D surface whose points *β*x, *β*y, *ω* are extracted from the calculated reflectivity *R*(*β*, *ϕ*, }*ω*). There, *β<sup>x</sup>* = *β*· cos *ϕ* and *β<sup>y</sup>* = *β*· sin *ϕ*, and they must match the BSW effective refractive index at any BSW energy }*ω*, i.e., <sup>q</sup> *β* 2 *<sup>x</sup>* + *β* 2 *<sup>y</sup>* = *βBSW*(}*ω*). Reflectivity values refer to the diffraction efficiency of the 0th order at different illumination conditions.

The periodic modulation along the x-direction results in a folding of the BSW dispersion and the formation of a forbidden band that is dispersed in energy. To facilitate the interpretation, the BSW folding and the corresponding energy gap at = 0 (corresponding to = 0, i.e., a direction parallel to the grating vector K) is presented in Figure 4b. The folding of the BSW dispersion due to diffraction at the −1 st order results in a gap opening at about 2 eV. As pointed out elsewhere [40,54], the energy gap width depends on the effective refractive index contrast introduced by the periodic modulation. More interestingly, along directions different than K, the gap shifts to higher energies. Band edges are identified at the intersection of the surface ( ) = ℏ ∙ (Λ ∙ ) −1 (first Brillouin zone boundary, K/2) with the full BSW dispersion, as shown in Figure 4c. The forbidden band can be better appreciated by projecting said BSW dispersion onto the The periodic modulation along the x-direction results in a folding of the BSW dispersion and the formation of a forbidden band that is dispersed in energy. To facilitate the interpretation, the BSW folding and the corresponding energy gap at *β<sup>y</sup>* = 0 (corresponding to *ϕ* = 0, i.e., a direction parallel to the grating vector K*DBR*) is presented in Figure 4b. The folding of the BSW dispersion due to diffraction at the −1 st order results in a gap opening at about 2 eV. As pointed out elsewhere [40,54], the energy gap width depends on the effective refractive index contrast introduced by the periodic modulation. More interestingly, along directions different than K*DBR*, the gap shifts to higher energies. Band edges are identified at the intersection of the surface *σ*(*βx*) = *π*}*c*·(Λ*DBR*·*βx*) −1 (first Brillouin zone boundary, K*DBR*/2) with the full BSW dispersion, as shown in Figure 4c. The forbidden band can be better appreciated by projecting said BSW dispersion onto the *βy*, }*ω* plane (Figure 4d). A maximum bandgap of about 150 meV at *β<sup>y</sup>* = 0 is found. For larger values of *β<sup>y</sup>*  and higher energies, the bandgap is narrowed down.

At each bandgap edge, two counterpropagating BSWs are produced, whose interference pattern has a spatial period equal to the DBR period ΛDBR. In Figure 5a,b, the normalized intensity distributions of the BSWs at the lower and the upper band are respectively shown. The field localization on either the ridges (high refractive index regions) or trenches (low refractive index regions) reflects the different energies associated to the two band edge modes. Calculations are performed at }*ω<sup>L</sup>* = 1.953 eV and }*ω<sup>U</sup>* = 2.1 eV for the lower and the upper band edge modes, respectively. In both cases, *β<sup>x</sup>* is taken on

the boundary of first Brillouin zone, i.e., *β<sup>x</sup>* = *πc*·(Λ*DBR*·*ωL*,*U*) −1 , while *β<sup>y</sup>* = 0. Field amplitudes are normalized to the incident radiation amplitude. (, ℏ) plane (Figure 4d). A maximum bandgap of about 150 meV at = 0 is found. For larger values of || and higher energies, the bandgap is narrowed down.

**Figure 4.** (**a**) Dispersion of a TE−polarized BSW represented as a 3D surface from the calculated reflectivity (, , ℏ) of the 1DPC patterned with the periodic DBR. Only the region beyond the light-line in air √ <sup>2</sup> + <sup>2</sup> ≥ 1 is considered; BSW dispersion on (**b**) the = 0 plane and (**c**) the boundary of the first Brillouin zone defined by the shaded surface = ℏ ∙ (Λ ∙ ) −1 ; (**d**) log (1 − (, , ℏ)) map on the surface, where = ∙ (Λ ∙ ) −1 . **Figure 4.** (**a**) Dispersion of a TE−polarized BSW represented as a 3D surface from the calculated reflectivity *R βx*, *βy*, }*ω* of the 1DPC patterned with the periodic DBR. Only the region beyond the light-line in air <sup>q</sup> *β* 2 *<sup>x</sup>* + *β* 2 *<sup>y</sup>* ≥ 1 is considered; BSW dispersion on (**b**) the *β<sup>y</sup>* = 0 plane and (**c**) the boundary of the first Brillouin zone defined by the shaded surface *σ* = *π*}*c*·(Λ*DBR*·*βx*) −1 ; (**d**) log 1 − *R βx*,*βy*, }*ω* map on the *<sup>σ</sup>* surface, where *<sup>β</sup><sup>x</sup>* <sup>=</sup> *<sup>π</sup>c*·(Λ*DBR*·*ω*) −1 . *Micromachines* **2023**, *14*, x FOR PEER REVIEW 6 of 13

bidden band created by the BSW folding. In Figure 6a, the BSW cavity mode corresponding to a spacer width = 500 nm is visible as a narrow band starting at about 2.03 eV. In comparison to the spacer mode in a conventional planar microcavity, the BSW cavity mode follows a dispersion curve depending on the momentum component that is trans-**Figure 5.** Intensity field distribution |〈| (, )| 2 〉| of (**a**) the lower band BSW, calculated at ℏ = 1.953 eV, = ∙ (Λ ∙ ) −1 , = 0; (**b**) the upper band BSW, calculated at ℏ = 2.1 eV, = ∙ (Λ ∙ ) −1 , = 0. **Figure 5.** Intensity field distribution < *Ey*(*x*, *z*) <sup>2</sup> >  of (**a**) the lower band BSW, calculated at }*ω<sup>L</sup>* = 1.953 eV, *β<sup>x</sup>* = *πc*·(Λ*DBR*·*ωL*) −1 , *β<sup>y</sup>* = 0; (**b**) the upper band BSW, calculated at }*ω<sup>U</sup>* = 2.1 eV, *β<sup>x</sup>* = *πc*·(Λ*DBR*·*ωU*) −1 , *β<sup>y</sup>* = 0.

verse to the direction of the refractive index modulation. In our case, since the pattern modulation is developing in the x-direction, the BSW cavity mode is energy dispersed as a function of . It is worth noting that, for || > 0.6, the cavity mode tends to merge with the lower band edge mode, thus loosing spatial localization. Instead, the BSW cavity mode is dispersionless along the direction (Figure 6b). The spectral position of the cavity mode can be modulated across the forbidden band by varying the spacer width *W* (Figure 6c). In Figure 6d, a cross-sectional view of the spatial near-field intensity of the BSW cavity mode at ℏ = 2.03 eV is illustrated. As in a Fabry–Perot planar microcavity, the energy of the cavity mode is preferentially localized in the defect region. However, the very low refractive index contrast experienced by the BSW prevents a strong confinement within the spacer only, with a significant penetration If a linear defect is introduced into the DBR, a cavity mode appears within the forbidden band created by the BSW folding. In Figure 6a, the BSW cavity mode corresponding to a spacer width *W* = 500 nm is visible as a narrow band starting at about 2.03 eV. In comparison to the spacer mode in a conventional planar microcavity, the BSW cavity mode follows a dispersion curve depending on the momentum component that is transverse tothe direction of the refractive index modulation. In our case, since the pattern modulation

**Figure 6.** Dispersion of a TE−polarized BSW on the 1DPC patterned with a cavity spacer = 500 nm surrounded by 180 DBR periods on each side, calculated on (**a**) the surface, (**b**) the =0 plane; (**c**) spectral position (purple circles) and corresponding spectra of the cavity mode across the

nm (pink line), 430 nm (green line), 450 nm (gray line), 470 nm (yellow line), 500 nm (black line), 550

nm (blue line), 600 nm (red line); (**d**) near-field intensity of the BSW cavity mode 〈|

−1

, 0) for several spacer values , i.e. 400

(, )| 2 〉 for

forbidden band calculated at (, ) = ( ∙ (Λ ∙ )

= 500 nm, at ℏ = 2.032 eV, = 0.

of the mode profile in the surrounding periodic corrugations.

is developing in the x-direction, the BSW cavity mode is energy dispersed as a function of *βy*. It is worth noting that, for *β<sup>y</sup>*  > 0.6, the cavity mode tends to merge with the lower band edge mode, thus loosing spatial localization. Instead, the BSW cavity mode is dispersionless along the *β<sup>x</sup>* direction (Figure 6b). Fabry–Perot planar microcavity, the energy of the cavity mode is preferentially localized in the defect region. However, the very low refractive index contrast experienced by the BSW prevents a strong confinement within the spacer only, with a significant penetration of the mode profile in the surrounding periodic corrugations.

The spectral position of the cavity mode can be modulated across the forbidden band by varying the spacer width *W* (Figure 6c). In Figure 6d, a cross-sectional view of the spatial near-field intensity of the BSW cavity mode at ℏ = 2.03 eV is illustrated. As in a

(, )| 2

〉| of (**a**) the lower band BSW, calculated at ℏ =

, = 0; (**b**) the upper band BSW, calculated at ℏ = 2.1 eV,

*Micromachines* **2023**, *14*, x FOR PEER REVIEW 6 of 13

**Figure 5.** Intensity field distribution |〈|

−1

) −1

, = 0.

1.953 eV, = ∙ (Λ ∙

= ∙ (Λ ∙ )

**Figure 6.** Dispersion of a TE−polarized BSW on the 1DPC patterned with a cavity spacer = 500 nm surrounded by 180 DBR periods on each side, calculated on (**a**) the surface, (**b**) the =0 plane; (**c**) spectral position (purple circles) and corresponding spectra of the cavity mode across the forbidden band calculated at (, ) = ( ∙ (Λ ∙ ) −1 , 0) for several spacer values , i.e. 400 nm (pink line), 430 nm (green line), 450 nm (gray line), 470 nm (yellow line), 500 nm (black line), 550 nm (blue line), 600 nm (red line); (**d**) near-field intensity of the BSW cavity mode 〈| (, )| 2 〉 for = 500 nm, at ℏ = 2.032 eV, = 0. **Figure 6.** Dispersion of a TE−polarized BSW on the 1DPC patterned with a cavity spacer *W* = 500 nm surrounded by 180 DBR periods on each side, calculated on (**a**) the *σ* surface, (**b**) the *β<sup>y</sup>* = 0 plane; (**c**) spectral position (purple circles) and corresponding spectra of the cavity mode across the forbidden band calculated at *βx*, *βy* = *πc*·(Λ*DBR*·*ω*) −1 , 0 for several spacer values *W*, i.e. 400 nm (pink line), 430 nm (green line), 450 nm (gray line), 470 nm (yellow line), 500 nm (black line), 550 nm (blue line), 600 nm (red line); (**d**) near-field intensity of the BSW cavity mode *Ey*(*x*, *z*) 2 for *W* = 500 nm, at }*ω* = 2.032 eV, *β<sup>y</sup>* = 0.

The spectral position of the cavity mode can be modulated across the forbidden band by varying the spacer width *W* (Figure 6c). In Figure 6d, a cross-sectional view of the spatial near-field intensity of the BSW cavity mode at }*ω* = 2.03 eV is illustrated. As in a Fabry–Perot planar microcavity, the energy of the cavity mode is preferentially localized in the defect region. However, the very low refractive index contrast experienced by the BSW prevents a strong confinement within the spacer only, with a significant penetration of the mode profile in the surrounding periodic corrugations.

Very recently, multilayers sustaining BSWs coated with organic/inorganic thin layers exhibiting excitonic resonances have been exploited as novel platforms to enhance light– matter interactions [12,13,55]. In this framework, 2D materials (e.g., Transition Metal Dichalcogenides, orTMD monolayers) are particularly attractive because of the possibility to operate polariton manipulation at room temperature [56,57]. Note that the TE-polarization of BSW facilitates the interaction with excitons whose transition dipole lays on the plane, as occurs, for example, in WS<sup>2</sup> monolayers [58]. Strong coupling effects between BSWs and excitons have been experimentally studied on flat or periodically patterned multilayers, wherein the resulting polariton field exhibits a spatially delocalized distribution.

When our flat 1DPC model is coated with a tungsten disulfide (WS2) monolayer (Figure 7a), the BSW dispersion exhibits the typical anti-crossing behavior illustrated

in Figure 7b. The upper and lower BSW polariton branches repel each other via a vacuum Rabi splitting Ω that can be estimated by means of a simple coupled harmonic oscillator model. The system Hamiltonian can be expressed as a 2 × 2 matrix [12] with non-zero off-diagonal elements H = *EBSW* Ω/2 Ω/2 *E<sup>X</sup>* , where *E<sup>X</sup>* is the exciton energy and *EBSW*(*βx*) is the dispersion of the uncoupled BSW. Hamiltonian eigenvalues are: E*up* = 1 2 *EBSW* + *E<sup>X</sup>* + q (*EBSW* − *EX*) <sup>2</sup> + Ω<sup>2</sup> ; E*l p* = <sup>1</sup> 2 *EBSW* + *E<sup>X</sup>* − q (*EBSW* − *EX*) <sup>2</sup> + Ω<sup>2</sup> , which correspond to the dispersions of the Upper and the Lower Polariton Branch (UPB and LPB), respectively. A fitting procedure is used to estimate Ω and *EX*. Worth specifying is that, to consider the dielectric loading effect introduced by the WS<sup>2</sup> monolayer, the values for *EBSW*(*βx*) calculated for the bare 1DPC need be adjusted (redshift) by an additional energy term *E*<sup>0</sup> that is also obtained from the fit. In Figure 7b, the blue and the red dashed lines represent the lower and upper polariton dispersion, respectively, the yellow line represents the exciton energy *E<sup>X</sup>* = 2.032 eV, and the green line represents the uncoupled BSW dispersion, redshift-corrected. From the fit, we obtain a Rabi splitting Ω = 52 meV, which is similar to previously published results on a different BSW platform [13] ure 7a), the BSW dispersion exhibits the typical anti-crossing behavior illustrated in Figure 7b. The upper and lower BSW polariton branches repel each other via a vacuum Rabi splitting Ω that can be estimated by means of a simple coupled harmonic oscillator model. The system Hamiltonian can be expressed as a 2 × 2 matrix [12] with non-zero offdiagonal elements H = ( Ω⁄2 Ω⁄2 ), where is the exciton energy and ( ) is the dispersion of the uncoupled BSW. Hamiltonian eigenvalues are: E = 1 2 ( + + √( − ) <sup>2</sup> + Ω<sup>2</sup>); E = 1 2 ( + − √( − ) <sup>2</sup> + Ω<sup>2</sup>), which correspond to the dispersions of the Upper and the Lower Polariton Branch (UPB and LPB),respectively. A fitting procedure is used to estimate Ω and . Worth specifying is that, to consider the dielectric loading effect introduced by the WS<sup>2</sup> monolayer, the values for ( ) calculated for the bare 1DPC need be adjusted (redshift) by an additional energy term <sup>0</sup> that is also obtained from the fit. In Figure 7b, the blue and the red dashed lines represent the lower and upper polariton dispersion, respectively, the yellow line represents the exciton energy = 2.032 eV, and the green line represents the uncoupled BSW dispersion, redshift-corrected. From the fit, we obtain a Rabi splitting Ω = 52 meV, which is similar to previously published results on a different BSW platform [13]

Very recently, multilayers sustaining BSWs coated with organic/inorganic thin layers exhibiting excitonic resonances have been exploited as novel platforms to enhance light– matter interactions [12,13,55]. In this framework, 2D materials (e.g., Transition Metal Dichalcogenides, orTMD monolayers) are particularly attractive because of the possibility to operate polariton manipulation at room temperature [56,57]. Note that the TE-polarization of BSW facilitates the interaction with excitons whose transition dipole lays on the plane, as occurs, for example, in WS<sup>2</sup> monolayers [58]. Strong coupling effects between BSWs and excitons have been experimentally studied on flat or periodically patterned multilayers, wherein the resulting polariton field exhibits a spatially delocalized distribution. When our flat 1DPC model is coated with a tungsten disulfide (WS2) monolayer (Fig-

*Micromachines* **2023**, *14*, x FOR PEER REVIEW 7 of 13

**Figure 7.** (**a**) Sketch of the z-cross section of the flat 1DPC topped with a uniform WS2 monolayer with an indefinite extension; (**b**) corresponding dispersion of BSW polariton; (**c**) sketch of the z−cross **Figure 7.** (**a**) Sketch of the z-cross section of the flat 1DPC topped with a uniform WS<sup>2</sup> monolayer with an indefinite extension; (**b**) corresponding dispersion of BSW polariton; (**c**) sketch of the z−cross section of the 1DPC patterned with the cavity topped with a WS<sup>2</sup> monolayer. The flake has a lateral size *Sf lake* = 1.6 µm and a thickness 0.65 nm; (**d**) corresponding dispersion of the BSW cavity polariton branches (UPBcav and LPBcav).

When the BSW cavity is considered, the BSW–exciton interaction occurs over a spatial region that is substantially limited by the transverse size of the cavity mode. The sketch in Figure 7c shows the case of a WS<sup>2</sup> monolayer deposited on top of the BSW cavity, extending over the spacer and the first surrounding DBR period, with a total size *Sf lake* = 1.05 µm and a thickness of 0.65 nm [59]. Similar to cavity modes in a planar microcavity, the cavity BSW hybridizes with the excitonic resonance, resulting in a mode dispersion splitting into upper and lower polariton branches (UPBcav and LPBcav), as shown in Figure 7d. The Rabi splitting estimated using the fitting procedure described above is found to be Ω = 32.1 meV. Given the half-widths of the BSW cavity mode γ*BSW* =6 meV and the A exciton of WS<sup>2</sup> γ*<sup>X</sup>* =11 meV [60,61], the condition for the occurrence of strong coupling, i.e., Ω > <sup>1</sup> 2 (γ*BSW* + γ*X*), is fulfilled. 7d. The Rabi splitting estimated using the fitting procedure described above is found to be Ω = 32.1 meV. Given the half-widths of the BSW cavity mode γ =6 meV and the A exciton of WS2 γ =11 meV [60,61], the condition for the occurrence of strong coupling, i.e., Ω > 1 2 (γ + γ ), is fulfilled.

section of the 1DPC patterned with the cavity topped with a WS2 monolayer. The flake has a lateral size = 1.6 m and a thickness 0.65 nm; (**d**) corresponding dispersion of the BSW cavity po-

When the BSW cavity is considered, the BSW–exciton interaction occurs over a spatial region that is substantially limited by the transverse size of the cavity mode. The sketch in Figure 7c shows the case of a WS<sup>2</sup> monolayer deposited on top of the BSW cavity, extending over the spacer and the first surrounding DBR period, with a total size = 1.05 m and a thickness of 0.65 nm [59]. Similar to cavity modes in a planar microcavity, the cavity BSW hybridizes with the excitonic resonance, resulting in a mode dispersion splitting into upper and lower polariton branches (UPBcav and LPBcav), as shown in Figure

*Micromachines* **2023**, *14*, x FOR PEER REVIEW 8 of 13

lariton branches (UPBcav and LPBcav).

Losses due to the WS<sup>2</sup> exciton resonance affect the mode quality of the two polariton branches. Figure 8a shows a calculated quality factor Q = *<sup>ω</sup>*0/*δω* for the bare BSW cavity mode, UPBcav, and LPBcav in multilayers made of 8 to 16 layers. For all three modes, resonance frequencies *ω*<sup>0</sup> are recorded on the corresponding dispersion curve, at *β<sup>y</sup>* = 0. While an increase in Q with the number of layers N<sup>L</sup> is ascribed to a decrease in leaky power tunnelling through the 1DPC, the absorption by WS<sup>2</sup> is responsible for the lower quality of both polariton branches, as compared to the bare cavity mode. Furthermore, the upper branch polariton is more severely attenuated, as also observed experimentally in an analogous system involving strong coupling between ZnO excitons and BSWs [12]. Intensity distributions calculated for UPBcav and LPBcav at *β<sup>y</sup>* = 0 are localized in the neighborhood of the spacer, similar to the bare cavity mode, with a substantially reduced enhancement due to losses (Figure 8b,c). Losses due to the WS<sup>2</sup> exciton resonance affect the mode quality of the two polariton branches. Figure 8a shows a calculated quality factor Q = 0 <sup>⁄</sup> for the bare BSW cavity mode, UPBcav, and LPBcav in multilayers made of 8 to 16 layers. For all three modes, resonance frequencies <sup>0</sup> are recorded on the corresponding dispersion curve, at = 0. While an increase in Q with the number of layers N<sup>L</sup> is ascribed to a decrease in leaky power tunnelling through the 1DPC, the absorption by WS<sup>2</sup> is responsible for the lower quality of both polariton branches, as compared to the bare cavity mode. Furthermore, the upper branch polariton is more severely attenuated, as also observed experimentally in an analogous system involving strong coupling between ZnO excitons and BSWs [12]. Intensity distributions calculated for UPBcav and LPBcav at = 0 are localized in the neighborhood of the spacer, similar to the bare cavity mode, with a substantially reduced enhancement due to losses (Figure 8b,c).

**Figure 8.** (**a**) Quality factor of upper (red diamonds) and lower (blue squares) BSW cavity polariton branches and the bare cavity mode (black circles) as a function of the overall number of layers N<sup>L</sup> in **Figure 8.** (**a**) Quality factor of upper (red diamonds) and lower (blue squares) BSW cavity polariton branches and the bare cavity mode (black circles) as a function of the overall number of layers N<sup>L</sup> in the 1DPC; intensity field distribution < *Ey*(*x*, *z*) <sup>2</sup> >  of (**b**) the lower cavity polariton LPBcav, calculated at }*ωLP* = 2 eV, *β<sup>x</sup>* = *πc*·(Λ*DBR*·*ωLP*) −1 , *β<sup>y</sup>* = 0; (**c**) the upper cavity polariton UPBcav, calculated at }*ωUP* = 2.04 eV, *β<sup>x</sup>* = *πc*·(Λ*DBR*·*ωUP*) −1 , *β<sup>y</sup>* = 0. Both intensity distributions refer to an N<sup>L</sup> = 12 1DPC.

It is interesting to note that the BSW hybridization presented above involves the cavity mode only. This is due to the limited size of the WS<sup>2</sup> monolayer, which overlaps mainly with the cavity spacer, where most of the energy is localized. Instead, if the whole patterned 1DPC is topped with an extended monolayer (Figure 9a), the BSW hybridization is found to occur also with the BSW band edge modes, as shown in Figure 9b,c, resulting in upper and lower band edge-polariton branches (UPBbe and LPBbe). Here, we find a Rabi splitting Ω = 37 meV for the BSW cavity polariton, and Ω = 37.6 meV for the lower BSW band edge polariton. The different spatial distribution of the band edge and the cavity modes enables us to further modulate the photonic mode–exciton interaction by controlling the size and position of the WS<sup>2</sup> monolayer transferred onto the BSW platform. upper and lower band edge-polariton branches (UPBbe and LPBbe). Here, we find a Rabi splitting Ω = 37 meV for the BSW cavity polariton, and Ω = 37.6 meV for the lower BSW band edge polariton. The different spatial distribution of the band edge and the cavity modes enables us to further modulate the photonic mode–exciton interaction by controlling the size and position of the WS<sup>2</sup> monolayer transferred onto the BSW platform.

(, )| 2

−1

It is interesting to note that the BSW hybridization presented above involves the cavity mode only. This is due to the limited size of the WS2monolayer, which overlaps mainly with the cavity spacer, where most of the energy is localized. Instead, if the whole patterned 1DPC is topped with an extended monolayer (Figure 9a), the BSW hybridization is found to occur also with the BSW band edge modes, as shown in Figure 9b,c, resulting in

−1

〉| of (**b**) the lower cavity polariton LPBcav, calcu-

, = 0. Both intensity distributions refer to an

, = 0; (**c**) the upper cavity polariton UPBcav, calcu-

*Micromachines* **2023**, *14*, x FOR PEER REVIEW 9 of 13

the 1DPC; intensity field distribution |〈|

N<sup>L</sup> = 12 1DPC.

lated at ℏ = 2 eV, = ∙ (Λ ∙ )

lated at ℏ = 2.04 eV, = ∙ (Λ ∙ )

**Figure 9.** (**a**) Sketch of the z−cross section of the 1DPC patterned with the cavity topped with an indefinitely large WS2 monolayer; corresponding dispersions of hybridized (**b**) BSW cavity polariton mode UPBcav, LPBcav, and (**c**) BSW band edge-polariton mode UPBbe, LPBbe. **Figure 9.** (**a**) Sketch of the z−cross section of the 1DPC patterned with the cavity topped with an indefinitely large WS<sup>2</sup> monolayer; corresponding dispersions of hybridized (**b**) BSW cavity polariton mode UPBcav, LPBcav, and (**c**) BSW band edge-polariton mode UPBbe, LPBbe.

#### **4. Discussion 4. Discussion**

Fostered by an increasing interest toward surface waves on dielectric photonic platforms, we investigated some aspects of the transverse confinement of BSWs in open Fabry–Perot microcavities. Despite the low effective refractive index contrast produced by the ultra-shallow patterning/modulations of the multilayer surface, BSWs can be resonantly confined within sub-micrometric regions. As an extension of the findings presented here, an omnidirectional confinement on the multilayer plane can be obtained by introducing axis-symmetric corrugations [40] or dielectric ridges [46]. The design of the surface modulation is addressed in conjunction with the multilayer design, as dielectric loading/unloading effects can significantly influence the BSW dispersion, resulting in momentum mismatch with the pattern. In the exemplary cases presented here, we designed an open Fabry–Perot cavity by employing materials and geometries that can be realistically handled in clean-room fabrication processes. The cavity mode follows a parabolic energy dispersion as a function of the transverse in-plane momentum (∝ ), thus suggesting a two-dimensional analogy to the case of 3D stacked planar microcavities. Cavity Q factors generally increase with the number of layers within the 1DPC due to reduced leakage losses. Although this feature may look advantageous when a BSW is near-field coupled from sources directly located on the 1DPC surface, it can be detrimental in cases of prism-coupling from external sources, because very little power is resonantly transferred from the glass substrate to the 1DPC top surface. Fostered by an increasing interest toward surface waves on dielectric photonic platforms, we investigated some aspects of the transverse confinement of BSWs in open Fabry– Perot microcavities. Despite the low effective refractive index contrast produced by the ultra-shallow patterning/modulations of the multilayer surface, BSWs can be resonantly confined within sub-micrometric regions. As an extension of the findings presented here, an omnidirectional confinement on the multilayer plane can be obtained by introducing axis-symmetric corrugations [40] or dielectric ridges [46]. The design of the surface modulation is addressed in conjunction with the multilayer design, as dielectric loading/unloading effects can significantly influence the BSW dispersion, resulting in momentum mismatch with the pattern. In the exemplary cases presented here, we designed an open Fabry–Perot cavity by employing materials and geometries that can be realistically handled in cleanroom fabrication processes. The cavity mode follows a parabolic energy dispersion as a function of the transverse in-plane momentum ∝ *β<sup>y</sup>* , thus suggesting a two-dimensional analogy to the case of 3D stacked planar microcavities. Cavity Q factors generally increase with the number of layers within the 1DPC due to reduced leakage losses. Although this feature may look advantageous when a BSW is near-field coupled from sources directly located on the 1DPC surface, it can be detrimental in cases of prism-coupling from external sources, because very little power is resonantly transferred from the glass substrate to the 1DPC top surface.

One of the main advantages in BSW-based platforms is the possibility of enhancing light–matter interaction on the surface of the photonic structure. This aspect is particularly useful in the framework of 2D materials such as TMD, that are often transferred on planar One of the main advantages in BSW-based platforms is the possibility of enhancing light–matter interaction on the surface of the photonic structure. This aspect is particularly useful in the framework of 2D materials such as TMD, that are often transferred on planar or structured surfaces upon exfoliation from bulk. Fostered by the potential applications of polariton handling at room temperature, we addressed the hybridization of a BSW cavity mode with the A exciton in a WS<sup>2</sup> monolayer deposited thereto. Strong coupling is found with both the cavity and the band-edge modes, depending on the spatial overlap with the interacting modal volumes, resulting in vacuum Rabi splitting in the range from 30 meV to 40 meV.

Our findings contribute to the scientific understanding of BSW-based platforms, particularly of polariton control at room temperature. We anticipate future opportunities for embedding TMDs on dielectric multilayers from the perspective of integrated quantum nanophotonic devices [62].

**Author Contributions:** N.M. has set up the computational model and run the calculations; T.-L.G. has performed deposition of TiO<sup>2</sup> and SiO<sup>2</sup> layers, S.P. and M.R. have conducted ellipsometric measurements and fitting procedures for the refractive index estimation; T.G. has provided the BSW cavity design; E.D. has coordinated the work. All authors have contributed to writing the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially funded by the "Dipartimento di Eccellenza 2018–2022" program from the Italian Ministry of Education, University and Research.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** Computational resources were provided by HPC@POLITO, a project of Academic Computing within the Department of Control and Computer Engineering at the Politecnico di Torino http://hpc.polito.it (accessed on 21 February 2023).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
