Germanium Monosulfide as a Natural Platform for Highly Anisotropic THz Polaritons

Terahertz (THz) electromagnetic radiation is key to access collective excitations such as magnons (spins), plasmons (electrons), or phonons (atomic vibrations), thus bridging topics between optics and solid-state physics. Confinement of THz light to the nanometer length scale is desirable for local probing of such excitations in low-dimensional systems, thereby circumventing the large footprint and inherently low spectral power density of far-field THz radiation. For that purpose, phonon polaritons (PhPs) in anisotropic van der Waals (vdW) materials have recently emerged as a promising platform for THz nanooptics. Hence, there is a demand for the exploration of materials that feature not only THz PhPs at different spectral regimes but also host anisotropic (directional) electrical, thermoelectric, and vibronic properties. To that end, we introduce here the semiconducting vdW-material alpha-germanium(II) sulfide (GeS) as an intriguing candidate. By employing THz nanospectroscopy supported by theoretical analysis, we provide a thorough characterization of the different in-plane hyperbolic and elliptical PhP modes in GeS. We find not only PhPs with long lifetimes (τ > 2 ps) and excellent THz light confinement (λ0/λ > 45) but also an intrinsic, phonon-induced anomalous dispersion as well as signatures of naturally occurring, substrate-mediated PhP canalization within a single GeS slab.

