1. Introduction

A wide branch of photonics in the mid-infrared range (MIR) is based on the ability to excite surface phonon polaritons (SPhPs) at the interface between a polar media and air or a transparent dielectric [1]. Developing photonic devices and platforms in the MIR is an active field of research motivated by applications, such as detection of traces of biomolecules [2] and harmful/explosive substances [3], passive radiative cooling [4] as well as devices for thermal imaging and for medical diagnostics [5]. However, SPhPs spectral range in polar materials, called Reststrahlen band (RB), vary significantly from one polar medium to another, therefore only offering limited range of operation frequencies. There is a strong need to create controllable and adaptable IR optical materials for frequency, directional and polarization tuning. Recently, van der Waals (vdW) [6] materials received wide attention due to their ability to excite SPhPs [2,3]. Furthermore, subwavelength patterning of the vdW polar dielectrics, allows hybridization and localization of SPhPs into ultra-small volumes with large confinement factor [7,8].

Among newly emergent class of van der Waals (vdW) polar materials, molybdenum trioxide (α-MoO₃) is attracting great attention for its birefringence in visible [9] and MIR [10,11]. It supports SPhPs along the three main crystallographic directions (between 550 and 1000 cm⁻¹) [4] and its strong anisotropy allows in plane and off plane hyperbolicity that has led to demonstration of MIR waveplates [12] and tunable [13] and broadband [14] absorption in far-field. In near-field studies, hyperbolic polaritons have recently been shown in α-MoO₃ flakes where the permittivity is negative along the [100] direction and positive along the [001] direction (from 818 cm⁻¹ to 974 cm⁻¹) [5]. Furthermore, the tuning of SPhPs wavefronts in α-MoO₃ has been experimentally demonstrated in some recent studies using two stacked MoO₃ flakes with a varying twist angle [15]. Excitation of SPhPs waves in twisted flakes of 2D α-MoO₃ has been proposed as an efficient and versatile way to get topological transition of SPhPS from hyperbolic to elliptical in the MIR [16]. Precise control of SPhPs properties as a function of twisting angle between the two flakes was experimentally demonstrated.

Controlling and engineering phonon polariton wavefronts and group velocities paves the way to achieve on-demand functionalities such as super-Planckian emission [17], and biosensing [18] and quantum nanophotonics [19]. Such control of SPhPs in twisted flakes plays an important role for near-field radiative heat transfer [16] where radiative heat transfer can be modulated by adjusting the relative rotation angle between the receiver and the emitter. In general, controlling optical properties and designing functional materials and devices in twisted layers is the topic of a newly emerging field called twistoptics [20,21] that is expected to create novel materials and exhibit interesting optical behavior.

Although using twisted layers is an exciting research direction, there are several drawbacks that might limit the wide adaption: a) SPhP dispersion can only be tuned passively; b) layer stacking with a desired angle requires precise alignment that is challenging with the existing alignment techniques for flakes with a few tens of nanometers thickness; c) the flakes often have irregular shape preventing a good superposition domain; d) flakes are usually limited to few 100 µm at most, therefore limiting potential use in large area or integrated/multifunctional devices. Therefore, having a similar degree of freedom in controlling SPhP dispersion within a single film, without the need to rotate, align and twist different layers will useful and less demanding for practical applications.

Here, we introduce an effective medium (metamaterial) approach for tuning and controlling the wavefronts of SPhPs in α-MoO₃ single flakes (or films). One-dimensional (1D) metamaterial approach has been already considered to expand the potential applications of hBN and obtain in plane hyperbolicity [22]. Due to its strong anisotropy and its natural hyperboliticy, α-MoO₃ is a better candidate and offer more rich physical phenomena than hBN. Extending the structural control in 2- and 3- dimensions will enable access to novel optical properties and unpredicted functionalities. This goal can be achieved by combining the tailored dispersion of artificial structures (metamaterials) and the peculiar properties of natural hyperbolic materials, as pointed out in [23]. One such example is the hybridization of SPhPs of different nature with non-radiative modes or longitudinal optical phonons aimed to the development of phonon polariton-based electrically pumped midinfrared emitters [24].

In this work, we propose nanostructured polar materials with subwavelength dielectric elements (i.e. air inclusions) arranged in random patterns to control SPhPs. As a proof of concept, we demonstrate precise control of the polariton dispersion through tailoring dielectric constants of an effective medium composed of α-MoO₃ matrix and randomly distributed, spherical air inclusions. We obtain new and almost non-dispersive bulk mode that enables hybridization with the intrinsic polariton modes of the crystal. By choosing the filling factor and controlling the geometrical parameters of the inclusions, the dispersion contours of the SPhPs modes remarkably switch from hyperbolic (open) to elliptical (closed) as a function of frequency and/or of the geometrical parameters. As a consequence, the resulting anisotropic polar metamaterial exhibits optical responses similar to the twisted flakes, including the arising of a mode characterized by collimated and fairly diffractionless propagation similarly to that observed at the topological transition angle (δ≈63°) in Ref.16 and at magic twist angle in Ref. 15.

