Analytical approximations for the dispersion of electromagnetic modes in slabs of biaxial crystals

Anisotropic crystals have recently attracted considerable attention because of their ability to support polaritons with a variety of unique properties, such as hyperbolic dispersion, negative phase velocity, or extreme confinement. Particularly, the biaxial crystal $\alpha$-MoO$_3$ has been demonstrated to support phonon polaritons, light coupled to lattice vibrations, with in-plane anisotropic propagation and unusually long lifetime. However, the lack of theoretical studies on electromagnetic modes in biaxial crystal slabs impedes a complete interpretation of the experimental data, as well as an efficient design of nanostructures supporting such highly anisotropic polaritons. Here, we derive the dispersion relation of electromagnetic modes in biaxial slabs surrounded by semi-infinite isotropic dielectric half-spaces with arbitrary dielectric permittivities. Apart from a general dispersion relation, we provide very simple analytical expressions in typical experiments in nano-optics: the limits of short polaritonic wavelength and/or very thin slabs. The results of our study will allow for an in-depth analysis of anisotropic polaritons in novel biaxial van der Waals materials.

Anisotropic crystals have recently attracted considerable attention because of their ability to support polaritons with a variety of unique properties, such as hyperbolic dispersion, negative phase velocity, or extreme confinement. Particularly, the biaxial crystal α-MoO3 has been demonstrated to support phonon polaritons, light coupled to lattice vibrations, with in-plane anisotropic propagation and unusually long lifetime. However, the lack of theoretical studies on electromagnetic modes in biaxial crystal slabs impedes a complete interpretation of the experimental data, as well as an efficient design of nanostructures supporting such highly anisotropic polaritons. Here, we derive the dispersion relation of electromagnetic modes in biaxial slabs surrounded by semi-infinite isotropic dielectric half-spaces with arbitrary dielectric permittivities. Apart from a general dispersion relation, we provide very simple analytical expressions in typical experiments in nano-optics: the limits of short polaritonic wavelength and/or very thin slabs. The results of our study will allow for an in-depth analysis of anisotropic polaritons in novel biaxial van der Waals materials.

I. INTRODUCTION
Anisotropic media have been a subject of fundamental and applied research in optics for several centuries since the earliest Bartholinus studies. Particularly, birefringence, responsible for the double refraction inside anisotropic crystals, is widely used nowadays in daily-life applications requiring polarization filtering as, for instance, sun glasses, liquid crystal displays, or scanning laser polarimetry (for monitoring glaucoma) [1]. In recent decades, the scientific interest to anisotropic optical phenomena has dramatically increased due to the design and fabrication of novel artificial materials (metamaterials) with a tailored optical response. Striking examples of the latter are photonic and plasmonic crystals [2][3][4] and metasurfaces [5,6], showing spectacular phenomena such as negative refraction [7], slow light [8], and superlensing [9,10], among others. Apart from these anisotropic artificial materials, a few years ago the concept of "atomic-scale" engineering with naturally anisotropic van der Waals (vdW) materials [11] was suggested, adding more scientific interest to this field. Currently, presenting one of the main strategies in low-dimensional optoelectronics, this concept has induced an intensive study of highly-confined anisotropic polaritons supported by vdW slabs and heterostructures [12,13]. The possibility of visualizing these polaritons in thin slabs of vdW crystals with the use of near-field microscopy [14][15][16][17][18] stimulates more and more experimental and theoretical studies in this direction.
From a theoretical point of view, in bulk uniaxial crystals (such as h-BN, SiC or layered metamaterials), characterized by two refractive indices, electromagnetic eigenmodes present ordinary and extraordinary waves. In many cases, the propagation of light along the boundaries of uniaxial crystals and inside the slabs can be straightforwardly analyzed analytically [14,19,20]. In stark contrast, biaxial crystals (such as α-MoO 3 or V 2 O 5 ) are characterized by three refractive indices and both electromagnetic eigenmodes are extraordinary. As a result, the understanding and analytical treatment of electromagnetic phenomena in biaxial media is significantly more complex than in the uniaxial case. In this context, a very recent study has reported on a rigorous analytical solution for the dispersion of surface waves on the boundaries of biaxial crystals [21]; however, up to now, studies on the electromagnetic modes in biaxial slabs have been mainly the subject of a numerical analysis [22][23][24] or some particular configurations, such as grounded crystal slabs with modes fixed propagation directions [25].
In this work, organized in a tutorial style, we present a detailed derivation of the dispersion relation of electromagnetic modes in a biaxial slab of a finite thickness (with arbitrary dielectric tensor), surrounded by two semi-infinite isotropic media with arbitrary dielectric permittivities. We assume that one of the principle crystal axes is perpendicular to the faces of the slab, while the mode propagates along the crystal slab at an arbitrary angle with respect to the other principal axes, which lays in a plane parallel to the faces of the slab. We show that our general dispersion relation successfully reduces to the known limiting cases, such as the case of a uniaxial slab or a semi-infinite crystal, among others. We manage to reduce the general dispersion relation to simple analytical expressions for short wavelength of the modes and small slab thicknesses, which are currently of great interest for the study of anisotropic polaritons in vdW slabs. To demonstrate the validity of our analytical approximations, we compare them to full-wave simulations, finding an excellent agreement.