by theoretical analysis, we provide a thorough characterization of the different in-plane hyperbolic and elliptical PhP modes in GeS.We find not only PhPs with long life times ( > 2 ps) and excellent THz light confinement ( 0 /> 45), but also an intrinsic, phononinduced anomalous dispersion as well as signatures of naturally occurring PhP canalization within one single GeS slab.
Polaritons are electromagnetic waves formed by light strongly coupled to collective excitations in matter. 1 The hybrid light-matter nature of polaritons offers a promising platform for the manipulation of the flow of light at the nanoscale. 2Notably, phonon polaritons in layered vdW materials such as hBN, -MoO3 or -V2O5 have recently attracted great interest [3][4][5] since, apart from featuring field confinement to the nanoscale, they naturally exhibit anisotropic (and particularly directional) propagation, ultra-long life times (of several ps), and low group velocities. 6Polaritons hold great promises in manifold applications, such as for instance in nanolasers, 7,8 polarization-sensitive detectors, 9 molecular sensors, 10 quantum nanophotonics, 11,12 hyper-lensing, 13,14 waveguiding 15 and nano-optoelectronics. 16 However, a major challenge on the road to such applications is presented by the PhPs exclusively residing in the polar material's reststrahlen bands (RB): within the spectral region between the transverse optical (TO) and longitudinal optical (LO) phonon modes, typically located in the mid-infrared (MIR) to THz part of the electromagnetic spectrum, where the negative sign of the permittivity enables the excitation of polariton modes. 17,18Thus, routes for spectral tunability (e.g.ion intercalation, 5 nano-structuring, 19 isotopic enrichment, 20 carrier photoinjection, 21 or dielectric environment 22,23 ) as well as novel materials with RBs covering complementary spectral bands are of great need.9][30][31] Furthermore, in contrast to in-plane hyperbolic PhPs in meta-materials, PhPs in natural vdW crystals exhibit significantly lower losses and are not limited by the precision of the fabrication. 6Yet, the palette of vdW materials supporting nanoscale-confined PhPs in the THz spectral range is still very scarce.Here, we add a new member to this palette by introducing the family of group-IV monochalcogenides semiconductor compounds MX (M = Ge, Sn; X = S, Se) as a rich platform for THz nanophotonics.Their layered orthorhombic crystal structure -similar to black phosphorous 32 -gives rise to a large anisotropy of their optical, vibrational, and electrical properties. 33All members of the material family show photovoltaic properties, 34 an enhanced thermoelectric effect, 34 large spin-orbit splitting, 35 high ionic dielectric screening, 36 and optical phonons well in the THz spectral range. 33In this study, we focus on the recently predicted 37 THz PhPs in the compound alpha-germanium(II) sulfide (-GeS, GeS) that exhibit an intriguing dispersion in the  = 6.0 -9.5 THz frequency range.Moreover, GeS has a direct bandgap of 1.6 eV, thus promising electric gating as a feasible approach for polaritonic control, 32,38 shows characteristic photoluminescence, 39 an outstanding Seebeck coefficient possibly allowing for photocurrent-nanoscopy 40 of PhPs, 41 ferroelectricity in twisted nanowires 42,43 and in the monolayer limit, 44 resistance to oxidation 34 and fascinating exciton polaritons at visible wavelengths. 45he objective of the present work is the comprehensive characterization of the rich THz PhP modes including their dispersion, volume field distribution, quality factors, life times, and electromagnetic field confinement at the nanometer length scale.To that end, we carry out polariton interferometry experiments by employing a free-electron laser (FEL) as a narrowband THz light source.Our results, supported by full-wave simulations and transfer matrix and analytical calculations, unveil THz PhPs with large quality factors (Q = 10), long life times ( > 2 ps), and deep subwavelength confinement (up to 0/45, where 0 is the incident free-space light wavelength), together with signatures of a natural PhP canalization regime.
The layered orthorhombic crystal structure of GeS (space group Pcmn) is depicted in Figure 1a.In analogy to that of black phosphorous, 32 it consists of covalently-bound layers stacked in [001]-direction with an armchair structure in the [100]-direction and zigzag structure in the [010]-direction. 33The difference in lattice constants (a = 4.29 Å, b = 3.64 Å, c = 10.42Å) is remarkable, in particular the large ratio a/b = 1.18 gives rise to a high structural anisotropy within the layers that is 2.5-times higher than in the likewise biaxial (x() ≠ y() ≠ z()) -MoO3 crystal (a/c = 1.072) 46 in which THz PhPs have been recently demonstrated. 24Our micro-Raman spectrum (Figure 1b) unveils four Raman peaks at Raman shifts of 112, 213, 240, and 270 cm −1 that can be readily attributed to the Ag 3 , B3g, Ag 1 , and Ag 2 phonon modes, respectively. 39Particularly, their polarization-dependence allows deducing the GeS crystal structure orientation of individual flakes.To that end, we find the maximum Raman intensity of the Ag 3 mode in Figure 1c (violet; that is parallel to the [100] crystal axis 47 ) to be aligned along the right edge of the GeS flake marked in the optical microscopy image in Figure 1d.The crystal is composed of layers of covalently bounded Ge (blue) and S (red) atoms with the vdW stacking direction along the [001] crystal direction.Similar to black phosphorous the layers show an armchair and zigzag geometry along the [100] and [010] crystal directions, respectively.The box marks the unit cell.b) Representative Raman spectrum for an incident linear polarization rotated by 4° relative to the [100] crystal axis.The characteristic polarization dependence of the Ag 3 mode at 122 cm -1 (violet) and B3g mode at 213 cm -1 (orange) can be used for unambiguous identification of the crystal axis orientation.c) Polar plot of the normalized Raman scattering intensities of the Ag 3 and B3g mode in (b) with  = 0° corresponding to the [100] crystal direction.The lobes of the two-fold rotational symmetric Ag 3 mode extend along the [100] crystal direction. 47The four-fold symmetric B3g mode is phase-shifted by 45° relative to the Ag 3 mode.d) Optical microscopy image of exfoliated GeS crystals on silicon.The dashed box marks the flake investigated in this work.e) Real (solid lines) and imaginary (dashed lines) parts of the complex permittivity .The permittivity in the THz regime is governed by four optical phonons and exhibits two (overlapping) in-plane reststrahlen bands RBy and RBx (shaded red and blue, respectively).The inset highlights the real part of  from 8 to 9 THz.
In addition to the Raman-active phonons, polar GeS exhibits several well-characterized, directional transversal optical (TO) phonons located in the THz spectral regime 33 that govern its dielectric permittivity  (Figure 1e).We define the coordinate system to align with the GeS crystallographic axes as x = [100], y = [010] and z = [001].At frequencies from 6 to 10 THz, a negative permittivity (Re(i) < 0, i = x, y, z) is found along different crystal axes within four RBs with two of them defined in the plane: RBy (TO,y = 6.06 and LO,y = 9.47 THz) and RBx (TO,x = 7.74 and LO,x = 9.65 THz) shaded red and blue, respectively (note the large overlap as marked by the hatching).In z-direction, GeS exhibits two out-of-plane TO phonons (TO,z,1 = 7.1 THz and (TO,z,2 = 8.4 THz) that spectrally overlap with the in-plane TO phonons, giving rise to an exotic, highly anisotropic optical response.
To experimentally study the excitation of PhPs in GeS within these RBs, we performed s-SNOM-based polariton interferometry applying a narrowband, tunable FEL. 24The experimental setup is sketched in Figure 2a: the pulsed THz radiation produced by the FEL (repetition rate 13 MHz, pulse duration > 5 ps) is focused on a metallized atomic force microscopy (AFM) tip that acts as a nanoantenna providing high k-vectors along with an enhanced, localized electric field.The polarized tip on top of the GeS flake launches PhPs that propagate away from the tip and are back reflected at edges of the 224 nm-thick flake.The electric field of the back-travelling PhPs is scattered by the same tip into the far-field, where it can be detected (see Supporting Information for details on the setup).By raster scanning the tip across the flake at a fixed incident frequency we get a spatial near-field (NF) S2  image of the polaritons' interference pattern. 48To ensure that our near-field images are recorded in an area with a homogeneous flake thickness and sharp edges, we restrict our s-SNOM measurements to the front-facing 90°-corner in Figure 2a (equivalent to the bottom right corner in Figure 1d).In order to corroborate the existence of propagating (Re[k] > Im[k]) PhPs in GeS, we carried out full-wave electromagnetic simulations at representative excitation frequencies within the in-plane reststrahlen bands: RBx and RBy.More specifically, we simulate the electromagnetic fields generated by a vertical point dipole, in analogy to the AFM tip.The result for a frequency of  = 7.33 THz within RBy is shown in Figure 2b (color plot), confirming an exotic polaritonic field distribution: starting from the exciting dipole located in the center of Figure 2b, a polariton propagates along the y (=[010])-direction with hyperbolic wave-fronts.Notably, no PhPs with wavevectors along the x (=[100])-direction are allowed, while the PhPs propagating along the y-direction have momenta  , 7.33 = (2.42 + 0.26) × 10 4 cm −1 ( , 7.33 = 2.60 µm) .The latter are more than 15-times higher than the free-space momentum  0 ≈ 1.54 × 10 3 cm −1 ( 0 = 40.9µm) , thus implying a considerable THz light confinement.The hyperbolic character of the excited PhPs is further supported by their analytically-calculated isofrequency curves (IFC, a section of the dispersion surface for a constant frequency) overlaid with the simulation in Figure 2b that exhibit the shape of open hyperbolas with its major axis aligned along the y-direction.Note that the opening angle of PhPs propagation in real space is well reproduced by the direction of the group velocity (vg, dashed arrow in Figure 2b) in momentum space, which is oriented perpendicular to the IFCs (compare dashed line parallel to vg).The IFCs were obtained using the recently derived equation for the polariton in-plane wavevector  2 =  ∥ 2 =   2 +   2 in a biaxial slab 15,49 with the slab thickness d, the permittivity of the superstrate (substrate) 1 (3), the mode quantization index l, the GeS permittivity tensor diagonal elements x, y, z, the angle between k and the x-axis , and with  = √  (  cos 2  +    2 ) ⁄ .For a frequency of  = 7.33 THz, this analytical model gives a polariton momentum in [010]-direction of  , 7.33 = (2.3 + 0.23) × 10 4 cm −1 ( , 7.33 = 2.73 µm) that excellently matches our numerical estimate.
The simulated field distribution Re(Ez) for a second frequency  = 8.57THz that is located within the overlap between GeS in-plane reststrahlen bands RBx and RBy, is displayed in Figure 2c.Here, the real parts of the corresponding permittivities Re(x) and Re(y) are both negative and, moreover, show a high degree of GeS crystal anisotropy Re(x)/Re(y) = 3.4.Consistently, the spatial distribution of Re(Ez) reveals an elliptically propagating polariton with largely different wave-vectors along the in-plane crystal directions.The simulation predicts the [100]crystal axis to be the prominent propagation direction with  , 8.57 = (3.31+ 0.93) × 10 4  −1 ( , 8.57 = 1.90 µm), accompanied by propagation along the [010]-direction with a about twofold higher momentum  , 8.57 = (6.49+ 4.1) × 10 4  −1 ( , 8.57 = 0.97 µm), albeit higher damping.These numerical findings are supported by the analytical IFCs in Figure 2c that are propeller-shaped with high (low) momenta along y (x)-direction.From Equation (1) we obtain the complex in-plane momenta  , 8.57 = (3.0+ 1.1) × 10  2c strongly imply a natural canalization along the [100]-direction at  = 8.579][30][31] A more detailed discussion of the natural canalization effect will be given below in the context of the dispersion as well as in the Supporting Information.
The optical near-field S2  image recorded at an illuminating frequency  = 7.33 THz (Figure 2d) shows polaritonic features matching the predicted in-plane hyperbolicity.The flake's nearfield response containing the polariton's Re(Ez) features characteristic fringes with a periodicity of half the polariton wavelength parallel to the horizontal flake edge that are caused by PhP propagation along [010]-direction. 4In contrast, similar fringes parallel to the vertical edge are absent, thus confirming the in-plane hyperbolic character of the PhPs at this frequency.In large contrast, at  = 8.57THz the recorded near-field image in Figure 2e shows polaritonic fringes parallel to both flake edges.In particular, the fringe spacing (0.5 p) parallel to the vertical edge is considerably larger than parallel to the horizontal edge.Moreover, while up to three distinct fringes are clearly visible decaying along the [100]-direction, the fringes decaying along the [010]-direction vanish quickly with distance from the flake edge with only a few fringes remaining barely visible.These observations are consistent with the corresponding complex momenta estimated from theory. ) and e), respectively.The black solid lines are fits to the profile using Eq. ( 2) where applicable.
Finally, to determine the experimental polariton momentum along the [100] and [010] crystal directions from the two near-field images, we extract line profiles along the marked positions starting from the flake edges (blue and red, respectively in Figures 2d,e).The width of the marked profiles represents the area used for averaging, which is necessary to suppress the high noise level due to the high-frequency power fluctuations of the THz source.Additionally, the profiles are normalized and off-set for better visibility.The red profile along the [010]-direction in Figure 2f at  = 7.33 THz clearly shows that PhPs decay exponentially with the distance to the GeS flake edge.This fitting (black curve in Figures 2d,e) was carried out using the common analytical description of an edge-reflected polaritonic wave 50 with arbitrary amplitude A and phase x0.The result of the fit yields the complex momentum  , . = (. + .) ×    − ( , . = . µ) that matches excellently to the values obtained from the simulation and analytical model.In the [100]-direction, no clear polaritonic features are visible, which is in line with the theoretical predictions, as well.The bright edge that appears in the corresponding NF image in Figure 2d is related to an edge effect due to far-field scattering off the GeS corner towards the substrate and hence does not appear in the profile for  > .
For the elliptical polariton at  = 8.57THz, both profiles in Figure 2g clearly show features of an exponentially decaying PhP electric field that can be fitted using Equation (2).We retrieve the momenta along the [100]-and [010]-direction to  , . = (. + .) ×    − ( , . = . µ) and  , . = (. + .) ×    − ( , . = . µ) , respectively.Overall, we find a good agreement between the experimental data and the analytical model for ( , ) along both directions.On the other hand, the experimental values for ( , ) fall short of the theoretical values consistently by a factor of 3-4.Generally, the experimental error in Im(k) is considerably higher than Re(k) (here, due to the small amount of fitted oscillations and potential edge-launched contributions), but on average the experimental Im(k) is expected to be larger than the theoretical values, for example due to the decreased phonon life times (and coherence) in a real crystal (related to phonon scattering processes).
By recording near-field images at various illuminating frequencies and fitting the extracted S2  profiles using Equation (2), we obtain the PhP dispersion (k) directly from the experimental data in the frequency range  = 6.0 -8.7 THz (dots in Figure 3a).The left (right) panel relates to PhPs propagating along the [100] ([010])-direction, starting at the TO frequency TO, where the permittivity becomes negative along the respective in-plane direction.Black curves show the corresponding PhP wavevectors Re[k()] calculated using Equation (1), in excellent agreement with the experiment (a minor adaptation to the GeS permittivity is carried out in RBx, see Supporting Information).Moreover, we evaluated the reflectivity rp(,k) of the layered air/GeS/Si-system via the transfer matrix formalism. 51The Im[rp(,k)] contains both the polariton dispersion and damping, with the positions of the maxima yielding the PhP dispersion and their width (FWHM) directly related to their damping.We find the Im[rp(,k)] (false color plot in Figure 3a) matching excellently the experimental and analytical data, thus further corroborating our results.
Along the [100]-direction we find a phonon polariton branch emerging on the dispersion plot at TO,[100] = 7.74 THz: The momentum Re(kx) increases with frequency  up to (  ) = . ×    − (i.e., with a positive group velocity) with the relatively large width of the peak of Im[rp(,k)] indicating considerable damping.Along the [010]-direction, in the frequency range  = 6.06 -7.9 THz, we comparably observe the polariton momentum Re(ky) to increase with frequency up to (  ) = . ×    − , accompanied by a smaller damping as compared to the polariton branch along the [100]-direction.However, at higher frequencies  > 7.9 THz, the behavior becomes much more intricate due to two consecutive back bending effects (anomalous dispersion) emerging in separate spectral areas.In general, back bending of a polariton's dispersion is well-known and can be due to several physical reasons: Firstly, surface plasmon or phonon polaritons show a similar effect at the spectral location where the area of negative permittivity ends. 52,53In this case, the bending appears near the light line and the polariton branch emerging in the area with Re() > 0 lacks interface confinement.Secondly, polaritons coupling to external excitations (such as phonons of the substrate 54 or nearby molecular resonances 10 ) have been reported to induce a back bending or anti-crossing to the otherwise monotonous polariton dispersion.Finally, in the present case of a highly anisotropic material we find the anomalous dispersion to be induced by the PhPs coupling to intrinsic phonons.Such precondition is also met in -MoO3, where a similar, considerably weaker effect was observed recently: 55 The weak [100]-phonon located at wavenumber TO = 998.7 cm -1 couples to the in-plane elliptical PhP (that is caused by the negative permittivity in vdW-direction).As for the present GeS, here the PhPs coupling to the strong (weak) z-phonon with TO,[001],1 = 7.1 THz (TO,[001],2 = 8.4 THz) causes the back bending around  = 8.5 THz (8 THz) (see Supporting Information).As seen in the reflectivity, the effect is accompanied by a substantial polariton damping that renders it challenging to observe experimentally with our current setup.Nevertheless, we measure the highest momenta (  ) = . ×    − in between the regions that exhibit anomalous dispersion.In a nutshell, the highly anisotropic permittivity of GeS governed by overlapping degenerate optical phonon modes in a narrow spectral regime introduces an exotic in-plane PhP dispersion.The latter features several back-bending effects and three characteristic areas with different polariton modes that will be discussed below.Area A, 7.74 -9.47 THz ( TO,[100] - LO,[010] ): In this spectral range, we find the permittivity Re(i) to be negative along all three crystal directions and, hence, the polariton is constituted by the coupled interface modes 2 (air/GeS and GeS/Si).For a representative frequency  = 8.35 THz, the in-plane IFC (Figure 3b), exhibits a propagating  =  mode (black curve).The shape relates to an in-plane elliptical PhPs in agreement with our experimental findings at  = 8.57THz (Figure 2e).In order to obtain the complete picture of the PhP field in three dimensions, we perform full-wave electromagnetic simulations of the layered system.Figure 3c displays a 3D representation of the vertical electric field component with the top face at a distance of 10 nm above the GeS slab, generated by the vertical electric dipole positioned z = 265 nm above the origin.The presented component Re(Ez) is directly linked to the experiment as it provides a valid numerical description of the signals measured in s-SNOM. 56As anticipated from the GeS permittivity, we find confined fields inside the GeS slab in both the x,z-and y,z-planes that are well confined to the GeS/Si interface and decay with z.Similarly, this interface PhP mode resides in the Si, featuring the same momentum ( ∥ ) as within the slab and decaying with decreasing z (into the substrate).Overall, the polariton along the [010]direction holds a considerably higher momentum Re(ky) > Re(kx) and a substantial damping Im(ky) > Im(kx), which is in agreement with our experimental findings at the frequency  = 8.57THz.Area B, 7.1 -7.74 THz ( TO,[001],1 - TO,[100] ): In this spectral range, GeS exhibits negative permittivity Re(y,z) < 0 along both, the [010] and [001] crystal directions, whilst Re(x) > 0, which leads to the excitation of in-plane hyperbolic polaritons. 49At  = 7.5 THz the IFCs in Figure 3d support the hyperbolic character with the l = 0 mode (black curve) forming an open hyperbola with its major axis oriented along the [010]-direction.Of particular interest is the  =  mode (red curve) that likewise holds a hyperbola, which, however, displays its major axis aligned along the [100] direction.While different shapes of IFCs (and thus different directional properties) for the  =  and  =  modes have been theoretically predicted, 49 so far their experimental evidence is lacking.The simulation in Figure 3e supports the findings from the IFCs: The x,y-plane in air shows a hyperbolic PhP propagating along the [010] crystal direction.The field in the y,z-plane demonstrates the coupling of the interface modes.Within the x,z-plane, however, the field in air and Si is dominated by the non-polaritonic near field of the exciting dipole.Intriguingly, inside the GeS slab we find a high momentum PhP that decays along [100] direction rapidly with lateral distance x from the exciting dipole, while the propagating along this direction is forbidden for the zero order ( = ) PhP slab mode, according to its IFC.We explain this observation by the propagation of the higher order ( = ) PhP slab mode with its IFC shown in Figure 3d by the red curve.Area C, 6.1 -7.1 THz ( TO,[010] - TO,[001],1 ): At these frequencies only the [010]-crystal axis exhibits negative permittivity Re(y) < 0, resulting in a type II hyperbolic optical response with hyperbolic PhP slab modes. 2 At  = 6.8 THz the IFCs (Figure 3f) show the peculiar in-plane hyperbolic shape with the hyperbola's major axis aligned along [100]-direction for both propagating modes ( = , ).Compared to  = 7.5 THz (Area B), we find a smaller Re(kx) as anticipated from the dispersion in Figure 3a.The simulated fields in Figure 3g within the x,yplane similar to that in Area B depict hyperbolic PhPs with low momentum ky propagating in [010]-direction in agreement with the the IFC of the  =  PhP mode (black curve in Figure 3f).The fields in the x,z-plane solely show the non-PhP near fields of the exciting dipole along [100]-direction as PhP propagation in this direction is forbidden.In contrast, the field distribution within the y,z-plane depicts a polariton with small damping propagating along the [010]-direction.The remarkably high momentum  =  mode predicted by Equation (1) cannot be easily observed.The primary difference between the PhPs in areas B and C is presented by the rotation of the hyperbolic  = 1 mode, as the  = 0 mode proves to be unaffected by the change of sign in Re(z).A more detailed analysis of the transition at the frequency TO,[001],1 can be found in the Supporting Information.
Ultimately, the aim of this work is to thoroughly characterize the PhPs in GeS to pave the way towards application in science and technology.To end, we determine the key GeS polaritonic properties (figure of merit, life time  and light confinement ) from our experiment as well as the analytical model and compare them to recent, established PhP-hosting materials.
The quality factor Q = (Re(k))⁄(Im(k)) presents a practical figure of merit (FOM) that (in real space) relates the polariton's wavelength to its decay length. 6For GeS, in Figure 4a we find FOMs of up to Q = 10 in RBy and Q = 3 in RBx along the [010] and [100] directions, respectively.The black curves are obtained directly from Equation (1) and describe well the experimental data.The data points distinctly falling below the model can be attributed to an increased bandwidth of the exciting FEL pulse that leads to an artificially increased polariton damping: in fact, an incident pulse with significant bandwidth according to the dispersion in Figure 3a excites a broad range of polaritons with various k-vectors.The super-position of such waves detected in the far-field exhibits a shortened propagation length due to destructive interference (more detailed explanation given in the Supporting Information).Note that the quality factor drops at the frequencies of the in-plane LO and TO as well as in the regions of the back-bending in the dispersion.Overall, the quality factors resemble those reported for -MoO3 (Q ≈ 7 -12) and are 2-times smaller than in naturally abundant hBN (Q ≈ 20) 20 and about 3-times higher than in -V2O5 (Q = 2.5). 5 The GeS PhPs life time ( = [  Im()] −1 with the group velocity   = 2  ⁄  ,  denoting the wavenumber) lies in the picosecond range (Figure 4b) as anticipated for PhPs. 57ithin RBy, life times of up to [010] = 2.3 ps can be found, while the life times in RBx are considerably smaller with [100] < 1.4 ps.The life times are thus comparable with those reported for PhPs in hBN (< 2 ps), 20 but shorter than in case of PhPs in -MoO3 (2 -8 ps) 4 and -V2O5 (3 -6 ps). 5 Note that in the same way as for the FOM, a higher excitation bandwidth can artificially decrease the experimental life time (see Supporting Information).Moreover, the large errors in the determination of life time as well as Q factor are related to the poor signal to noise ratio within the experiment, which is appositely backed up by our simulations: For GeS, we find a smaller overall polaritonic field Re(Ez) per exciting field strength as compared to -MoO3 and hBN, for example.Moreover, it is worth mentioning that following the common definition of the life time (propagation length L = Im(k) -1 divided by the group velocity vg) can erroneously lead to negative values (in the regions of anomalous dispersion, where  ⁄  < 0) as for instance in the dispersion curve along [010]-direction in Figure 3a.In the anomalous dispersion region, one has to use a different (more general) determination of the life time, based on the eigenmode analysis in the space of complex frequency and real wavevector. 10inally, we calculate the thickness-dependent light confinement  =  ⁄  0 =  0 ⁄  (that is, the ratio of the polariton in-plane momentum k in respect to the incident, free-space light k0) at different frequencies in RBy.As presented in Figure 4c, the experimental values of  follow very well the ~1/d-dependence anticipated from Equation (1) (solid curves).The squared data point is taken from the dispersion of the 224 nm-thick flake in Figure 3a at a slightly shifted frequency of  = 7.33 THz.We find in our experiment the highest field confinement of  = 47 in the 224 nm-thick flake at  = 8.57THz, whereas considerably larger values are expected for thinner GeS flakes.
Eventually, we want to address the apparent PhP canalization indicated by the numerical simulation and analytical model presented in Figure 2e.Similar effects have been studied mainly theoretically in the context of plasmonic and phononic metamaterials 19,58,59 and, to our knowledge, have not been experimentally demonstrated in a single layer of a natural material yet.In GeS, the canalization is caused by the anisotropic permittivity that leads to the propellershaped, almost flat IFCs together with a large damping along the [010]-direction, resulting in propagating polaritons with nearly parallel group velocities.Hence, the wave fronts lose their hyperbolic character as observed in the Re(Ez) image.Direct experimental observation of this (nearly) diffraction-less propagation will require an antenna for efficiently launching these polaritons in order to probe this electric field distribution emitted by an effective point source.Hence, a dedicated, frequency-dependent study of this canalization effect will be presented in a future work.
In conclusion, we extensively explored the properties of THz phonon polaritons in the highly anisotropic semiconducting vdW material -GeS, finding different confined polaritonic modes between frequencies of 6 and 9 THz.Of particular interest are the anomalous dispersion and the anticipated natural canalization effect that both originate from the interplay of the highly directional RBs.For these reasons, GeS promises to become a feasible, versatile platform for THz light confinement in the future.Moreover, we believe the work presented here will inspire novel research on THz PhPs: While on one hand, the material family presents a toolbox for THz PhP engineering (for example via stacking and twisting), on the other hand, GeS holds the promise of tuning the PhPs via electrical gating (i.e.PhP-electron interaction).Direct control of the charge carrier concentration moreover enables the study of plasmon-phonon coupling with the goal to actively control the anisotropic polaritons.Lastly, the large thermoelectric effect motivates investigation of the thermoelectric properties of the PhPs that could potentially be probed via photocurrent-nanoscopy 26 .To that end, the scarcity of tunable THz sources currently presents the only limitation.

