Finite element modeling of electromagnetic properties in photonic bianisotropic structures

Abstract

Given a constitutive relation of the bianisotropic medium, it is not trivial to study how light interacts with the photonic bianisotropic structure due to the limited available means of studying electromagnetic properties in bianisotropic media. In this paper, we study the electromagnetic properties of photonic bianisotropic structures using the finite element method. We prove that the vector wave equation with the presence of bianisotropic is self-adjoint under scalar inner product. we propose a balanced formulation of weak form in the practical implementation, which outperforms the standard formulation in finite element modeling. Furthermore, we benchmark our numerical results obtained from finite element simulation in three different scenarios. These are bianisotropy-dependent reflection and transmission of plane waves incident onto a bianisotropic slab, band structure of bianisotropic photonic crystals with valley-dependent phenomena, and the modal properties of bianisotropic ring resonators. The first two simulated results obtained from our modified weak form yield excellent agreements either with theoretical predictions or available data from the literature, and the modal properties in the last example, i.e., bianisotropic ring resonators as a polarization-dependent optical insulator, are also consistent with the theoretical analyses.

Graphical abstract

Keywords

bianisotropic / finite element method / adjoint

1 Introduction

Recent progress in chiral photonics, i.e., bianisotropic metamaterials [1–7] and bianisotropic metasurfaces [8–15], has significantly advanced our understanding of light transport in complex photonic structures and spurred numerous applications. Indeed, bianisotropic metamaterials have been used to realize novel resonators [1,2,6] and construct photonic topological insulators [4,5,7]. Bianisotropic metasurfaces have been used to manipulate the polarization of light [8,9] and to achieve exotic refraction and transmission of light [10–14]. In metamaterials or metasurfaces, it is common to use optically complicated structures inside each unit cell together with a very large number of inclusions to obtain certain functionalities. Thus the complete multiple scattering of light among all the inclusions is hardly traceable. Therefore, under proper approximation, the chiral metamaterials [15–20] or chiral metasurfaces [21,22] can be described using effective constitutive parameters (permittivity and permeability), which in principle accounts for the bianisotropic response of the complicated unit cell. Notably, the grading effect in a wide range of metasurfaces, e.g., with a rotated orientation of each unit cell structure, can also be included in the spatial dependent bianisotropic constitutive parameters in the effective medium description.

However, the available means of calculating photonic bianisotropic medium is limited within analytical methods, which are usually complicated and only apply to simple structures, i.e., interfaces or slab structures [23–26]. A generic numerical approach that calculates the electromagnetic properties for bianisotropic structure with arbitrary shapes is lacking. In this paper, we intend to fill this gap and propose a finite element method (FEM) to simulate the optical properties of bianisotropic media.

The paper is organized as follows: In Section 2, we derive the balanced formulation of the weak form for the finite element implementation. In Sections 3−5, we illustrate and validate our numerical approach via three examples, i.e., bianisotropic slab structure, bianisotropic photonic crystal, and bianisotropic ring resonator respectively. Finally, Section 6 concludes the paper.

2 Model

The constitute relations of the generic bianisotropic materials are given by \mathit{D}=\overline{\mathit{ϵ}}\mathit{E}+{\overline{\mathit{\chi }}}_{eh}\mathit{H}, \mathit{B}=\overline{\mathit{\mu }}\mathit{H}+{\overline{\mathit{\chi }}}_{he}\mathit{E}, where \mathit{D} and \mathit{B} are electric displacement vector and the magnetic induction intensity respectively, \mathit{E} and \mathit{H} are electric and magnetic field respectively, \overline{\mathit{\epsilon }}={\epsilon }_{\mathrm{0}}{\overline{\mathit{\epsilon }}}_r and \overline{\mathit{\mu }}={\mu }_{\mathrm{0}}{\overline{\mathit{\mu }}}_r are permittivity and permeability, respectively, {\overline{\mathit{\chi }}}_{eh}=\sqrt{{\epsilon }_{\mathrm{0}}{\mu }_{\mathrm{0}}}{\overline{\mathit{\chi }}}_{eh}^r and {\overline{\mathit{\chi }}}_{he}=\sqrt{{\epsilon }_{\mathrm{0}}{\mu }_{\mathrm{0}}}{\overline{\mathit{\chi }}}_{he}^r are coupling constants. As for the source-free Maxwell’s equations with time-harmonic dependence of {\mathrm{e}}^{\mathrm{i}\omega t}, i.e., \nabla \times \mathit{E}=-\mathrm{i}\omega \mathit{B}, \nabla \times \mathit{H}=\mathrm{i}\omega \mathit{D}, {\mathrm{i}}^{\mathrm{2}}=-\mathrm{1}, one reformulated Maxwell’s equations using the normalized field, i.e., \mathit{e}(\mathit{r})=\mathit{E}(\mathit{r}) and \mathit{h}(\mathit{r})=\sqrt{{\displaystyle \frac{{\mu }_{\mathrm{0}}}{{\epsilon }_{\mathrm{0}}}}}\mathit{H}(\mathit{r}), as follows