2. Effective medium approach

In order to model the metamaterial as an effective medium, we developed an extension of the classical homogenization techniques such as Maxwell Garnett theory, in combination with transfer matrix method (TMM) for anisotropic, dissipative [25] layered structures. This numerical method facilitates the modeling of lossy anisotropic materials, to retrieve the emissivity of films composed of oriented, subwavelength, ellipsoidal inclusions, randomly dispersed into a polar matrix [26]. The shape effects of the dielectric inclusions are taken into account by the depolarization factors calculated along the three axes of the ellipsoidal inclusion in the three orthogonal directions, whose values depend on the ratio between ellipsoid axes [27]. Assuming spherical shape of air inclusions, the inclusions content ratio in the effective medium is taken into account through the filling factor parameter ($f$) and the effective dielectric tensor components ${\varepsilon _{eff,j}}$ are:

For such an effective medium, the typical polariton resonances can be shifted and tuned by changing structural parameters such as filling factor and depolarization factors. The dielectric constants of corresponding metamaterial can be tailored to obtain negative values in specific ranges of the IR spectrum, along one or even two crystal directions, giving rise to uniaxial or biaxial hyperbolic metamaterial (HMM) response. It is worthy of note that the effective medium approach, i.e. a random distribution of inclusions into a polar matrix (or vice versa), relies heavily on the single scatterer response. Therefore, we restricted our analyses to inclusions with filling factor values up to 0.3 [22]. Within this limit, percolation coupling effects among inclusions due to contact can be neglected and the homogenized effective medium behavior is still due to the single element response.

Once the effective dielectric permittivity values are calculated along the three main directions, the equivalent permittivity is introduced as a uniform layer into TMM to carry out a set of numerical simulations of complex reflectivity values as a function of the wavelength and the parallel component of the wavevector with respect to the interface. The imaginary part of the complex reflectivity reproduces SPhP’s dispersion behavior [28].

With the introduction of randomly dispersed air voids into α-MoO₃ matrix results in an effective dielectric permittivity with additional resonances [29]. Depending on the filling factor and depolarization factors, this effect can produce a split in the material RBs, creating an additional band with positive permittivity. Overall, in the effective medium dispersion two close RBs appear. The interesting added feature is due to the fact that the split between the two newly created bands is close to the upper limit of the original RB (close to the longitudinal phonon frequency). As an example, in Fig. 1 we plotted the dispersion of dielectric permittivity tensor components of bulk α-MoO₃ (Fig. 1(a)) layer and a metamaterial composed of α-MoO₃ with spherical air inclusions (Fig. 1(b)) with a filling factor of 0.3. We focus our attention on the ɛ_yy and ɛ_xx components of the dielectric permittivity tensor. In particular, the wide RB related to the negative sign of the ɛ_yy in the bulk dispersion (from 545 cm⁻¹ to 851 cm⁻¹) now has an additional artificial resonance that splits the metamaterial dispersion to two RB bands (545 cm⁻¹ to 820 cm⁻¹ and 835 cm⁻¹ to 851 cm⁻¹). Moreover, the second RB displays a new effective ω_TOy frequency (835 cm⁻¹) which is closer to the ω_TOx frequency (820 cm⁻¹).

 

Fig. 1. (a) Schematics of the effective medium under investigation composing of spherical air inclusions dispersed into α-MoO₃ film. (b) Dispersion plot of the real part of ɛ_yy (blue) and ɛ_xx (red) components of the dielectric permittivity tensor components for bulk α-MoO₃ layer and (c) for a metamaterial with spherical air inclusions having a filling factor 0.3. Highlighted frequency bands correspond to the bulk RB (b) and effective medium RBs (c) along the y-axis where ɛ_yy takes negative values.

Download Full Size | PDF

The effect of such band splitting on the SPhP dispersion is to create SPhPs dispersion branches along zx and zy incident planes in the same frequency range with similar effective linear momentum. In order to further support this effect, we evaluated and compared SPhP’s dispersion by calculating the imaginary part of the TM reflectivity for two different structures. First, 100 nm thick layer of α-MoO₃ on top of a glass substrate and then, an effective medium layer of the same thickness obtained by adding spherical subwavelength inclusions into the α-MoO₃ film and a filling factor of 0.3.