Scattering-type scanning near-field optical microscopy (s-SNOM):
The near-field images were recorded applying a (modified) commercial near-field microscope (Neaspec GmbH, Germany) conjoined with the free-electron laser located at Helmholtz-Zentrum Dresden-Rossendorf (HZDR), Germany.By illuminating the oscillating (Ω = 250 kHz), metallized s-SNOM-tip in the vicinity of the sample surface, the excited tip acts as an antenna providing a strong, localized electric field at its apex.The confined field interacts with the sample volume and hence, its local optical response becomes imprinted in the back-scattered signal S. In order to separate the sample's near-field optical response from the dominant farfield background, the non-linear distance-dependence of the NF contribution (as compared to the linear dependence of the far-field) is exploited: 60 Through employing a lock-in amplifier we obtain individual components of the scattered signal at multiples of the cantilever oscillation frequency SnΩ (n = 1, 2, 3…) and find effective background suppression in the components with n ≥ 2. 61 Throughout this work, a self-homodyne detection scheme was applied and the backscattered optical signal was demodulated at n = 2.The optical signal at the frequencies  = 6 -9 THz was recorded using a gallium-doped germanium photoconductive detector by QMC Instruments Ltd., UK.

Free-electron laser:
Light sources in the THz spectral regime that are suitable for s-SNOM application currently present a major limitation to near-field optical investigation of collective excitations in condensed matter physics.Established table-top solutions such as gas-lasers or quantum cascade lasers are either restricted by the small range of accessible frequencies or the lack of sufficient spectral power density. 62In contrast, light emission of relativistic electrons can be exploited in large-scale facilities (namely synchrotrons or free electron lasers) to provide either broadband (in the first case) or continuously tuneable, narrowband THz radiation (in the latter case).While synchrotron infrared nanospectroscopy currently is operational at frequencies down to > 9.6 THz, 25,63 FELs in particular have been successfully applied in s-SNOM from 1.3 -30 THz. 24,64 this work, we apply the free-electron laser FELBE at the ELBE Center for High Power

