Prediction of hyperbolic exciton-polaritons in monolayer black phosphorus

Hyperbolic polaritons exhibit large photonic density of states and can be collimated in certain propagation directions. The majority of hyperbolic polaritons are sustained in man-made metamaterials. However, natural-occurring hyperbolic materials also exist. Particularly, natural in-plane hyperbolic polaritons in layered materials have been demonstrated in MoO3 and WTe2, which are based on phonon and plasmon resonances respectively. Here, by determining the anisotropic optical conductivity (dielectric function) through optical spectroscopy, we predict that monolayer black phosphorus naturally hosts hyperbolic exciton-polaritons due to the pronounced in-plane anisotropy and strong exciton resonances. We simultaneously observe a strong and sharp ground state exciton peak and weaker excited states in high quality monolayer samples in the reflection spectrum, which enables us to determine the exciton binding energy of ~452 meV. Our work provides another appealing platform for the in-plane natural hyperbolic polaritons, which is based on excitons rather than phonons or plasmons.

The complex sheet optical conductivity (σ) of monolayer BP can be determined from the reflection contrast spectrum through the following expression 1 : where 0 is the vacuum impedance, is the refractive index of the substrate ( = 1.39 for the PDMS substrate). For few-layer BP (layer number >1), the infrared extinction spectra were obtained in the spectral range from 0.36 eV to 1.36 eV. The relationship between the sheet optical conductivity σ and the extinction spectra reads 1 : Generally, the sheet optical conductivity σ, excluding Drude response of possible free carriers, can be modeled by a superposition of Lorentzian oscillators: Here, is the spectral weight, is the inter-band transition frequency, and is the resonance width of the ℎ transition resonance. However, to ensure the validity of fitting results, we rewrite the Supplementary Equation (3): where is the sheet thickness ( = 0.53 nm for monolayer BP), 0 is the vacuum permittivity. For 2L-5L BP, the thickness is assumed to vary linearly with layer number (M): = 0.53 nm. Moreover, is the permittivity, where denotes the oscillator strength. In the limit of → 0, giving the static dielectric constant.
Combining Supplementary Equation (5) (5) into Supplementary Equation (4). Given the asymmetric line-shape of the 1s resonance state 2, 3 , multi-oscillators were utilized, which gives an almost perfect fitting result, as shown in Supplementary Fig. 1a. The color curves there indicate the contributions from each oscillator.
The corresponding fitted dielectric function is presented in Supplementary Fig. 1b. In our analysis, at the zero frequency ( → 0), the fitting results yield a (0) = 12, which agrees well with the theoretical 4,5 and experimental values determined by the far-infrared interference spectra 6 . This validates our fitting results.
For the polarization along the zigzag direction, we rewrite the ( ): The contributions from higher frequency oscillators beyond our measurement range are all accounted by the background dielectric constant b . Keeping in mind of the anisotropy of the BP, for the response along the zigzag direction, we adjusted the parameter b ( b = 3.9) to achieve (0) = 10 (ref. 7). Using Supplementary Equation (7), Equation (4), and Equation (2)

Supplementary Note Ⅲ. Derivation of the HEPs dispersion in anisotropic BP
In this note, we derive the dispersion relation for the EPs in anisotropic BP, which is placed between vacuum ( vac = 1 ) and a PDMS substrate ( sub = 1.93 ) . We also determine the propagating direction of HEPs from the hyperbola asymptote angles.
We utilize the loss function to determine the polaritons dispersion. The loss function −Im(1/ ), is related to the imaginary part of the inverse of the random phase approximation (RPA) dielectric function [14][15][16][17] . The q-dependent dielectric function has the following form in 2D case for the high symmetry directions: To quantify the basic properties of exciton-polaritons, the quality factor Q = ∆ were calculated 18,19 , where and ∆ are the resonance energy and linewidth of the excitonpolaritons. We plot the loss function values versus frequency at a fixed wave vector in Supplementary Fig. 6. The Q values were then extracted from the peak energy and linewidth ∆ of each curve. At a near zero wave vector, the Q factor can reach as high as 35.52, which is reasonably good. As the wave vector increases, the Q factor decreases. Next, let us consider the characteristics of HEPs that can be supported by a hyperbolic material. The low symmetry of BP allows for two independent in-plane components of the optical conductivities, σ AC ( ) and σ ZZ ( ). We assume that the wave vector of HEP is at an angle with respect to the armchair axis, in which case the wave vector is no longer along the principal axes, but along a certain direction. With the magnitude of = √ AC 2 + ZZ 2 , = cos , ZZ = sin , and the sheet optical conductivity can be rewritten as: σ = σ AC cos 2 + σ ZZ sin 2 ,