A Cartesian coordinate system is defined such that the x and y axes correspond to the [100] and [010] crystal directions of the α-MoO₃ film, while φ is the propagation angle of SPhPs with respect to the x-axis, i.e. with respect to [100] crystal direction, and is defined as positive in the counterclockwise direction (Fig. 2(a-c)). As a benchmark, we first investigated SPhPs in a homogeneous α-MoO₃ layer (Fig. 2(d-f)). We show that for φ = 0° (xz incident plane - Fig. 2(a)) the SPhPs are excited in the x-Reststrahlen band (frequency range from 820 cm⁻¹ to 972 cm⁻¹, Fig. 2(d)). When φ =90° (Fig. 2(b)), SPhPs exist in the y-Reststrahlen band, (between 545 cm⁻¹ and 851 cm⁻¹, Fig. 2(e)). Using an intermediate angle, i.e. for φ =45° both x- and y- surface waves modes can be excited (Fig. 2(f)). The green and the red lines in the plots represent the ω_LOy (851 cm⁻¹) frequency and the ω_TOx (820 cm⁻¹) frequency, respectively. Within this region, both x- and y- surface waves modes are allowed since both ɛ_xx and ɛ_yy have negative values. We note that despite of an overlap between the RB of ɛ_xx and ɛ_yy, for the homogenous film SPhPs in the xz- and in the yz planes are mainly excited in two well-separated frequency bands.

 

Fig. 2. (a-c) Cartesian coordinate system with x and y axes along the [100] and [010] crystal directions of the α-MoO₃ film where φ is the propagation angle of SPhPs with respect to the x-axis. (d-f) SPhPs calculated for a homogeneous α-MoO₃ layer for φ set to (d) 0°, (e) 90° and (f) 45°, respectively. (g-i) SPhPs calculated for a system composed of α-MoO₃ layer and randomly dispersed spherical air-inclusions, with filling factor set to 0.3.

Download Full Size | PDF

3. Hybridization of surface modes

When spherical air-inclusions are added into the polar matrix, the split in the y-RB produces two different “families” of SPhPs (Fig. 2(h)). We note that the upper band family frequency range matches the RB of the xz SPhPs (Fig. 2(g)). This way, by rotating the plane of incidence it is possible to obtain hybrid modes. In particular, it is interesting to note that at the edge of the newly generated y’-RB, there is an almost dispersionless bulk mode corresponding to an effective ω_TOy frequency which is artificially engineered (see Fig. 1(b)). In Figs. 3(a) and 3(b) we show a zoomed-in version of the SPhPs dispersion for the same systems: homogenous layer (Fig. 3(a)) and effective medium (Fig. 3(b)). The coupling between the two modes is evident from the anti-crossing point reported in Fig. 3(b).

 

Fig. 3. (a) Details of Fig. 2(f): SPhPs calculated for φ set to 45° for a homogeneous α-MoO₃ layer and (b) details of Fig. 2(i): SPhPs calculated for φ set to 45° for a system composed of α-MoO₃ layer and randomly dispersed spherical air-inclusions, with filling factor set to 0.3.

Download Full Size | PDF

Different surface modes corresponding to the two branches forming for frequencies below and above 838 cm⁻¹, respectively, display completely different wavefronts and, consequently, propagating modes. Following these considerations, we calculated the isofrequency curves highlighting how surface phonon polaritons propagate inside the frequency range where x- (lattice-originated) and y’- (inclusions-induced) surface modes are overlapping.

In order to check reliability of the effective medium approach and to compare with fully numerical finite difference time domain (Lumerical FDTD 3D solver) simulations, we focused our attention to a 1 µm thick MoO₃ film on a SiO₂ substrate. We chose subwavelength (100 nm radius) spherical inclusions with a filling factor of 0.3 in an extended spatial range (we considered a patterned area of 4µm x 4µm). Thinner films would have required inclusion of smaller spheres and finer discretization meshes and eventually significantly large amount of inclusions to cover the same area. Our choice is a nice trade-off between effective medium validity ranges and computational resources needed. However, thicker films support multiple modes at the same frequency. Indeed, isofrequency curves (evaluated with the TMM method, displayed in Fig. 4) clearly show multiple modes and the corresponding fields spatial patterns (calculated with the FDTD method considering a point dipole excitation and depicted in Fig. 5) are not as clean as they would appear for single surface mode excitation. Nevertheless, it is still possible to see and analyze the fingerprint of mode hybridization.