Figure 1 .
Figure 1.Material properties of germanium sulfide (-GeS).a) Crystal structure of -GeS.The crystal is composed of layers of covalently bounded Ge (blue) and S (red) atoms with the vdW stacking direction along the [001] crystal direction.Similar to black phosphorous the layers show an armchair and zigzag geometry along the [100] and [010] crystal directions, respectively.The box marks the unit cell.b) Representative Raman spectrum for an incident linear polarization rotated by 4° relative to the [100] crystal axis.The characteristic polarization dependence of the Ag 3 mode at 122 cm -1 (violet) and B3g mode at 213 cm -1 (orange) can be used for unambiguous identification of the crystal axis orientation.c) Polar plot of the normalized Raman scattering intensities of the Ag 3 and B3g mode in (b) with  = 0° corresponding to the [100] crystal direction.The lobes of the two-fold rotational symmetric Ag 3 mode extend along the [100] crystal direction.47The four-fold symmetric B3g mode is phase-shifted by 45° relative to the Ag 3 mode.d) Optical microscopy image of exfoliated GeS crystals on silicon.The dashed box marks the flake investigated in this work.e) Real (solid lines) and imaginary (dashed lines) parts of the complex permittivity .The permittivity in the THz regime is governed by four optical phonons and exhibits two (overlapping) in-plane reststrahlen bands RBy and RBx (shaded red and blue, respectively).The inset highlights the real part of  from 8 to 9 THz.