(1)\begin{array}{l}\begin{array}{rr}\hfill \left(\nabla \times +\mathrm{i}{k}_{\mathrm{0}}{\overline{\mathit{\chi }}}_{he}^r\right)\mathit{e}(\mathit{r})+\mathrm{i}{k}_{\mathrm{0}}{\overline{\mathit{\mu }}}_r\mathit{h}(\mathit{r})& \hfill =\mathrm{0},\\ \hfill \left(\nabla \times -\mathrm{i}{k}_{\mathrm{0}}{\overline{\mathit{\chi }}}_{eh}^r\right)\mathit{h}(\mathit{r})-\mathrm{i}{k}_{\mathrm{0}}{\overline{\mathit{\epsilon }}}_r\mathit{e}(\mathit{r})& \hfill =\mathrm{0.}\end{array}\end{array}

Eliminating the magnetic field \mathit{h}(\mathit{r}), one arrives at the vector wave equation given as follows

(2)\begin{array}{l}\left(\nabla \times -\mathrm{i}{k}_{\mathrm{0}}{\overline{\mathit{\chi }}}_{eh}^r\right)\left[{\displaystyle \frac{\mathrm{1}}{{\overline{\mathit{\mu }}}_r}}\left(\nabla \times +\mathrm{i}{k}_{\mathrm{0}}{\overline{\mathit{\chi }}}_{he}^r\right)\right]\mathit{e}(\mathit{r})-k_{\mathrm{0}}^{\mathrm{2}}{\overline{\mathit{\epsilon }}}_r\mathit{e}(\mathit{r})=\mathrm{0.}\end{array}

Once the electric \mathit{e} field known, the magnetic field \mathit{h} is given by \mathit{h}(\mathit{r})=-{\displaystyle \frac{\mathrm{1}}{\mathrm{i}{k}_{\mathrm{0}}}}{\overline{\mathit{\mu }}}_r^{-\mathrm{1}}\nabla \times \mathit{e}(\mathit{r})-{\overline{\mathit{\mu }}}_r^{-\mathrm{1}}{\overline{\mathit{\chi }}}_{he}^r\mathit{e}(\mathit{r}). As for Eq. (2), the vectorial wave equation is essentially determined by operator \mathit{L}, i.e., \mathit{L}=\left(\nabla \times -\mathrm{i}{k}_{\mathrm{0}}{\overline{\mathit{\chi }}}_{eh}^r\right)\left[{\displaystyle \frac{\mathrm{1}}{{\overline{\mathit{\mu }}}_r}}\left(\nabla \times +\mathrm{i}{k}_{\mathrm{0}}{\overline{\mathit{\chi }}}_{he}^r\right)\right]-k_{\mathrm{0}}^{\mathrm{2}}{\overline{\mathit{\epsilon }}}_r. Under scaler inner product, we show that the operator L is indeed self-adjoint, which is described as follows

(3)\begin{array}{l}(\mathit{F},L\mathit{E})=(L\mathit{F},\mathit{E}),\end{array}

see proof in Supplementary Material A. Importantly, we can further proof that