Starting with the analysis of the TMM results, displayed in Fig. 4, we show that while ω slightly varies between ω_TOx and ω_LOy i.e. from 820 to 850 cm⁻¹, the wavefront of phonon polaritons highlights a topological transition allowed in a narrow frequency range. In particular, in Fig. 4 we show that for ω= 828 cm⁻¹ (Fig. 4(d)) the polaritonic mode launched at ω>ω_TOx, where real(ɛ_xx) is negative (Fig. 1(b)), has a typical hyperbolic behavior along x-direction. However, within this frequency range, the real(ɛ_yy) of the effective medium (metamaterial) goes through a transition from positive to negative values, although it displays an anomalous dispersion due to the newly introduced resonance. As a consequence, a frequency can be found ω≈838 cm⁻¹ (Fig. 4(d)), where the phonon polariton wavefront switches from hyperbolic to oval, appearing as slightly closed fringes (Fig. 4(e)). This frequency corresponds to the minimum of real(ɛ_yy), i.e. where real(ɛ_yy) switches from anomalous to normal dispersion. At 848 cm⁻¹ (Fig. 4(f)) the dispersion switches back to hyperbolicity in x-direction. As a comparison we plot the isofrequency curves calculated at the same frequencies for the unprocessed α-MoO₃ film which shows hyperbolic dispersion along x for all frequency values (Fig. 4(a-c)).

 

Fig. 4. Isofrequency curves calculated for the unprocessed MoO₃ film (a-c) at three different frequencies inside the x-RB. When spherical air inclusions are added into MoO₃ matrix (d-f) the isofrequency curves, calculated at the same frequencies, show that wavefronts remarkably change within the newly introduced RB.

Download Full Size | PDF

From the numerical results summarized in Fig. 4, it is evident that the metamaterial approach allows SPhPs dispersion tuning and it is possible to determine an equivalent topological transition frequency, (ω=838 cm⁻¹) where the directional phonon polaritons, are collimated and are characterized by an almost diffractionless propagation. At this specific frequency at which the topological transition occurs, the dispersion must necessarily flatten, as its topology changes from closed to open fringes, leading to diffractionless and low-loss directional SPhPs canalization.

In Figs. 5(a-f) we report FDTD simulations of the field patterns (E_z, i.e. the normal component of the electric field with respect to the planar interface) evaluated for the three frequencies under analysis. We compare the fields evaluated in a plane placed 30 nm above the MoO₃ surface when a z-oriented dipole source located 100 nm above the surface is used as a source. Figures 5(a-c) illustrate the field profile for a uniform MoO₃ film while the corresponding fields for the case of the porous film are shown in Fig. 5(d-f). We notice that the field profiles for the uniform film do not sensibly vary in the selected frequency range (828–848 cm⁻¹). Indeed, in this range both in plane dielectric permittivities of the MoO₃ are negative, the medium is indefinite and no evidence of hyperbolic surface phonons results from numerical simulations (Fig. 5(a-c)). On the other hand, for the porous material it is possible to notice the formation of a hybrid-hyperbolic surface phonon polariton (Fig. 5(d). As previously mentioned, numerical simulations are strongly affected by scattering losses because accurate modeling of the porous material on large areas require heavy computational resources. However, beside the scattering related noise, it is clearly visible that for higher frequencies with respect to the characteristic frequency (838 cm⁻¹) the energy flow is mostly collimated along the horizontal (x-) direction, as expected form the isofrequency curves of Fig. 4(e-f), while the signature of the hybrid hyperbolic surface polariton becomes evident in Fig. 5(d). The field profile at 828 cm⁻¹ displays the typical hyperbolic behavior with the energy flowing mainly in a “X” pattern according to the isofrequency results depicted in Fig. 4(d).

 

Fig. 5. FDTD simulations of the fields’ profiles evaluated for the three frequencies under investigation (a-c) for the uniform α-MoO₃ film and (d-f) for the porous film, respectively.

Download Full Size | PDF

It is worthy of note that further functionalities could be exploited when addressing the different parameters of the effective medium. Tuning and varying inclusion content along layer thickness allows one to obtain strongly asymmetric emissivity features paving the way for applications such as strongly directional emissivity filters and thermal logical elements [14]. The topological transition can also be induced by changing the inclusions’ dielectric constant, i.e. using a medium whose refractive index might depend on an analyte concentration [30]. This opens an interesting perspective for sensing applications.

4. Conclusions

In conclusion, we propose an effective medium metamaterial approach to control and tune the surface phonon polariton dispersion relation and wavefronts. Our method can be further extended to other types of vdW materials supporting SPhPs. Our results highlight how effective medium theory can be applied to model and predict exotic wavefronts of polaritons in anisotropic materials and demonstrate that nanostructuring van der Waals materials can form a highly variable and compact platform for hyperbolic infrared metasurface devices and circuits. Despite the actual technological restrains and feasibility of porous vdW metamaterials, this simple but versatile approach suggests a path to enable hybridization of polaritons. Different patterning with respect to spherical pores could be explored in the future with the aim to develop large area patterning techniques. This way it could be possible to develop several patterned areas on the same polar material matrix to conceive and implement different functionalities in distinct areas of the same platform [18]. Finally, this approach allows for further tuning of inclusions dielectric constant values and consequently dispersion curve engineering of the surface phonon polaritons for sensing applications.