Figure 2 .
Figure 2. Near-field imaging of THz PhPs in GeS.a) Schematic of the experimental setup.The AFM tip excited by the FEL's THz radiation launches PhPs that propagate across a 224 nm-thick GeS slab.b),c) Simulated real-space field distributions Re(Ez) at the GeS/air interface (false color plot, scale bar) overlaid with the isofrequency curves in momentum space (solid

Figure 3 .
Figure 3. Dispersion and field distributions of PhPs in GeS.a) Dispersion (k) along the [100] (left side) and [010] (right side) crystal directions for a 224 nm-thick GeS slab.The dots represent the experimental data extracted from near-field profiles and the black curve corresponds to eq. (1) with  = 0 and = 0, /2.The false color plot presents the imaginary part of the reflection coefficient rp(,k) calculated via the transfer-matrix formalism. 51b),d),f) IFCs calculated for three characteristic frequencies  = 8.35, 7.5, and 6.8 THz within the different reststrahlen bands, respectively.The black (red) curves relate to the propagating  = 0 ( = 1) modes.c),e),g) 3D representation of the simulated field distributions Re(Ez) at the same three frequencies for a 224 nm-thick GeS slab between air and an Si substrate.The top face relates to the x,y-plane as observed in the experiment and the right (left) face corresponds to the y,z-(x,z-) face.