(4)\begin{array}{l}{\int }^{\mathrm{}}\mathrm{d}V\left(\left[{\displaystyle \frac{\mathrm{1}}{{\overline{\mathit{\mu }}}_r}}\left(\nabla \times +\mathrm{i}{k}_{\mathrm{0}}{\overline{\mathit{\chi }}}_{he}^r\right)\right]\mathit{F}\cdot {\left[\nabla \times -\mathrm{i}{k}_{\mathrm{0}}{\overline{\mathit{\chi }}}_{eh}^r\right]}^{\mathrm{T}}\mathit{E}-k_{\mathrm{0}}^{\mathrm{2}}{\overline{\mathit{\epsilon }}}_r\mathit{F}\cdot \mathit{E}\right)=(\mathit{F},L\mathit{E})\\ =(\mathit{L}F,\mathit{E})={\int }^{\mathrm{}}\mathrm{d}V\left(\left[{\displaystyle \frac{\mathrm{1}}{{\overline{\mathit{\mu }}}_r}}\left(\nabla \times +\mathrm{i}{k}_{\mathrm{0}}{\overline{\mathit{\chi }}}_{he}^r\right)\mathit{E}\right]\cdot \left[\nabla \times -\mathrm{i}{k}_{\mathrm{0}}{\left({\overline{\mathit{\chi }}}_{eh}^r\right)}^{\mathrm{T}}\right]\mathit{F}-k_{\mathrm{0}}^{\mathrm{2}}{\overline{\mathit{\epsilon }}}_r\mathit{E}\cdot \mathit{F}\right),\end{array}

which is coined as the balanced weak form. This balanced weak form contains a few implications that deserves further discussion: (1) In practical implementation of FEM, the Galakin procedures is usually adopted. In the Galakin procedures, the test function \mathit{F} is by default selected as the basis functions that are used to interpolate the electric field. Thus, it is essential to keep the test function space and the expansion function spaces undergone the same operations, i.e., transformation in exactly the same manner; (2) Applying standard FEM procedure leads to the unbalanced formulation of weak form in bianisotropic media as used in reference [27], which gives rise to artificial numerical errors due to the spatial differentiation of constitutive parameters, i.e., giving rise to the unusual imaginary part of the propagation constant of bianisotropic waveguides; (3) In this paper, the balanced formulation is adopted, it turns to be extremely important in the finite element modeling in bianisotropic medium and overcome the problem in unbalanced formulation from the comparison between the standard weak form and balanced weak form formulation that is not shown in the paper.

3 Reflection and transmission of light in layered bianisotropic medium

First, we benchmark our finite element model of light reflection and transmission against the analytical calculations. We consider a slab geometry in which light is incident in normal direction, as shown in Fig. 1(a). The slab in dark gray surrounded with air in light gray has {\epsilon }_r=\mathrm{4}, {\mu }_r=\mathrm{1} and the magnetoelectric coupling constants

{\overline{\mathit{\chi }}}_{he}^r=-{\left({\overline{\mathit{\chi }}}_{eh}\right)}^{\mathrm{T}}=\left(\begin{array}{lll}\mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{i}{\chi }_{\mathrm{33}}\end{array}\right).

The slab has a finite thickness L=\mathrm{1}\mathrm{\hspace{0.17em}}\mathrm{m} in horizontal direction, and the incident light contains in-plane (p) polarization and out-of-plane (s) polarization. The practical implementation is realized by modifying the original weak form in COMSOL Multiphysics [28] into our balanced weak form (the right term in Eq. (4)) for the bianisotropic medium.

The numerical calculation of the {\chi }_{\mathrm{33}} dependent reflection/transmission spectrum at vacuum wavelength \lambda =\mathrm{1}\mathrm{\hspace{0.17em}}\mathrm{m} from our reformulated weak form are perfectly coincides with semi-analysis method (see details in Supplementary Material B) displayed in Figs. 1(b)−1(e). As a side remark, there is no difference between p and s polarization for normal incidence at {\chi }_{\mathrm{33}}=\mathrm{0}. As for nonzero {\chi }_{\mathrm{33}}, the reflection/transmission of p and s polarization differ from each other. This is reasonable since, inside the bianisotropic slab, the p and s polarization are not eigen-polarization and gets mixed due to {\chi }_{\mathrm{33}}, and the true eigen-polarization essentially has different effective refractive index. As a result, the incident p and s polarization can be decomposed into the two eigen-polarizations that has different refractive index depending on {\chi }_{\mathrm{33}}, which leads to the variation of reflection and transmission against {\chi }_{\mathrm{33}}, as confirmed both by our numerical simulation and analytical calculations.

