Perfectly matched layer for biaxial hyperbolic materials

Zixun Ge; Zixun Ge; Zixun Ge; Sicen Tao; Sicen Tao; Sicen Tao; Huanyang Chen; Huanyang Chen

doi:10.1364/OE.483094

1. Introduction

Hyperbolic materials, whose electric or magnetic tensor has a component along one principal axis with an opposite sign to the other principal component, display a unique hyperbolic dispersion curve. They have been extensively investigated from artificial hyperbolic metamaterials [1–6] to natural van der Waals (vdW) crystals [7–12], especially the new type of polar semiconductor, $\mathrm{\alpha}$-phase molybdenum trioxide ($\mathrm{\alpha} - \textrm{Mo}{\textrm{O}_\textrm{3}}$) [9–12], which supports highly confined in-plane hyperbolic phonon polaritons (PhPs) in the mid-infrared range. Due to their unique dispersion properties, hyperbolic materials usually exhibit high-momentum propagation, enhanced density of states, strong light confinement and long lifetimes [1,4–12], and have a great many applications, including, but not limited to super-resolution imaging [2], negative refraction [3], enhanced spontaneous emission [5], polariton manipulation [8–12]. Other brilliant effects such as the topological transition of PhPs from elliptic to hyperbolic dispersion can be achieved by twisted bilayer of $\mathrm{\alpha} - \textrm{Mo}{\textrm{O}_\textrm{3}}$ flakes [13,14] or dynamically tailored graphene/$\mathrm{\alpha} - \textrm{Mo}{\textrm{O}_\textrm{3}}$ heterostructures [15]. Similar effect also occurs in metric signature transition of space-times [16–19], which is akin to the exchange of space and time, giving the ideas of mimicking space-times with hyperbolic materials, and vice versa, providing more possibilities for hyperbolic wave control [18,19]. In a word, the increasing studies have shown the significance and potential of hyperbolic materials, however, there is an elusive issue remained in the numerical simulations of hyperbolic materials at frequency domain in the finite element software such as COMSOL Multiphysics (hereafter referred to as COMSOL), that is, a large amount of reflection and thus turbid images will be generated because the built-in perfectly matched layer (PML) is no longer applicable.

As a boundary condition, PML is an important component in the full-wave simulations to truncate the infinite domains into finite domains without influencing the propagation of electromagnetic waves [20–22]. Any electromagnetic wave transmitted from the computational domain can be attenuated in the PML without reflection, as schematically shown in Fig. 1(b), which is a zoom-in region of the simulation model in Fig. 1(a) (indicated by the dotted box). The model we consider in this paper is a rectangular computational domain with PML boundary conditions. However, as demonstrated in Fig. 1(c), when the background material in the computational domain is a hyperbolic material with permittivity and permeability parameters of $\{{{\mu_x} = 1,{\mu_y} ={-} 1,{\varepsilon_z} = 1} \}$ for two-dimensional (2D) transverse electric (TE) polarization, the result by COMSOL will show a messy field pattern. We next analyze the reason for this phenomenon.

Fig. 1. (a) The schematic model, where the number represents the domain number of the PML and the region bounded by the dotted line is (b). (b) The wave near the interface transmits from the computational domain into the domain 4 of the PML and decays exponentially. (c) The field pattern ($Re({{E_z}} )$) of a point source located at (0,0) in hyperbolic medium $\{{{\mu_x} = 1,{\mu_y} ={-} 1,{\varepsilon_z} = 1} \}$ by using the COMSOL’s built-in PML. (d) When the material of the computational domain has some losses, such as $\{{{\mu_x} = 1 + 0.001i,{\mu_y} ={-} 1 + 0.2i,{\varepsilon_z} = 1} \}$, the field pattern ($Re({{E_z}} )$) of the point source at (0,0) shows an incomplete image at a working frequency of 0.3 GHz.

Download Full Size | PDF

In the following, and without loss of generality, we focus on the 2D TE polarization in a Cartesian coordinate. The relative permittivity and permeability tensors $({\overleftrightarrow {{\varepsilon_r}},\overleftrightarrow {{\mu_r}}} )$ of biaxial materials can be expressed as:

(1)$$\overleftrightarrow {{\varepsilon _r}} = \left[ {\begin{array}{{ccc}} {{\varepsilon _x}}&0&0\\ 0&{{\varepsilon _y}}&0\\ 0&0&{{\varepsilon _z}} \end{array}} \right]\; ,\overleftrightarrow {{\mu _r}} = \left[ {\begin{array}{{ccc}} {{\mu _x}}&0&0\\ 0&{{\mu _y}}&0\\ 0&0&{{\mu _z}} \end{array}} \right]$$

Considering harmonic time dependence ${e^{i\omega t}}$ and supposing $\sqrt { - 1} = i$, we start from solving Maxwell’s equations {$\nabla \times \vec{H} = i\omega \overleftrightarrow {{\varepsilon _r}}{\varepsilon _0}\vec{E};\; \nabla \times \vec{E} ={-} i\omega \overleftrightarrow {{\mu _r}}{\mu _0}\vec{H}$}. For 2D TE polarization (${E_x} = {E_y} = {H_z} = 0$), the curl equations will be simplified as:

(2)$$\left\{ {\begin{array}{{c}} {\frac{{\partial {H_y}}}{{\partial x}} - \frac{{\partial {H_x}}}{{\partial y}} = i\omega {\varepsilon_z}{\varepsilon_0}{E_z}}\\ {\frac{{\partial {E_z}}}{{\partial y}} ={-} i\omega {\mu_x}{\mu_0}{H_x}}\\ {\frac{{\partial {E_z}}}{{\partial x}} = i\omega {\mu_y}{\mu_0}{H_y}} \end{array}\; } \right.\; \; $$

where the contributing material components are {${\mu _x}$, ${\mu _y}$, ${\varepsilon _z}$}. The out-of-plane electric field ${E_z}$ can be expressed as ${E_0}{e^{ - i({{k_x}x + {k_y}y} )}}$, we substitute it into the above formulas, and the dispersion relation of 2D TE waves in this material will be obtained as:

(3)$$\frac{{k_x^2}}{{{\mu _y}{\varepsilon _z}}} + \frac{{k_y^2}}{{{\mu _x}{\varepsilon _z}}} = k_0^2$$

where $k_0^2 = {\omega ^2}{\varepsilon _0}{\mu _0}$ being the vacuum wavevector. For transverse magnetic (TM) polarization (${H_x} = {H_y} = {E_z} = 0$), the contributing material components are {${\varepsilon _x}$, ${\varepsilon _y}$, ${\mu _z}$}, and Eq. (3) will become $\frac{{k_x^2}}{{{\varepsilon _y}{\mu _z}}} + \frac{{k_y^2}}{{{\varepsilon _x}{\mu _z}}} = k_0^2$, which only exchanges the dielectric and permeability parameters. Therefore, we just take TE mode as an example to study the feasibility of the designed PML. If ${\mu _x}{\varepsilon _z}\ast {\mu _y}{\varepsilon _z} < 0$, the isofrequency contour in this medium would be a hyperbola with one of the wave vector components being imaginary. Suppose that the field propagate in the computational domain is a plane wave ${e^{i\omega t - i{k_x}x - i{k_y}y}}$. The algorithm (complex coordinate stretching) of traditional PML makes the wave attenuate by adding an attenuation term $\alpha $ to the coordinate of propagation direction of the plane wave. For example, in the x direction, $\textrm{x}$ is replaced by x$+ \mathrm{i\alpha }(x )$, then the plane wave becomes ${e^{{k_x}\alpha (x )}}\ast {e^{i\omega t - i{k_x}x - i{k_y}y}}$, so that the wave will decay exponentially in the PML when ${k_x}\alpha (x )< 0$. However, this mapping will not work if the wave number k is an imaginary number. In terms of hyperbolic media, the waves will be expressed as ${e^{i\ast Im({{k_x}} )\alpha (x )}}\ast {e^{i\omega t + Im({{k_x}} )x - i{k_y}y}}$, ${k_x} = i\ast Im({{k_x}} )$ . Because one component of the wave vectors is imaginary, such as the ${k_x}$ component above, the original propagating term will become an exponential increasing term, and the added attenuation term $\alpha $ makes the wave attenuate no more but propagate until reaching the outer boundary, such as perfect electric conductor (PEC) here, and bounce back to the computational domain, thus resulting in reflections and scatterings.

In the past, different ways have been proposed to redesign PML for hyperbolic media [23]. For example, the general tools for first order hyperbolic systems were proposed [24]. A new construction of hyperbolic PML was designed to broaden the applicability [25]. Although these methods have precise numerical arguments, they are too specialized and tedious for direct application in simulations. On the other hand, the way of factitiously adding loss term in the material of the computational domain is well-known to all, thereinto, some of the cases are due to the loss of the material itself, like $\alpha - \textrm{Mo}{\textrm{O}_\textrm{3}}$ [10,11]. Although a clearer field pattern can be obtained in this way, like that in Fig. 1(d), the obvious blurred asymptote and incomplete and gradually faded image are inevitable, and besides, to make the waves attenuate to zero before reaching the boundary for avoiding reflection, larger losses will be required.

Bearing the advantages of flexible and high-efficiency in light field manipulation, transformation optics (TO) as a marvelous theoretical tool has been widely used in device designs [26–28]. The coordinate stretching method [22,29] can be well applied in TO to obtain the desired medium due to the form invariance of Maxwell’s equations, and is also used in non-reflection designs of PML [22,29,30]. But the mappings they suggested are only valid for traditional cases, whereas we will show in this paper that the mapping we used can perfectly accommodate hyperbolic materials.

In this paper, we utilize an imaginary coordinate mapping together with a complex coordinate stretching to design the material parameters of PML based on TO theory. The specific factors in the mapping and stretching are determined by analyzing the parameters of the computed hyperbolic materials. Through the analysis of the transformed plane wave expression, the attenuation term is designed to reach a relatively optimal effect. Our simulations of different cases show clear and complete field patterns, consistent with theoretical results and almost without reflections. In addition, the effects of the PML show a robustness against any thickness variation or source position changes, which provides more practicability for hyperbolic wave manipulation.

2. Design of PML

We first address our design according to TO theory by defining a coordinate mapping, which is derived in reverse from the calculated material. Assume the orthogonal coordinates in virtual space and physical space are $({X(x ),Y(y )} )$, and $({x,y} )$, respectively, and the material in virtual space is air, the TE mode material parameters of physical space are set as $\{{\mu_x^{\prime},\mu_y^{\prime},\varepsilon_z^{\prime}} \}$. The expressions of $\{{\mu_x^{\prime},\mu_y^{\prime},\varepsilon_z^{\prime}} \}$ can be calculated based on the TO theory as [26,28]

(4)$$\left\{ {\begin{array}{{c}} {\mu_x^{\prime} = \frac{{{\alpha_x}}}{{{\alpha_y}}}}\\ {\mu_y^{\prime} = \frac{{{\alpha_y}}}{{{\alpha_x}}}}\\ {\varepsilon_z^\mathrm{{\prime}} = \frac{1}{{{\alpha_x}{\alpha_y}}}} \end{array}} \right.$$

where ${\alpha _x} = dx/dX(x )\; $ and ${\alpha _y} = dy/dY(y )\; $ represent the Jacobi relations between virtual space and physical space. The dispersion relation of virtual space is $k_x^2 + k_y^2 = k_0^2$. From Eq. (4), the anisotropic refractive index components are:

(5)$$\left\{ {\begin{array}{{c}} {n_x^2 = \mu_y^{\prime}\varepsilon_z^\mathrm{{\prime}} = \frac{{{\alpha_y}}}{{{\alpha_x}}}\ast \frac{1}{{{\alpha_x}{\alpha_y}}} = 1/\alpha_x^2 = m}\\ {n_y^2 = \mu_x^{\prime}\varepsilon_z^{\prime} = \frac{{{\alpha_x}}}{{{\alpha_y}}}\ast \frac{1}{{{\alpha_x}{\alpha_y}}} = 1/\alpha_y^2 = n} \end{array}} \right.$$

where we set $\mu _y^{\prime}\varepsilon _z^\mathrm{{\prime}} = m,\; \mu _x^{\prime}\varepsilon _z^{\prime} = n$. Then, the dispersion relation of Eq. (3) can be rewritten as:

(6)$$\frac{{k_x^2}}{m} + \frac{{k_y^2}}{n} = k_0^2$$

Based on the complex coordinate stretching theory [22,29], the way of designing PML is to add an attenuation term to the coordinate. In this paper, we use such a coordinate mapping and stretching equations for PML:

(7)$$X(x )= \sqrt m x + A(x );Y(y )= \sqrt n y + B(y )$$

where the mappings $X(x )= \sqrt m x$, $Y(y )= \sqrt n y$ obey the anisotropic Fermat’s principle [19], resulting in the dispersion relation which is consistent with Eq. (6), and the stretching expressions $A(x )$ and $B(y )$ make the waves attenuate in x direction and y direction in the PML, respectively. While in the computational domain, $A(x )$ and $B(y )$ are equal to zero. Suppose that the electromagnetic fields propagate in the form of plane waves, which can be expressed as ${e^{i\omega t - i{k_x}X(x )- i{k_y}Y(y )}}$. As mentioned before, the wave vector component ${k_x}\; \; or\; \; {k_y}$ here can be imaginary when the physical space medium is hyperbolic medium with $m\ast n < 0$, which means that the mappings we proposed should contain an imaginary mapping along one coordinate, and the attenuation terms $A(x )$ or $B(y )$ cannot simply be treated as the pure imaginary expressions like $A(x )= i\sigma (x )$, $\sigma (x )\in R$. For the sake of later discussion, we use the notations of wave vectors $({k_x^{\prime},k_y^{\prime}} )$ instead of $({{k_x},{k_y}} )$ in physical space, where $({k_x^{\prime},k_y^{\prime}} )\; \in R$. Substitute the coordinate mapping and stretching equations (7) into the plane wave expression ${e^{i\omega t - i{k_x}X(x )- i{k_y}Y(y )}}$, we have:

(8)$${e^{ - ik_x^{\prime}\frac{{A(x )}}{{\sqrt m }}}}{e^{ - ik_y^{\prime}\frac{{B(y )}}{{\sqrt n }}}}{e^{i\omega t - ik_x^{\prime}x - ik_y^{\prime}y}}$$

where $k_x^{\prime} = {k_x}\sqrt m $, $k_y^{\prime} = {k_y}\sqrt n $. As shown from Eq. (8), the exponential terms $exp({ - ik_x^{\prime}x} )$ and $exp({ - ik_y^{\prime}y} )$ are the propagation terms. In order to make the wave decay in PML, the exponential terms $exp\left( { - ik_x^{\prime}\frac{{A(x )}}{{\sqrt m }}} \right)$ and $exp\left( { - ik_y^{\prime}\frac{{B(y )}}{{\sqrt n }}} \right)$ should be the attenuation terms. Next, we set the attenuation terms as $A(x )= \frac{{i\sqrt {|m |} }}{{\sqrt m }}{\sigma _x}(x )$ and $B(y )= \frac{{i\sqrt {|n |} }}{{\sqrt n }}{\sigma _y}(y )$ to represent the general expressions for different hyperbolic materials, where ${\sigma _x}(x )$ and ${\sigma _y}(y )$ are the attenuation functions that control the intensity of the attenuation terms, and ${\sigma _x}(x )\; ,\; {\sigma _y}(y )\; \in R$. We continue to discuss the role of these attenuation terms in PML design.

The model studied in this paper is a rectangle domain, as shown in Fig. 1(a), where a and b are the distances from the center of the computational domain to the inner boundary of PML, d is the thickness of PML layer. We set the center of the model as the origin of the coordinate. On account of the particularity of the rectangular model, PML is divided into eight domains as shown in Fig. 1(a). We need to keep the material parameters and plane wave energy continuous at the interfaces $|x |= a$ and $|y |= b$ for less reflections purpose. We discuss the selection of attenuation functions by taking ${\sigma _x}(x )$ as an example. Firstly, the attenuation function ${\sigma _x}(x )$ needs to satisfy the continuous condition:

(9)$${\sigma _x}(x ),\partial {\sigma _x}(x )/\partial x = \left\{ {\begin{array}{{c}} {{\sigma_x}(x ),\partial {\sigma_x}(x )/\partial x\; \; \; |x |> a}\\ {0\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; |x |= a} \end{array}} \right.$$

Besides, ${\sigma _x}(x )$ also needs to meet the requirement that $exp\left( { - ik_x^{\prime}\frac{{A(x )}}{{\sqrt m }}} \right)$ is a continuous monotone decreasing function in the interval of $|x |\in [{a,a + d} ]$. Here, we set it as $|{{\sigma_x}(x )} |= \rho \ast {\left( {\frac{{|x |- a}}{d}} \right)^2}$, where the constant $\rho $ is used to control the intensity of the attenuation function, the thickness d as the denominator keeps the attenuation of electromagnetic waves robust in different PML thickness. One question is, why do we choose quadratic function? On the one hand, quadratic function is enough to make the waves vanish in the PML before reaching the outer PEC boundary. On the other hand, quadratic function does not change drastically with increasing propagation distance of wave in the range of PML compared with the cubic or higher-degree function, which is conducive to reducing reflection or scattering caused by material discontinuities. In addition, the sign choice of the functions ${\sigma _x}(x )$ and ${\sigma _y}(y )$ should meet the conditions $- ik_x^{\prime}\frac{{A(x )}}{{\sqrt m }} < 0\; \; and\; - ik_y^{\prime}\frac{{B(y )}}{{\sqrt n }} < 0$, which can also be written as: $k_x^{\prime}\frac{{\sqrt {|m |} }}{m}{\sigma _x}(x )\le 0$ and $k_y^{\prime}\frac{{\sqrt {|n |} }}{n}{\sigma _y}(y )\le 0$. Combining all the above factors, we reach the final coordinate mapping and stretching expressions:

(10)$$\left\{ {\begin{array}{{c}} {X(x )= \sqrt m x \mp \frac{{\rho i\sqrt {|m |} }}{{\sqrt m }}{{\left( {\frac{{|x |- a}}{d}} \right)}^2}}\\ {Y(y )= \sqrt n y \mp \frac{{\rho i\sqrt {|n |} }}{{\sqrt n }}{{\left( {\frac{{|y |- b}}{d}} \right)}^2}} \end{array}} \right.$$

The sign “$\mp$” chooses “$+ $” (“$- $”) when the wave propagates along the positive (negative) direction of the coordinate. After adding up all the above considerations, we have the following Jacobi relations of the whole domain:

(11)$$\left\{ {\begin{array}{{c}} {{\alpha_x} = \frac{{dx}}{{dX}} = \frac{1}{{\sqrt m - \frac{{2\rho i\sqrt {|m |} }}{{\sqrt m }}\left( {\frac{{|x |- a}}{{{d^2}}}} \right)sgn({|x |- a} )}}}\\ {{\alpha_y} = \frac{{dy}}{{dY}} = \frac{1}{{\sqrt n - \frac{{2\rho i\sqrt {|n |} }}{{\sqrt n }}\left( {\frac{{|y |- b}}{{{d^2}}}} \right)sgn({|y |- b} )}}} \end{array}} \right.$$

where $sgn(x )$ is the step function, for $x > 0$, $sgn(x )= 1$, and for others $sgn(x )= 0$.

This PML can be easily applied to the electromagnetic simulation software, such as COMSOL, where the whole material parameters are set as Eq. (4), and $({{\alpha_x},{\alpha_y}} )$ are presented by Eq. (11) and $({m,n} )$ can be expressed by Eq. (5). It is worth mentioning that, it does not require to set both x and y attenuation in every domain [23]. In domains of 2,4,6,8, there are only attenuation perpendicular to the corresponding direction of the domains. As for the diagonal domains 1,3,5,7, because there are electromagnetic waves from adjacent PML domains, the PML needs attenuation along both x-axis and y-axis the same as those in adjacent domains. All expressions are available after the input studied material parameters $\{{\mu_x^{\prime},\mu_y^{\prime},\varepsilon_z^{\prime}} \}$ and the geometrical parameters of the rectangular model $({a,b,d} )$. For TM polarization, since the calculation process is exactly the same as that for TE polarization, we only need to replace $\{{\mu_x^{\prime},\mu_y^{\prime},\varepsilon_z^{\prime}} \}$ in Eq. (4) and (5) with $\{{\varepsilon_x^{\prime},\varepsilon_y^{\prime},\mu_z^{\prime}} \}$, which means $\{{m = \varepsilon_y^{\prime}\mu_z^{\prime};n = \varepsilon_x^{\prime}\mu_z^{\prime}} \}$. And the field patterns are changed to illustrate the magnetic field $Re({H_{z}} )$.

3. Simulation results

We now put our designed PML into COMSOL to verify its accuracy. The field patterns of a point source under different variables are illustrated, that is, material parameters, point source positions, and PML thicknesses, respectively. The working frequency is 0.3 GHz. The constant ρ can be slightly adjusted according to the specific material. Because when ρ=0.5, the simulation versatility is higher, we choose ρ=0.5 for all the simulations in our paper. The model of the computational domain is a square with side length of 10 meters, i.e., $a = b = 5m$ (could be arbitrary unit).

Firstly, the thickness of PML is fixed to $d = 1$m, the electric field patterns ($Re({{E_z}} )$) of a point source at (0, 0) in the materials of $\{{{\mu_x} = 1,{\mu_y} ={-} 1,{\varepsilon_z} = 1} \}$, $\{{{\mu_x} = 2,{\mu_y} ={-} 1,{\varepsilon_z} = 1} \}$ and $\{{{\mu_x} = 1,{\mu_y} ={-} 2,{\varepsilon_z} = 1} \}$ by COMSOL are shown as Fig. 2(a), (b), and (c), respectively. In contrast with the field pattern by the way of adding loss in the computational domain (Fig. 1(d)), the field pattern with our PML in Fig. 2(a) is obviously more complete. Although there are slight ripples near the point source, they can be weakened by increasing the mesh or improving the accuracy of the software to keep them within tolerance. To better illustrate the effect of the PML, we analytically present the wave patterns of analytical results for comparison. Helmholtz equation is solved to obtain the corresponding solution in cylindrical coordinates, which is the first kind of zero-order Hankel function $\textrm{}H_0^{(1 )}\; ({{k_0}\mathrm{\ast }r} )$. The radial coordinate is expressed as $r = \sqrt {{{({Re[{X(x )} ]} )}^2} + {{(Re[Y(y )])}^2}} $. After the coordinate mappings $X(x )= \sqrt m x,Y(y )= \sqrt n y$, the solution is changed to be $H_0^{(1 )}\; \left( {{k_0}\mathrm{\ast }\sqrt {m{x^2} + m{y^2}} } \right)$. The corresponding results with the same material parameters are showing in Fig. 2 (d), (e) and (f), respectively, reflecting the coherence of the field patterns under the effect of our designed PML. Besides, we can get a full and clear wave pattern and asymptotic boundary without enlarging the size of the computational domain or enhancing the strength of the attenuation term.

Fig. 2. Electric field patterns ($Re({{E_z}} )$) from the point source at (0,0) with materials (a) [(d)] $\{{{\mu_x} = 1,{\mu_y} ={-} 1,{\varepsilon_z} = 1} \}$, (b) [(e)] $\{{{\mu_x} = 2,{\mu_y} ={-} 1,{\varepsilon_z} = 1} \}$ and (c) [(f)] $\{{{\mu_x} = 1,{\mu_y} ={-} 2,{\varepsilon_z} = 1} \}$ by COMSOL (or analytically). The thickness of the PML is 1 m.

Download Full Size | PDF

We further investage the influence of the point source position and the PML’s thickness on the effect to hyperbolic materials. We keep the material unchanged as $\{{{\mu_x} = 1,{\mu_y} ={-} 1,{\varepsilon_z} = 1} \}$, when the thickness of the PML is $d = 1$m, the field patterns of point source at (-5,0) and (5,-5) are shown in Fig. 3(a) and (b). Then we fix the point souce at (5,-5), and demonstrate the field patterns with the thickness of PML changing from 2 m to 0.5 m, as shown in Fig. 3(c) and (d). As a result, all the above images display good waveforms with no reflection. It indicates that our PML can be well applied in most of the cases with robustness under the variation of the point source location and the thickness of the PML.

Fig. 3. Electric field patterns ($Re({{E_z}} )$) of material $\{{{\mu_x} = 1,{\mu_y} ={-} 1,{\varepsilon_z} = 1} \}$, with the thickness of PML being 1 m from the point source at (a)[(b)] (0,-5) [ (5,-5)] by COMSOL, and with the thickness of PML (c) [(d)] being 2 m (0.5 m) from the point source at (5,-5) by COMSOL.

Download Full Size | PDF

4. Conclusion

In this paper, a new TO-based PML is proposed for biaxial hyperbolic materials through an imaginary coordinate mapping and complex coordinate stretching. The attenuation term is clearly analyzed from the plane wave expressions. Compared with the method of adding loss in the computed materials components, our PML provides clear and complete wave patterns. Besides, the PML can be compatible with different material parameters, various thickness of PML and changed point source locations. The proposed PML can be regarded as a platform for hyperbolic material simulations, the design methods can be further extended for PML in cylindrical or spherical coordinates, and other fields for wave simulations, such as acoustics, water waves, etc.

Funding

National Natural Science Foundation of China (92050102); National Key Research and Development Program of China (2020YFA0710100); Natural Science Foundation of Jiangxi Province (20224ACB201005).Fundamental Research Funds for the Central Universities (20720200074, 20720220033, 20720220134).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. Z. Jacob, L. V. Alekseyev, and E. Narimanov, “Optical hyperlens: far-field imaging beyond the diffraction limit,” Opt. Express 14(18), 8247 (2006). [CrossRef]

2. Z. Liu, H. Lee, Y. Xiong, C. Sun, and X. Zhang, “Far-field optical hyperlens magnifying sub-diffraction-limited objects,” Science 315(5819), 1686 (2007). [CrossRef]

3. J. Yao, Z. Liu, Y. Liu, Y. Wang, C. Sun, G. Bartal, A. M. Stacy, and X. Zhang, “Optical negative refraction in bulk metamaterials of nanowires,” Science 321(5891), 930 (2008). [CrossRef]

4. Z. Jacob, J.-Y. Kim, G. V. Naik, A. Boltasseva, E. E. Narimanov, and V. M. Shalaev, “Engineering photonic density of states using metamaterials,” Appl. Phys. B 100(1), 215–218 (2010). [CrossRef]

5. M. A. Noginov, H. Li, Y. A. Barnakov, D. Dryden, G. Nataraj, G. Zhu, C. E. Bonner, M. Mayy, Z. Jacob, and E. E. Narimanov, “Controlling spontaneous emission with metamaterials,” Opt. Lett. 35(11), 1863–1865 (2010). [CrossRef]

6. A. Poddubny, I. Iorsh, P. Belov, and Y. Kivshar, “Hyperbolic metamaterials,” Nat. Photonics 7(12), 948–957 (2013). [CrossRef]

7. D. Lee, S. So, G. Hu, M. Kim, T. Badloe, H. Cho, J. Kim, H. Kim, C.-W. Qiu, and J. Rho, “Hyperbolic metamaterials: fusing artificial structures to natural 2D materials,” eLight 2(1), 1 (2022). [CrossRef]

8. S. Dai, Z. Fei, Q. Ma, S. Rodin, M. Wagner, A. S. Mcleod, M. K. Liu, W. Gannett, W. Regan, K. Watanabe, T. Taniguchi, M. Thiemens, G. Dominguez, A. H. Castro Neto, A. Zettl, F. Keilmann, P. Jarillo-Herrero, M. M. Fogler, and D. N. Basov, “Tunable phonon polaritons in atomically thin van der Waals crystals of boron nitride,” Science 343(6175), 1125–1129 (2014). [CrossRef]

9. P. Li, I. Dolado, F. J. Alfaro-Mozaz, F. Casanova, L. E. Hueso, S. Liu, J. H. Edgar, A. Y. Nikitin, S. Vélez, and R. Hillenbrand, “Infrared hyperbolic metasurface based on nanostructured van der Waals materials,” Science 359(6378), 892–896 (2018). [CrossRef]

10. W. Ma, P. Alonso-González, S. Li, A. Y. Nikitin, J. Yuan, J. Martín-Sánchez, J. Taboada-Gutiérrez, I. Amenabar, P. Li, S. Vélez, C. Tollan, Z. Dai, Y. Zhang, S. Sriram, K. Kalantar-Zadeh, S.-T. Lee, R. Hillenbrand, and Q. Bao, “In-plane anisotropic and ultra-low-loss polaritons in a natural van der Waals crystal,” Nature 562(7728), 557–562 (2018). [CrossRef]

11. Z. Zheng, N. Xu, S. L. Oscurato, M. Tamagnone, F. Sun, Y. Jiang, Y. Ke, J. Chen, W. Huang, W. L. Wilson, A. Ambrosio, S. Deng, and H. Chen, “A mid-infrared biaxial hyperbolic van der Waals crystal,” Sci. Adv. 5(5), eaav8690 (2019). [CrossRef]

12. F. Sun, W. Huang, Z. Zheng, N. Xu, Y. Ke, R. Zhan, H. Chen, and S. Deng, “Polariton waveguide modes in two-dimensional van der Waals crystals: an analytical model and correlative nano-imaging,” Nanoscale 13(9), 4845–4854 (2021). [CrossRef]

13. G. Hu, Q. Ou, G. Si, Y. Wu, J. Wu, Z. Dai, A. Krasnok, Y. Mazor, Q. Zhang, Q. Bao, C.-W. Qiu, and A. Alù, “Topological polaritons and photonic magic angles in twisted α-"Mo” “O” _"3” bilayers,” Nature 582(7811), 209–213 (2020). [CrossRef]

14. J. Duan, N. Capote-Robayna, J. Taboada-Gutiérrez, G. Álvarez-Pérez, I. Prieto, J. Martín-Sánchez, A. Y. Nikitin, and P. Alonso-González, “Twisted nano-optics: Manipulating light at the nanoscale with twisted phonon polaritonic slabs,” Nano Lett. 20(7), 5323–5329 (2020). [CrossRef]

15. Y. Zeng, Q. Ou, L. Liu, C. Zheng, Z. Wang, Y. Gong, X. Liang, Y. Zhang, G. Hu, Z. Yang, C. W. Qiu, Q. Bao, H. Chen, and Z. Dai, “Tailoring Topological Transitions of Anisotropic Polaritons by Interface Engineering in Biaxial Crystals,” Nano Lett. 22(10), 4260–4268 (2022). [CrossRef]

16. I. I. Smolyaninov and E. E. Narimanov, “Metric signature transitions in optical metamaterials,” Phys. Rev. Lett. 105(6), 067402 (2010). [CrossRef]

17. S. Fumeron, B. Berche, F. Santos, E. Pereira, and F. Moraes, “Optics near a hyperbolic defect,” Phys. Rev. A 92(6), 063806 (2015). [CrossRef]

18. S. Dehdashti, A. Shahsafi, B. Zheng, L. Shen, Z. Wang, R. Zhu, H. Chen, and H. Chen, “Conformal hyperbolic optics,” Phys. Rev. Res. 3(3), 033281 (2021). [CrossRef]

19. S. Tao, T. Hou, Y. Zeng, G. Hu, Z. Ge, J. Liao, S. Zhu, T. Zhang, C. W. Qiu, and H. Chen, “Anisotropic Fermat’s principle for controlling hyperbolic van der Waals polaritons,” Photonics Res. 10(10), B14–B22 (2022). [CrossRef]

20. J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of computational physics 114.2 (1994): 185-200.

21. Z. S. Sacks, D. M. Kingsland, R. Lee, and J. F. Lee, “A perfectly matched anisotropic absorber for use as an absorbing boundary condition,” IEEE Trans. Antennas Propag. 43(12), 1460–1463 (1995). [CrossRef]

22. W. C. Chew and W. H. Weedon, “A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates,” Microw. Opt. Technol. Lett. 7(13), 599–604 (1994). [CrossRef]

23. K. Duru and G. Kreiss, “The perfectly matched layer (PML) for hyperbolic wave propagation problems: A review,” arXiv, arXiv:2201.03733(2022).

24. D. Appelö, T. Hagstrom, and G. Kreiss, “Perfectly matched layers for hyperbolic systems: general formulation, well-posedness, and stability,” SIAM J. Appl. Math. 67(1), 1–23 (2006). [CrossRef]

25. T. Hagstrom, “A new construction of perfectly matched layers for hyperbolic systems with applications to the linearized Euler equations,” Mathematical and Numerical Aspects of Wave Propagation WAVES 2003. Springer, Berlin, Heidelberg, 2003. 125-129.

26. J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science 312(5781), 1780–1782 (2006). [CrossRef]

27. U. Leonhardt, “Optical conformal mapping,” Science 312(5781), 1777–1780 (2006). [CrossRef]

28. H. Chen, C. T. Chan, and P. Sheng, “Transformation optics and metamaterials,” Nat. Mater. 9(5), 387–396 (2010). [CrossRef]

29. L. Knockaert and D. De Zutter, “On the stretching of Maxwell's equations in general orthogonal coordinate systems and the perfectly matched layer,” Microw. Opt. Technol. Lett. 24(1), 31–34 (2000). [CrossRef]

30. L. Zhao and A. C. Cangellaris, “A general approach for the development of unsplit-field time-domain implementations of perfectly matched layers for FDTD grid truncation,” IEEE Microwave and Guided Wave Letters 6(5), 209–211 (1996). [CrossRef]

Perfectly matched layer for biaxial hyperbolic materials

Abstract

1. Introduction

2. Design of PML

3. Simulation results

4. Conclusion

Funding

Disclosures

Data availability

References

Data availability

Cited By

Figures (3)

Equations (11)

Optics Express