Figure 4 :
Figure 4: Characteristic properties of PhPs in GeS.a) Experimental quality factors as a function of frequency for polaritons propagating along [100] (blue) and [010] (red) crystal directions.The two black curves are the analytical estimates obtained from eq. (1) for the inplane crystal axes (= 0, /2 and l = 0).b) Polariton life time as a function of frequency along the in-plane crystal axes.Again, the black curves are obtained from Equation (1).c) Thicknessdependence of the polariton field confinement at four different frequencies within RBy.The experimental data taken at  = 7.32 THz for a set of different GeS flakes follow the well-known ~1/d behavior in eq.(1) (solid lines).The squared data point was taken at  = 7.33 THz on the d = 224 nm-thick flake.

Radiation
Sources at HZDR, Germany, capable of generating coherent THz radiation over the spectral range of 1.2 -60 THz with a repetition rate of 13 MHz.Particularly the U100 FEL oscillator provides the required brightness to launch and detect PhPs in the 6 -9 THz spectral regime.The spectral bandwidth of individual pulses was minimized by slightly detuning the cavity, resulting in values of about 0.5 -0.9%FWHM and transform-limited pulse durations of > 5 ps.The implied pulse spectral diagnostic was performed applying a Czerny-Turner type scanning grating spectrometer (Princeton Instruments SP-300i).