II. INFINITE BIAXIAL CRYSTAL
Let us consider an infinite, nonmagnetic biaxial medium with dielectric permittivity tensorε. The coordinate system {x, y, z} is chosen is such a way thatε is diagonal [see Fig. 1 To accurately decompose the electromagnetic fields in the biaxial medium, we need to define appropriate basis vectors. To that end, we follow a standard procedure, as for example in Ref. [21]. Namely, we represent the electric and magnetic fields in the biaxial medium in the form of plane waves: where e and h are unknown dimensionless field basis vectors, E 0 and H 0 are arbitrary field amplitude coefficients, ω is the angular frequency, k is the wave vector, and r is the radius vector.
Substituting E from Eq. (2) into Eq. (3), we obtain a linear homogeneous system of equations for the three components of the unknown basis vector, e: where k 0 = ω/c is the free space wave vector, q x,y = k x,y /k 0 are the in-plane components of the normalized wave vector, and q z is the out-of-plane component of the normalized wave vector, so that k z = ±iq z k 0 . The +(−) sign must be taken for the wave propagating along (opposite to) the z-axis, while ∆ i are defined as Our choice of the dependence of the basis vectors upon the coordinates (the propagation along z is treated differently from the propagation along x and y) is dictated by the geometry of the problem (see Fig. 1), consistently with the standard waveguide theory.
The system (4) has nontrivial solutions only when det (M) = 0, that gives the well-known Fresnel's equation for biaxial media [26,27]: It is the quadratic equation in terms of the squared z-component of the wave vector, q z . Its solutions q ez and q oz read as with D being the discriminant In Eq. (7), the sign "+" and "−" correspond to the labels "o", and "e", respectively. Substituting Eq. (7) into the system (4), we find all three components of the two eigenvectors e. Since the system (4) is homogeneous, one of the components of the eigenvectors must be fixed. Without loss of generality, fixing the y-component to e y = q x for the root "o" and to e y = q y for the root "e", we find where the factor 1/q (q being the normalized in-plane wave vector, ) stands for the normalization, and From Eq. (6) we can easily find the asymptotes of the isofrequency curves for large q x,y . Tending both q x and q y to infinity, and setting q z = 0, we find In thin vdW slabs these asymptotes yield the direction of the propagation of the polaritonic "rays", excited by localized sources [15][16][17]. In the particular case of a uniaxial crystal (with the axis C pointing parallel to the z-axis, C Oz), ε x = ε y = ε ⊥ and ε z = ε , the derived basis vectors (9) can be straightforwardly transformed to the basis vectors for the ordinary and extraordinary waves. Taking into account that the z-components of the wave vectors (7) are reduced to the well-known expressions for the ordinary and extraordinary waves we find that ∆ 1 = 0, ∆ 2 = ε ⊥ + q 2 ez and obtain the basis vectors In case of an isotropic medium, ε ⊥ = ε = ε, the z-components of the wave vectors degenerate q 2 oz = q 2 ez = q 2 z = q 2 − ε and the basis vectors (9) reduce to the ones for the s-and p-polarized waves (e o → e s and e e → e p ):

A. General form of the dispersion relation
Let us first represent the electric fields above (z > 0) and below (z < −d) the slab, in the isotropic media "1" and "3", respectively. In these regions we can take the fields in the form of the s-and p-polarized plane waves. For compactness, from now on, we will use Dirac notation, in which the s-and p-polarization basis vectors read where q 1,3z = q 2 x + q 2 y − ε 1,3 > 0 is the out-of-plane component of the normalized wave vector. Here and in the definition of the p-polarization basis vector, |p 1,3 ± , the +(−) sign should be taken for the wave propagating along (opposite to) the z-axis, while in case of the s-polarization |s 1,3 + and |s 1,3 − are degenerated. For convenience, we also introduce the in-plane subvectors of the vectors given by Eq. (15): The fields of the mode propagating along the slab in the upper and lower media can be compactly written as the sum of the s-and p-polarized plane waves: with unknown amplitudes a 1,3 β . In contrast, the electric fields inside the biaxial slab (0 > z > −d) should be represented with the help of the basis vectors found in Section II as where |γ ± denotes |o ± and |e ± , being the polarization basis vector in the biaxial slab The factors a 2↑ γ and a 2↓ γ represent the unknown amplitudes of the plane waves travelling along and opposite to the z-axis, respectively. Analogously to isotropic regions, we can introduce the in-plane subvectors |o and |e of the vectors |o ± , |e ± , respectively. These subvectors can be compactly written as where |u x = (1, 0) T , and To find the magnetic fields we use Maxwell's equation, In case of a plane wave it simplifies to where q i is the normalized wave vector of a corresponding plane wave. Then, we can apply the boundary conditions, which imply the continuity of the in-plane components of both electric and magnetic fields on the faces of the film (at z = 0 and at −d): where the subscript "t" in Eqs. (23) stands for the in-plane subvectors. According to Eq. (22), the in-plane subvectors of the magnetic field can be written as H t = e z × q i × E.
Using the field representation (17), (19), we can rewrite the boundary condition at z = 0 in Eqs. (23) in a more explicit way: where q 1,3± = (q x , q y , ±iq 1,3z ) T and q γ± = (q x , q y , ±iq γz ) T . To simplify Eq. (25), let us introduce auxiliary threedimensional vectors, |β 1,3 ′ ± and |γ ′ ± : where β = s, p and γ = o, e. Calculating the vector products in Eq. (26) we obtain the following explicit relations for the in-plane two-dimensional subvectors |s 1,3 ′ and |p 1,3 ′ : being Y 1,3 β the admittances for the s-and p-polarized waves: Since according to Eq. (26), the z-component of the three-dimensional vectors |γ ′ ± is 0, we keep the same notation for their two-dimensional in-plane subvectors: |o ′ ± and |e ′ ± , which read explicitly as with auxiliary vectors As a result, using definition (26) and Eqs. (27), we obtain a simple form of the Eq. (25): If we multiply (24) and (31) by β| [here, only the exponential of the bra-vector should be complex conjugated, for example, s| = (−q y q x )e −ikxx−iky y ] and taking into account that β|β ′ = δ β,β ′ , we get the following system of equations corresponding to the boundary condition at the interface z = 0: Analogously, for the interface z = −d, with the help of the auxiliary vectors |γ ′ ± we find γ=o,e a 2↓ γ β|γ e qγzk0d + a 2↑ γ β|γ e −qγz k0d − a 3 Equations (32) and (33) form a system of eight linear equations with eight unknowns. By defining ξ γ↓ = e qγzk0d and ξ γ↑ = e −qγz k0d with γ = o, e (for the waves propagating along and opposite to z-axis, respectively), we have: Using the explicit expressions for the vectors |s and |p , |o and |e , and |o ′ ± , |e ′ ± , given by Eqs. (16), (20), and (29), respectively, the scalar products in Eq. (34) can be explicitly calculated as where , and η 4 = 1 + q 2 x c2 q 2 . The homogeneous system (34) has non-trivial solutions only when its determinant equals to zero. The zeros of the determinant yield the dispersion relation for the modes in the biaxial slab. In general, the dispersion relation can be analyzed numerically, but in the limit of a small slab thickness, k 0 d ≪ 1, as well as in the short-wavelength limit (large values of q), it can be written in a compact analytical form, as will be shown below, in Sections V and VI, respectively. Before considering these interesting limits, we will ensure that our dispersion relation analytically reproduces some known examples.

IV. VERY THICK SLABS: SURFACE MODES AT THE BIAXIAL CRYSTAL BOUNDARIES
Let us consider now another extreme case, assuming that the thickness of the slab tends to infinity, d → ∞. Then our dispersion relation should split into two independent dispersion relations describing surface modes at the interfaces between the biaxial crystal and two isotropic media with dielectric permittivities ε 1 and ε 3 . To obtain these dispersion relations in a simple analytical form, we multiply the third and fifth columns of the determinant of the system (34) by ξ o↑ and ξ e↑ , respectively. Then tending d → ∞ in the determinant, and assuming that both q oz and q ez have a non-vanishing real part, we see that all the matrix elements proportional to ξ o,e↑ (third and fifth elements in the four first rows and fourth and sixth ones in the four last rows) vanish. As a result, the 8 × 8 determinant becomes a product of the two determinants 4 × 4, each of them describing the surface modes at the 1-2 (z = 0) and 2-3 (z = −d) interfaces. Without loss of generality, let us consider only one of these determinants 4 × 4, corresponding to the interface 1-2. Zeroing the determinant, we have Then, using the Gauss method, we can reduce the dimension of the matrix to 2 × 2, as To write the dispersion relation in a compact form, we express ε z from Frensel's equation for biaxial slabs (6) as a function of q oz : Then, using the identities q 2 = ε 1 + q 2 1z and q 2 = q 2 x + q 2 y , we substitute ε z from Eq. (44) into the expressions for Y 1 s,p in Eq. (28) and scalar products given by Eq. (35), appearing in the elements of the matrix in Eq. (43). After some algebraic operations, Eq. (43) reproduces the dispersion relation for surface waves on boundaries of biaxial crystals, derived in Ref. [21]: A. Uniaxial crystal with the axis perpendicular to the interface Consider a uniaxial crystal with the axis C along the z-axis, thus directed perpendicularly to the interface of the crystal. Defining, as before, ε x = ε y = ε ⊥ and ε z = ε , and taking into account that according to Eq. (12), ε ⊥ = q 2 − q 2 oz , Eq. (45) simplifies as: Since Re(q oz ) > 0 and Re(q 1z ) > 0, and therefore q oz (q z + q oz ) = 0, we can divide Eq. (46) by q oz (q z + q oz ). Then it transforms to: Assuming that q ez + q oz = 0, we obtain the dispersion relation for the surface wave on a boundary of a uniaxial crystal: where q 1z = q 2 − ε 1 and q ez = ε ⊥ ε q 2 − ε ⊥ . Deriving q from this equation, the dispersion relation takes the well-known form (see e.g. Ref. [28]): In the isotropic case, ε = ε ⊥ = ε, Eq. (48) simplifies to the dispersion relation for surface waves at the interface between two isotropic media:

B. Uniaxial case with in-plane crystal axis
Consider a uniaxial crystal with the C axis along the y-axis, thus lying in the plane of the interface of the crystal. Redefining here ε x = ε z = ε ⊥ and ε y = ε , Eq. (7) for q o,ez transforms to (similar to Eq. (12)): Substituting the expression for q ez from Eq. (51) into Eq. (45) and dividing it by q ez , we get the famous dispersion relation for Dyakonov surface waves [19]:

V. ULTRATHIN SLAB LIMIT
Recently, polaritons in ultra-thin slabs and monolayers (for instance, plasmon polaritons in a monolayer graphene [29] or hyperbolic phonon polaritons in thin slabs of polar dielectrics, such as h-BN [30]) have attracted particularly high attention. Therefore, the limit of a vanishing slab thickness d → 0 is of a great practical interest. Let us illustrate how our general dispersion relation, given by the determinat of the system (34), can be simplified for ultra-thin slabs. Analogously to the methodology used for isotropic slabs [31], we can approximate the slab of a finite thickness by a two-dimensional conductive sheet, with the effective conductivity,σ, given byσ = ωdε 4πi . To that end, let us assume that all the components of the tensorε are large, i.e. |ε i | ≫ 1 (i = x, y, z). Then, retaining in Eq. (7) the first non-vanishing terms depending upon q x and q y in the expressions for the normalized z-components of the wave vectors, q o,ez , Eq. (7) can be greatly simplified: where "+" and "−" should be taken for the labels "o", and "e", respectively. Equations (53) then further simplify to Using Eq. (54), we find c 1 = − εz q 2 y and c 2 = −1, so that the scalar products (35) can be written as Additionally, assuming the small thickness of the slab, we can simplify the elements of the matrix in Eq. (34) by expanding the exponentials ξ γ into the Taylor series in k 0 d and retaining the first non-vanishing terms. We have ξ γ↓ = e qγzk0d = 1 + q γz k 0 d and ξ γ↑ = e −qγzk0d = 1 − q γz k 0 d. To simplify the determinant of the system (34), we sum up its third and fifth columns with the fourth and sixth columns, respectively, and then subtract the fourth and sixth columns (both multiplied by the factor 1/2), from the third and fifth columns, respectively. Then using row operations, we eliminate two first and two last columns in the obtained determinant by the Gauss method: we multiply first (second) row to Y 1 s (Y 1 p ) and subtract it from the third (fourth) row and analogously for fifth (sixth) and seventh (eighth) rows. As a result, we get the following equation: where α x,y = 2πσx,y c = k0dεx,y 2i are the normalized 2D effective conductivity components. Using the smallness of k 0 d, on the one hand and the assumed large values of the components of the tensorε, on the other hand, the determinant (56) can be further significantly simplified. Namely, a more detailed analysis (which we omit here) shows that the elements proportional to k 0 d (the first and third elements of the third and fourth rows) yield the contribution of a higher order of smallness and thus can be neglected. As a result, the determinant (56) factorizes into a product of the two determinants of the sub-matrices 2 × 2: Taking into account that the first determinant in Eq. (57) gives − εz q 2 εx−εy εx−εz = 0, the dispersion relation is given by the vanishing of the the second determinant: This dispersion relation, written for biaxial slabs of a small but nonzero thickness (the effective conductivities α x,y are thickness-dependent), has been used for the analysis of hyperbolic phonon polaritons in thin slabs of α-MoO 3 [16]. Nevertheless, to our knowledge, it has not been consistently derived for a nonvanishing slab thickness up to now. For 2D anisotropic sheets (of zero thickness) Eq. (58) is exact and had been reported in Refs. [32,33]. As expected, the asymptotes of the dispersion relation (58) (q x,y → ∞) coincide with the ones, following from the Fresnel equations [compare with Eq. (11)]: In case of an isotropic 2D sheet α x = α y = α and Eq. (58) splits into two independent equations describing the dispersion of the TE and TM modes [34] propagating along the sheet: To demonstrate the validity of the simplified dispersion relation (58), we compare in Fig. 3 the refractive indices of a mode found from (58) (solid curves) to those found from full-wave simulations (points). Figures 3(a) and 3(b) represent the result for two different illustrative sets of parameters: in Fig. 3(a) only one of the in-plane dielectric permittivities is negative, ε x < 0 (while ε y , ε z > 0), and in Fig. 3(b) both in-plane permittivities are negative, ε x , ε y < 0 (while ε z > 0). Both in Figs. 3(a) and 3(b) the propagation of the mode at different angles ϕ (see Fig.  2) is considered. As expected, the agreement between the analytical approximation and rigorous numeric simulations improves for smaller values of k 0 d, although in Fig. 3(b) the agreement is good in the whole shown range of k 0 d. Impressively, the agreement between the analytical and numerical results is in general excellent for all the shown propagation directions, even though neither k 0 d is very small nor the values of ε x , ε y , and ε z are very large, as it was initially assumed for the derivation of (58).

VI. THE LIMIT OF A LARGE REFRACTIVE INDEX OF THE MODES
The general dispersion relation given by the zeroing of the determinant in Eq. (34) can be also greatly simplified under the assumption of large refractive indices of the modes, |q| ≫ 1. Such simplification is similar to the one made for the dispersion of the modes in uniaxial crystal slabs [14]. For large q the expressions for the z-components of the wave vectors inside the slab can be approximated as: where we have retained the second-order (the highest-order) term in q both in q 2 ez and q 2 oz , as well as the zeroth-order term in the expression for q 2 oz , to avoid uncertainty in ∆ 1 (since ∆ o x = q 2 x and ∆ z = q 2 oz = q 2 , we have ∆ 1 = 0 0 ). Substituting the Eqs. (61) into Eqs. (10), (21), we have where we have neglected all small amendments in the expressions for each constant. Then we obtain the simplified expressions for the scalar products (35): We can verify that Eq. (71) transforms into the dispersion of modes in a uniaxial slab (with the axis C along the zaxis), setting ε x = ε y = ε ⊥ and ε z = ε . Then defining ψ = −ρ = −i ε ε ⊥ , and taking into account that ε1,3ρ ε = ε1,3 ψε ⊥ , we reproduce the dispersion relation, used for the analysis of hyperbolic phonon polaritons in h-BN crystal slabs [14]: The same results can be straightforwardly derived from the exact Eq. (39), in the limit of large q. Equation (71) also reduces to Eq. (72) when the propagation of the mode coincides either with the x axis (in this case we should set ε x = ε ⊥ and ε z = ε ) or with the y-axis (in this case we should set ε y = ε ⊥ and ε z = ε ). In the two latter particular cases, anisotropic polaritons in α-MoO 3 slabs were studied via Eq. (72) in Ref. [18].
(73) Fig. 4 shows the isofrequency curves extracted from both Eq. (71) (red and black curves) and full-wave numeric simulations (dots). For the parametric sets a), c), d) corresponding to the volume modes, the isofrequency curves for the two lowest modes are shown: l = 0 (black discontinuous curve and black dots) and l = 1 (red continuous curve and red dots). In contrast, for the set of parameters b), the mode exponentially decays inside and outside the slab (so that it has a surface wave character) and therefore only the solution with l = 0 makes sense. In all panels of Fig. 4 we see an excellent agreement between the numeric simulations and analytical approximations for large q (q 10) and even very reasonable agreement for q comparable to 1. This agreement, particularly for the case of small and moderate values of q, unambiguously evidences that Eq. (71) can be used in a wide space of parameters for the characterization of diverse modes in both natural and artificial biaxial crystal slabs.

VII. CONCLUSIONS
To summarize, we have presented an analytical derivation of the electromagnetic modes that can be guided along biaxial crystal slabs. We have provided simple expressions for the dispersion of the modes in the limit of an ultra-thin slab and for the case of large k-vectors of the modes. Both limits are currently of great importance for studying highly-confined anisotropic polaritons in vdW biaxial crystal slabs and, particularly, for the interpretation of the state-of-the-art near-field experiments.