4 Band structure of bianisotropic photonic crystals

Bianisotropic metamaterials and metasurfaces have been used to realize photonic topological insulator or photonic valley Hall effect. Dong and his coworkers [29] explored the photonic valley in honeycomb photonic crystals, in which the inversion symmetry is broken via the material bianisotropy displayed in Fig. 2(a). The unit cell in this honeycomb lattice has two different rods with radii r in purple and blue, respectively, and lattice constant a. These two rods have uniaxial permittivity and permeability {\bar{\epsilon }}_r={\bar{\mu }}_r=\mathrm{diag}({\epsilon }_{//},{\epsilon }_{//},{\epsilon }_{\perp }), and nonzero bianisotropic tensor

{\overline{\mathit{\chi }}}_{eh}^r={\overline{\mathit{\chi }}}_{he}^r=\left(\begin{array}{lll}\mathrm{0}& \mathrm{i}\kappa {\epsilon }_{//}& \mathrm{0}\\ -\mathrm{i}\kappa {\epsilon }_{//}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right),

where the parameter \kappa in rods 1 and 2 have opposite sign. With the modified weak form, we calculated band structure of this photonic lattice in Fig. 2(b), where the dispersion curves for two pseudo-spins (spin-up and spin-down) of light on the high symmetry line of Brillouin zone are shown and have perfect agreement with the data from Ref. [29]. The two spins have almost the same band structure in low frequency, but the spin-up state has a smaller bandgap near K and a larger bandgap around K\mathrm{\text{'}} than that of the spin-down state. As a result, in the range of frequency from 0.375 to 0.435, the bandgap vanishes at K valley for the spin-up state; the same is true for spin-down state at K\mathrm{\text{'}} valley. Thus, the two spin states propagate in different directions determined by the associated valley index, as illustrated with two arrows in Fig. 2(a).

5 Bianisotropic ring resonator

We consider a bianisotropic slab waveguide as shown in Fig. 3(a), which extends infinitely along x and z direction and is sandwiched by air. The bianisotropic slab has a relative permittivity {\epsilon }_r=\mathrm{1}, relative permeability {\mu }_r=\mathrm{1}, and magnetoelectric coupling tensor

{\bar{\chi }}_{eh}^r=-{\left({\bar{\chi }}_{he}^r\right)}^{\mathrm{T}}=\left(\begin{array}{lll}\mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{i}{\chi }_{\mathrm{33}}\end{array}\right).

Light propagates along x, the effective refractive index of the fundamental mode increases for larger {\chi }_{\mathrm{33}}, as can be seen from the simplified wave equation in the presence of {\overline{\mathit{\chi }}}_{eh}^r in Supplementary Material C. The fundamental mode is linearly polarized for vanished {\chi }_{\mathrm{33}}, and is elliptically polarized for nonzero {\chi }_{\mathrm{33}}. By time reversal symmetry, the polarization of light will be converted into its time-reversal partner, i.e., the helicity flipping sign, for opposite propagating direction. As the slab waveguide is rolled to form a ring resonator shown in Fig. 3(b), the x, y and z components of an electric field in an ordinary slab waveguide corresponds to the azimuthal (E_{\theta }), radial (E_r) and z components. Therefore, light propagates along the azimuthal direction in the ring resonator is the same as it transports along x in the slab waveguide. Thus, light propagates clockwise or counter-clockwise, the polarization of light will have opposite helicity. Consequently, the ring resonator could function as a polarization depended insulator as shown in Figs. 3(c)−3(f). The clockwise and counter-clockwise modes in the ring resonator have opposite helicity of polarization in electric field basis (E_z\mp \mathrm{1.72}\mathrm{i}{E}_{\theta }) at vacuum wavelength \mathrm{1.0615}\mathrm{\hspace{0.17em}}\mathrm{\mu m}. A slab waveguide with the same permittivity and permeability supports transverse electric and transverse magnetic modes, which span the degenerate basis for the polarization. In Figs. 3(c) and 3(d), the polarized optical beam E_z+\mathrm{1.72}\mathrm{i}{E}_y is excited at the left and the right ports of the slab waveguide respectively, it couples with counter-clockwise mode E_z+\mathrm{1.72}\mathrm{i}{E}_{\theta } in the ring resonator with positive direction propagation (from left to right) while it propagates through the slab without any coupling with the ring from right to left. In contrast, its time-reversal partner, i.e., E_z-\mathrm{1.72}\mathrm{i}{E}_y, propagates without coupling in a positive direction but couples with clockwise mode E_z-\mathrm{1.72}\mathrm{i}{E}_{\theta } in negative propagation as shown in Figs. 4(e) and 4(f). Therefore, the polarized beams E_z\pm \mathrm{1.72}\mathrm{i}{E}_y are insulated in positive and negative directions, respectively. Furthermore, the polarization ellipticity of the modes in the ring can be manipulated by changing {\chi }_{\mathrm{33}}, this structure provides a platform to isolate or pick out any ellipse polarized beam in a simple way. Evidently, our numerical results are fully consistent with the earlier polarization-dependent optical insulation based on the bianisotropic ring resonator.

6 Conclusion

We propose a numerical method based on FEM approach to study the optical properties of complex bianisotropic structures. We first benchmark our method in reflection and transmission of a bianisotropic slab, and further validate our numerical calculation of the band structure in a bianisotropic photonic crystal that plays a relevant role in topological photonics. Our results have perfectly coincided either with the analytical results or the literature data. Finally, we illustrate that our method can be used to design a bianisotropic ring resonator, which functions as a polarization-dependent optical insulator.

References

[1] Gansel J K, Thiel M, Rill M S, Decker M, Bade K, Saile V, von Freymann G, Linden S, Wegener M. Gold helix photonic metamaterial as broadband circular polarizer. Science, 2009, 325(5947): 1513–1515

[2] Kriegler C E, Rill M S, Linden S, Wegener M M. Bianisotropic photonic metamaterials. IEEE Journal of Selected Topics in Quantum Electronics, 2010, 16(2): 367–375

[3] Soukoulis C M, Wegener M. Past achievements and future challenges in the development of three-dimensional photonic metamaterials. Nature Photonics, 2011, 5(9): 523–530

[4] Guo Q, Gao W, Chen J, Liu Y, Zhang S. Line degeneracy and strong spin-orbit coupling of light with bulk bianisotropic metamaterials. Physical Review Letters, 2015, 115(6): 067402

[5] Slobozhanyuk A P, Khanikaev A B, Filonov D S, Smirnova D A, Miroshnichenko A E, Kivshar Y S. Experimental demonstration of topological effects in bianisotropic metamaterials. Scientific Reports, 2016, 6(1): 22270

[6] Moritake Y, Tanaka T. Bi-anisotropic Fano resonance in three-dimensional metamaterials. Scientific Reports, 2018, 8(1): 9012

[7] Yu Y Z, Kuo C Y, Chern R L, Chan C T. Photonic topological semimetals in bianisotropic metamaterials. Scientific Reports, 2019, 9(1): 18312

[8] Pfeiffer C, Grbic A. Bianisotropic metasurfaces for optimal polarization control: analysis and synthesis. Physical Review Applied, 2014, 2(4): 044011

[9] Pfeiffer C, Zhang C, Ray V, Guo L J, Grbic A. Polarization rotation with ultra-thin bianisotropic metasurfaces. Optica, 2016, 3(4): 427–432

[10] Pfeiffer C, Zhang C, Ray V, Guo L J, Grbic A. High performance bianisotropic metasurfaces: asymmetric transmission of light. Physical Review Letters, 2014, 113(2): 023902

[11] Odit M, Kapitanova P, Belov P, Alaee R, Rockstuhl C, Kivshar Y S. Experimental realisation of all-dielectric bianisotropic metasurfaces. Applied Physics Letters, 2016, 108(22): 221903

[12] Epstein A, Eleftheriades G V. Arbitrary power-conserving field transformations with passive lossless omega-type bianisotropic metasurfaces. IEEE Transactions on Antennas and Propagation, 2016, 64(9): 3880–3895

[13] Asadchy V S, Díaz-Rubio A, Tretyakov S A. Bianisotropic metasurfaces: physics and applications. Nanophotonics, 2018, 7(6): 1069–1094

[14] Li J, Shen C, Díaz-Rubio A, Tretyakov S A, Cummer S A. Systematic design and experimental demonstration of bianisotropic metasurfaces for scattering-free manipulation of acoustic wavefronts. Nature Communications, 2018, 9(1): 1342

[15] Li J, Díaz-Rubio A, Shen C, Jia Z, Tretyakov S A, Cummer S. Highly efficient generation of angular momentum with cylindrical bianisotropic metasurfaces. Physical Review Applied, 2019, 11(2): 024016

[16] Chen X, Wu B I, Kong J A, Grzegorczyk T M. Retrieval of the effective constitutive parameters of bianisotropic metamaterials. Physical Review E, 2005, 71(4): 046610

[17] Li Z, Aydin K, Ozbay E. Determination of the effective constitutive parameters of bianisotropic metamaterials from reflection and transmission coefficients. Physical Review E, 2009, 79(2): 026610

[18] Ouchetto O, Qiu C W, Zouhdi S, Li L W, Razek A. Homogenization of 3-D periodic bianisotropic metamaterials. IEEE Transactions on Microwave Theory and Techniques, 2006, 54(11): 3893–3898

[19] Hasar U C, Muratoglu A, Bute M, Barroso J J, Ertugrul M. Effective constitutive parameters retrieval method for bianisotropic metamaterials using waveguide measurements. IEEE Transactions on Microwave Theory and Techniques, 2017, 65(5): 1488–1497

[20] Zhao J, Jing X, Wang W, Tian Y, Zhu D, Shi G. Steady method to retrieve effective electromagnetic parameters of bianisotropic metamaterials at one incident direction in the terahertz region. Optics & Laser Technology, 2017, 95: 56–62

[21] Shaltout A, Shalaev V, Kildishev A. Homogenization of bi-anisotropic metasurfaces. Optics Express, 2013, 21(19): 21941–21950

[22] Albooyeh M, Tretyakov S, Simovski C. Electromagnetic characterization of bianisotropic metasurfaces on refractive substrates: General theoretical framework. Annalen der Physik, 2016, 528(9-10): 721–737

[23] Lindell I V, Sihvola A H, Viitanen A J, Tretyakov S A. Electromagnetic Waves in Chiral and Bi-Isotropic Media. Boston: Artech House on Demand, 1994

[24] Serdiukov A, Semchenk I, Tertyakov S, Sihvola A. Electromagnetics of Bi-Anisotropic Materials-Theory and Application. Singapore: Gordon and Breach, 2001

[25] Sersic I, Tuambilangana C, Kampfrath T, Koenderink A F. Magnetoelectric point scattering theory for metamaterial scatterers. Physical Review B, 2011, 83(24): 245102

[26] Peng L, Zheng X, Wang K, Sang S, Chen Y, Wang Y G. Layer-by-layer design of bianisotropic metamaterial and its homogenization. Progress In Electromagnetics Research., 2017, 159: 39–47

[27] Xu J, Wu B, Chen Y. Elimination of polarization degeneracy in circularly symmetric bianisotropic waveguides: a decoupled case. Optics Express, 2015, 23(9): 11566–11575

[28] Multiphysics COMSOL. 5.2: a finite element analysis, solver and simulation software.

[29] Dong J W, Chen X D, Zhu H, Wang Y, Zhang X. Valley photonic crystals for control of spin and topology. Nature Materials, 2017, 16(3): 298–302

RIGHTS & PERMISSIONS

Higher Education Press