Discontinuities in open photonic waveguides: Rigorous 3D modeling with the finite element method

TL;DR

This paper introduces a comprehensive 3D finite element method approach for analyzing discontinuities in open photonic waveguides, including complex structures like plasmonic inclusions, with applications in infrared sensing.

Contribution

It presents a rigorous methodology for modeling arbitrary discontinuities in open waveguides using full vector Maxwell's equations and finite element analysis.

Findings

Computed leaky modes and their scattering behavior.

Validated the approach with numerical examples in infrared sensing.

Compared modal analysis with transmission in structured waveguides.

Abstract

In this paper, a general methodology to study rigorously discontinuities in open waveguides is presented. It relies on a full vector description given by Maxwell's equations in the framework of the finite element method. The discontinuities are not necessarily small perturbations of the initial waveguide and can be very general, such as plasmonic inclusions of arbitrary shapes. The leaky modes of the invariant structure are first computed and then injected as incident fields in the full structure with obstacles using a scattered field approach. The resulting scattered field is finally projected on the modes of the invariant structure making use of their bi-orthogonality. The energy balance is discussed. Finally, the modes of open waveguides periodically structured along the propagation direction are computed. The relevant complex propagation constants are compared to the transmission obtained for a finite number of identical cells. The relevance and complementarity of the two approaches are highlighted on a numerical example encountered in infrared sensing. Open source models allowing to retrieve most of the results of this paper are provided.

Equations

$$\varepsilon_{r,2D}(x,y)=\left\{\begin{array}[]{l}\varepsilon_{r,g}\mbox{ in the guide,}\\ \varepsilon_{r,l}\mbox{ in the low index region,}\\ \varepsilon_{r,s}\mbox{ in the substrate,}\\ \varepsilon_{r,t}\mbox{ in the superstrate}\end{array}\right..$$

$$\varepsilon_{r,3D}(x,y,z)=\left\{\begin{array}[]{l}\varepsilon_{r,g}\mbox{ in the guide,}\\ \varepsilon_{r,l}\mbox{ in the low index region,}\\ \varepsilon_{r,s}\mbox{ in the substrate,}\\ \varepsilon_{r,t}\mbox{ in the superstrate},\\ \varepsilon_{r,d}(x,y,z){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mbox{ in the obstacles }}\\ \hskip 14.22636pt{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mbox{ of the modified segment}}\end{array}\right..$$

$$\mathbf{curl}\left({\mathbf{\mu }}_{r,2D}^{-1}\mathbf{curl}\mathbf{E}\right)={\mathbf{\epsilon }}_{r,2D}{\left(\frac{{\omega }_0}{c}\right)}^2\mathbf{E}$$

$$\begin{array}{rl}& \mathbf{curl}\left({\mathbf{\mu }}_{r,2D}^{-1}\mathbf{curl}{\mathbf{e}}_k\right)-k_0^2{\mathbf{\epsilon }}_{r,2D}{\mathbf{e}}_k\\ +i{\beta }_k& \left[\hat{\mathbf{z}}\times \left({\mathbf{\mu }}_{r,2D}^{-1}\mathbf{curl}{\mathbf{e}}_k\right)+\mathbf{curl}\left({\mathbf{\mu }}_{r,2D}^{-1}\hat{\mathbf{z}}\times {\mathbf{e}}_k\right)\right]\\ +(i{\beta }_k{)}^2& \hat{\mathbf{z}}\times \left({\mathbf{\mu }}_{r,2D}^{-1}\hat{\mathbf{z}}\times {\mathbf{e}}_k\right)=0,\end{array}$$

$$-\mathbf{curl}\left[{\mathbf{\mu }}_{r,3D}^{-1}\mathbf{curl}{\mathbf{E}}^{\mathrm{tot}}\right]+k_0^2{\mathbf{\epsilon }}_{r,3D}{\mathbf{E}}^{\mathrm{tot}}=0.$$

$$-\mathbf{curl}\left[{\mathbf{\mu }}_{r,3D}^{-1}\mathbf{curl}{\mathbf{E}}^{\mathrm{d}}\right]+k_0^2{\mathbf{\epsilon }}_{r,3D}{\mathbf{E}}^{\mathrm{d}}=k_0^2({\mathbf{\epsilon }}_{r,2D}-{\mathbf{\epsilon }}_{r,3D}){\mathbf{E}}^{\mathrm{inc}}.$$

$$\begin{split}P_{in}&=\int_{\Gamma_{in}}\mathrm{\mathbf{S}}^{\mathrm{inc}}\cdot\mathbf{n}_{\Gamma}\,\mathrm{d}S,\\ P_{tr}&=\int_{\Gamma_{out}}\mathrm{\mathbf{S}}^{\mathrm{tot}}\cdot\mathbf{n}_{\Gamma}\,\mathrm{d}S,\\ P_{ref}&=\int_{\Gamma_{in}}\mathbf{S}^{d}\cdot\mathbf{n}_{\Gamma}\,\mathrm{d}S\mbox{ and }\\ P_{abs}&=\frac{\varepsilon_{0}\,\omega_{0}}{2}\int_{\Omega_{d}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathrm{Im}\{\varepsilon_{r,d}\}}\,|\mathrm{\mathbf{E}}^{\mathrm{tot}}|^{2}\,\mathrm{d}\Omega,\end{split}$$

$$T=\frac{P_{tr}}{P_{in}}\text{mbox},R=-\frac{P_{ref}}{P_{in}}\text{mbox}andA=\frac{P_{abs}}{P_{in}}.$$

$${\int }_S{\mathbf{e}}_j\times {\mathbf{h}}_k\cdot \hat{\mathbf{z}}\mathrm{d}S={\int }_S{\mathbf{e}}_k\times {\mathbf{h}}_j\cdot \hat{\mathbf{z}}\mathrm{d}S=A_k{\delta }_{kj}\text{mbox},where$$

$$\left\{\begin{array}[]{l}t_{k}=\displaystyle\int_{S_{out}}\mathrm{\mathbf{E}}^{\mathrm{tot}}\times\mathbf{h}_{k}\cdot\hat{\mathbf{z}}\,\mathrm{d}S/A_{k}\\[11.38109pt] r_{k}=\displaystyle\int_{S_{in}}\mathrm{\mathbf{E}}^{\mathrm{d}}\times\mathbf{h}_{k}\cdot\hat{\mathbf{z}}\,\mathrm{d}S/A_{k}\end{array}\right.$$

$$\mathbf{curl}\left[{\mathbf{\mu }}_{r,\#}(\mathbf{x}{)}^{-1}\mathbf{curl}\mathbf{E}\right]={\mathbf{\epsilon }}_{r,\#}(\mathbf{x},{\omega }_0){\left(\frac{{\omega }_0}{c}\right)}^2\mathbf{E}.$$

$$\mathbf{E}={\mathbf{E}}_{\#}(x,y,z)e^{-i({\omega }_{0}t-\gamma z)},$$

$$-\mathbf{curl}\left[{\mathbf{\mu }}_{r,\#}^{-1}\mathbf{curl}\left[{\mathbf{E}}_{\#}{e}^{i\gamma z}\right]\right]+k_0^2{\mathbf{\epsilon }}_{r,\#}(\mathbf{x}){\mathbf{E}}_{\#}{e}^{i\gamma z}$$

$$\mathrm{div}\left[{\mathbf{\epsilon }}_{r,\#}(\mathbf{x}){\mathbf{E}}_{\#}{e}^{i\gamma z}\right]$$

$$\begin{array}{rl}-\mathbf{curl}\left[{\mathbf{\mu }}_{r,\#}^{-1}\mathbf{curl}{\mathbf{E}}_{\#}\right]& +k_0^2{\mathbf{\epsilon }}_{r,\#}(\mathbf{x}){\mathbf{E}}_{\#}\\ -i\gamma & \hat{\mathbf{z}}\times \left[{\mathbf{\mu }}_{r,\#}^{-1}\mathbf{curl}{\mathbf{E}}_{\#}\right]\\ -i\gamma & \mathbf{curl}\left({\mathbf{\mu }}_{r,\#}^{-1}\hat{\mathbf{z}}\times {\mathbf{E}}_{\#}\right)\\ -(i\gamma {)}^2& \hat{\mathbf{z}}\times \left({\mathbf{\mu }}_{r,\#}^{-1}\hat{\mathbf{z}}\times {\mathbf{E}}_{\#}\right)=0.\end{array}$$

$$\begin{array}{rl}-& {\int }_{\Omega }\left[{\mathbf{\mu }}_{r,\#}^{-1}\mathbf{curl}{\mathbf{E}}_{\#}\right]\cdot \mathbf{curl}\overline{\mathbf{W}}\mathrm{d}\Omega \\ +& {\int }_{\Omega }{k}_0^2{\mathbf{\epsilon }}_{r,\#}(\mathbf{x}){\mathbf{E}}_{\#}\cdot \overline{\mathbf{W}}\mathrm{d}\Omega \\ -i\gamma & {\int }_{\Omega }\hat{\mathbf{z}}\times \left[{\mathbf{\mu }}_{r,\#}^{-1}\mathbf{curl}\mathbf{E}\right]\cdot \overline{\mathbf{W}}\mathrm{d}\Omega \\ -i\gamma & {\int }_{\Omega }\left({\mathbf{\mu }}_{r,\#}^{-1}\hat{\mathbf{z}}\times {\mathbf{E}}_{\#}\right)\cdot \mathbf{curl}\overline{\mathbf{W}}\mathrm{d}\Omega \\ +(i\gamma {)}^2& {\int }_{\Omega }\left({\mathbf{\mu }}_{r,\#}^{-1}\hat{\mathbf{z}}\times {\mathbf{E}}_{\#}\right)\cdot (\hat{\mathbf{z}}\times \overline{\mathbf{W}})\mathrm{d}\Omega \\ -& {\int }_{\partial \Omega }\left[{\mathbf{n}}_{|\partial \Omega }\times \left({\mathbf{\mu }}_{r,\#}^{-1}\mathbf{curl}{\mathbf{E}}_{\#}\right)\right]\cdot \overline{\mathbf{W}}\mathrm{d}S\\ -i\gamma & {\int }_{\partial \Omega }\left[{\mathbf{n}}_{|\partial \Omega }\times \left({\mathbf{\mu }}_{r,\#}^{-1}\hat{\mathbf{z}}\times {\mathbf{E}}_{\#}\right)\right]\cdot \overline{\mathbf{W}}\mathrm{d}S\\ & =0\end{array}$$

$$\mathrm{div}\left({\mathbf{\epsilon }}_{r,\#}{\mathbf{E}}_{\#}{e}^{i\gamma z}\right)=0=\mathrm{div}({\mathbf{\epsilon }}_{r,\#}{\mathbf{E}}_{\#})+i\gamma \hat{\mathbf{z}}\cdot ({\mathbf{\epsilon }}_{r,\#}{\mathbf{E}}_{\#})$$

$$\begin{array}{rl}{\int }_{\Omega }& \left[\mathrm{div}\left({\mathbf{\epsilon }}_{r,\#}{\mathbf{E}}_{\#}\right)+i\gamma \hat{\mathbf{z}}\cdot \left({\mathbf{\epsilon }}_{r,\#}{\mathbf{E}}_{\#}\right)\right]\overline{\phi }\mathrm{d}\Omega =0\\ =-{\int }_{\Omega }& {\mathbf{\epsilon }}_{r,\#}{\mathbf{E}}_{\#}\cdot \overline{\mathbf{grad}\phi }\mathrm{d}\Omega +i\gamma {\int }_{\Omega }\hat{\mathbf{z}}\cdot \left({\mathbf{\epsilon }}_{r,\#}{\mathbf{E}}_{\#}\right)\overline{\phi }\mathrm{d}\Omega ,\end{array}$$

$$\begin{array}{rl}-& {\int }_{\Omega }{\mathbf{\mu }}_{r,\#}^{-1}\mathbf{curl}{\mathbf{E}}_{\#,k}\cdot \mathbf{curl}\overline{\mathbf{W}}\mathrm{d}\Omega \\ +& {\int }_{\Omega }{k}_0^2{\mathbf{\epsilon }}_{r,\#}(\mathbf{x}){\mathbf{E}}_{\#,k}\cdot \overline{\mathbf{W}}\mathrm{d}\Omega \\ -i\gamma & {\int }_{\Omega }\hat{\mathbf{z}}\times \left({\mathbf{\mu }}_{r,\#}^{-1}\mathbf{curl}\mathbf{E}\right)\cdot \overline{\mathbf{W}}\mathrm{d}\Omega \\ -i\gamma & {\int }_{\Omega }\left({\mathbf{\mu }}_{r,\#}^{-1}\hat{\mathbf{z}}\times {\mathbf{E}}_{\#,k}\right)\cdot \mathbf{curl}\overline{\mathbf{W}}\mathrm{d}\Omega \\ +(i\gamma {)}^2& {\int }_{\Omega }\left({\mathbf{\mu }}_{r,\#}^{-1}\hat{\mathbf{z}}\times {\mathbf{E}}_{\#,k}\right)\cdot (\hat{\mathbf{z}}\times \overline{\mathbf{W}})\mathrm{d}\Omega \\ +& {\int }_{\Omega }{\mathbf{\epsilon }}_{r,\#}\mathbf{grad}{\phi }_k\cdot \overline{\mathbf{W}}\mathrm{d}\Omega \\ +i\gamma & {\int }_{\Omega }{\mathbf{\epsilon }}_{r,\#}{\phi }_k\hat{\mathbf{z}}\cdot \overline{\mathbf{W}}\mathrm{d}\Omega =0\end{array}$$

$$\begin{array}{rl}& {\int }_{\Omega }{\mathbf{\epsilon }}_{r,\#}{\mathbf{E}}_{\#,k}\cdot \overline{\mathbf{grad}{\phi }_k}\mathrm{d}\Omega \\ -i\gamma & {\int }_{\Omega }\hat{\mathbf{z}}\cdot \left({\mathbf{\epsilon }}_{r,\#}{\mathbf{E}}_{\#,k}\right)\overline{{\phi }_k}\mathrm{d}\Omega =0\end{array}$$

Full text

1 Introduction

The study of discontinuities is an old research topic in waveguide studies due to its importance for practical applications in many areas of physics. One must cite the seminal contribution of Schwinger for the development of variational methods in the forties [1] and the results obtained by Lewin [2].

These methods, often complex and specific, do not generally consider the exact solutions of Maxwell’s equations and rely on specific configurations, hypotheses, initial guesses for the solution forms. During the last two decades, the versatile Finite-Difference Time-Domain (FDTD) method allowed the study of waveguide discontinuities, including 3D ones, taking into account the full set of Maxwell’s equations [3, 4]. Nevertheless, the computational resources both in terms of memory and time requirements are huge when realistic 3D photonic devices are considered with a uniform square grid, especially nanophotonic ones with high quality factors. As for harmonic methods, Fourier modal methods [5, 6] also allow to tackle discontinuities in waveguides but they are restricted to geometries with straight walls. In acoustics and optics, coupled modal-finite element techniques have been successfully used in varying cross-section waveguides [7, 8].

Two types of non-Hermitian eigenvalue problems arise in open nanophotonic structures. When considering resonators in the general case, no particular ansatz can be guessed for the electromagnetic field and a natural eigenvalue is the (complex) frequency. Diffraction gratings are a special yet frequent case where the Bloch theorem applies ; the Bloch variable arises as a wavenumber and a corner stone of the dispersion relation of gratings consists in looking for complex frequencies for a given real wavenumber. Recent benchmarks of the numerical methods cited in the previous paragraph can be found in Refs. [9, 10] for this type of non-Hermitian eigenvalue problems. But when considering guiding structures as it is the case in this article, a more natural eigenvalue is the propagation constant at a given real frequency. When coupling a laser – indeed operating at a real frequency – into a waveguide, the relavant quantities are the light velocity and attenuation which are directly related to the real and imaginary parts of the eigenvalues (i.e. complex propagation constants) of the modes of the leaky waveguide.

This article addresses the numerical characterization of open waveguides with discontinuities. We demonstrate that it can be carried out efficiently with adapted formulations of the finite element method (FEM) which has already proven its efficiency and versatility in many field of computational electrodynamics [11].

We can state four main advantages of the FEM-based method: i) curved geometries are naturally treated using high order mesh elements and corresponding shape functions, ii) conforming non-uniform meshing is now a standard for mesh generators which is particularly relevant when rapid and strong permittivity changes must be tackled, iii) the domain decomposition method, now available in several FEM solvers, allows the treatment of large scale 3D problem, iv) and the possibility to reuse the inverse matrix for several incident modes propagating in the invariant structure – a subtlety detailed later which is a key advantage for the optimization multi-mode guides. This is especially worthy when the simulations are performed within a topology optimization frame [12].

The practical context motivating this theoretical and numerical study is the design of efficient plasmonic waveguides for infrared sensing [13] since the mid-IR spectral domain is known to be the molecular fingerprint region, due to the fact that most molecule including pollutants have intense fundamental vibrational bands in this spectral range. The device configuration is fully integrated and based on a ridge waveguide upon which metallic scattering nano-objects will ensure the coupling between the guided modes and superstrate of the device. Chalcogenide glasses are chosen for the main layers due to their high transparencies for infrared wavelengths [14, 15]. Ultimately, the metallic scatterers are planned to be functionalized in order to react to the targeted chemical species. The sensing property relies on the subsequent modification of the guidance of the full structure.

With this application in mind, we present a general framework to study rigorously discontinuous waveguides using a full vector description given by Maxwell’s equations in the framework of the finite element method. The discontinuity can be very general and is not necessarily a small perturbation of the initial waveguide. The full structure under investigation is made of 3 segments: The input one is a uniform waveguide invariant along its main propagation axis, the intermediate one (called “modified segment” in Fig. 1) contains the opto-geometrical modifications of the waveguide, the output one is again an invariant waveguide. In order to model the response of the resulting 3D guiding structure, we adopt a scattered field formulation consisting of three sequential steps, the output of first step being the input of the second one, the output of the first and second being the input of the third one.

First, we determine the leaky modes of the unperturbed 2D waveguide for a fixed frequency corresponding to the freespace wavelength of interest. The example chosen to illustrate our method consists in a ridge waveguide made of chalcogenide layers on a silicon substrate, assumed to be invariant along its propagation axis. We use our usual vector FEM method with the Galerkin approach to solve the relevant eigenvalue problem [16, 11, 17]. This first step provides both the propagation constants (eigenvalues) and the associated modes profiles (eigenvectors).

Second, these guided modes are used as incident fields for the full 3D problem in the modified segment. The electromagnetic problem to solve for this second step is then a mere scattering problem [18]. It is possible to define a proper energy balance (transmission and reflection, absorption taking place into the obstacles, extra radiation losses) allowing to fully evaluate the impact of the modified segment on the energy propagation.

Third, outside of the modified region, the total field is expanded on a fixed number of leaky modes of the output segment of the full structure. A special attention is paid to the coupling efficiency into the mode initially injected after crossing the modified segment. Our method allows to compute all the required energy-related quantities to investigate quantitatively the behavior of the full structure, notably the impact of the modified segment, and to take into account the way it is excited by the selected input propagating mode.

Note that our approach differs from the one exposed in Ref. [11] where total field formulations making use of port boundary conditions are applied to closed discontinuous waveguides, whereas it is proposed here to use a general scattered field formulation to deal with open discontinuous waveguides.

Finally, we also compute the modes of waveguides infinitely periodically structured along the propagation axis and compare the relevant complex propagation constants to the transmission obtained with a finite number of identical cells. After deriving the formulation, the relevance and complementarity of the two approaches are highlighted on a numerical example.

2 Direct problem

In this section, a direct – as opposed to modal – scattering approach is introduced. A typical and realistic structure is sketched in Fig. 1: A $z$-invariant dielectric rectangular waveguide (of width $w$ and thickness $h_{g}$ in blue) is deposited on a low index spacer (of thickness $h_{l}$, in green) lying on a semi-infinite substrate (in purple). The $z$-invariance of this guiding structure is locally broken, by adjunction of a finite number of obstacles. These obstacles can be in practice any bounded modification of permittivity: Ellipsoidal patches above the guiding layer labelled ➊ in Fig. 1, holes ➋ in the guiding layer, obstacles or resonators next to the waveguide ➌ or even a combination of all ➍…Note that the method applies irrespectively of the number of layers of the $z$-invariant structure and that the obstacles can be arbitrarily shaped and located in (or above) the structure. It is shown how the obstacles (more generally the modified waveguide segment) perturb a mode propagating in the $z$-invariant structure. A first step consists in the numerical computation of the modes of the invariant structure, which are used in a second step as incident fields for the full 3D structure.

2.A Obtaining the incident fields

The classical guiding $z$-invariant structure is characterized by its permittivity function defined by parts as:

[TABLE]

Between $z_{min}$ and $z_{max}$, one can now break the $z$-invariance by a local modification of the permittivity function which leads to a 3D scattering problem, which in turn can be characterized by its permittivity function defined by parts:

[TABLE]

The starting point consists in computing the modes of an annex problem formed by the $z$-invariant structure solely. This is a very classical problem [19, 17] where one introduces the ansatz $\mathbf{E}(\mathbf{x})=\mathbf{e}(\mathbf{x})e^{-i(\omega_{0}t-\beta z)}$ in the source-free Helmholtz equation:

[TABLE]

for a given real angular frequency $\omega_{0}$. Note that the relative permittivity and permeability in (3) are tensors fields since cartesian PMLs adapted to each domain with infinite extension (superstrate, low index region and substrate) are used to damp the radially blowing leaky modes of this open structure [17]. The domains $\Omega$ involved in all the formulations of the paper are correspond to geometrical domains surrounded by appropriate PMLs of finite thicknesses.

It results in a quadratic non-Hermitian eigenvalue problem amounting to find non trivial pairs $(\beta_{k},\mathbf{e}_{k})\in\mathbb{C}\times H^{1}(\Omega_{2D},\mathbf{curl})$ such that :

[TABLE]

with $k_{0}:=\omega_{0}/c$.

This equation can be solved using a mixed finite element formulation involving edge elements for the discretization of the transverse component $(e_{x},e_{y})$ coupled to a nodal basis for the (continuous) longitudinal component $e_{z}$. The rather lengthy details of the resulting weak formulation can be found in Refs. [20, 17].

Throughout the paper, the following numerical values are considered for the $z$-invariant waveguide [21]: The operating freespace wavelength $\lambda_{0}=$7.7\text{,}\mathrm{\SIUnitSymbolMicro m}, $\varepsilon_{r,g}=7.1824$ ($\mbox{Se}_{4}$), $\varepsilon_{r,l}=6.2001$ ($\mbox{Se}_{2}$), $\varepsilon_{r,s}=11.69024481$ (silicon) [[22](#bib.bib22)], $\varepsilon_{r,t}=1$ (air), $w_{g}=$14\text{\,}\mathrm{\SIUnitSymbolMicro m}, $h_{g}=$2.2\text{,}\mathrm{\SIUnitSymbolMicro m} and $h_{l}=$5.3\text{\,}\mathrm{\SIUnitSymbolMicro m}. All the materials are assumed to be non-magnetic: $\boldsymbol{\mu}_{r,2D}=\underline{\underline{I}}$ (except in the PMLs where $\boldsymbol{\mu}_{r,2D}$ takes the appropriate value), where $\underline{\underline{I}}$ is the $3\times 3$ identity tensor.

The modes of this structure associated with eigenvalues with lowest imaginary parts are depicted in Fig. 2(a-f), sorted in ascending attenuation (i.e. $\mbox{Im}\{\beta_{1}\}$ is the smallest). The black and white colormaps show the norm of the electric eigenfields $|\mathbf{e}_{k}|$ in Fig. 2(a-f) and the orange arrows represent the real part of $\mathbf{e}_{k}$, allowing to distinguish a TE-like mode from a TM-like one. These six modes are also those with electric eigenfield most confined into the core region of the structure (ridge).

Finally, all geometries and conformal meshes have been obtained using the Gmsh software [23] and all the finite element formulations in this article are implemented thanks to the flexibility of the finite element software GetDP [24]. Open source models allowing to retrieve most of the results of this article are provided [25].

2.B Computation of the scattered field

One can now use any of these 2D modes $\mathbf{E}_{k,2D}:=\mathbf{e}_{k}\,e^{i(\beta_{k}z-\omega_{0}t)}$ as an incident field $\mathrm{\mathbf{E}}^{\mathrm{inc}}$ on the obstacles and look for $\mathrm{\mathbf{E}}^{\mathrm{tot}}$, the total field solution of the source-free Helmholtz equation:

[TABLE]

Let us define the scattered field as $\mathrm{\mathbf{E}}^{\mathrm{d}}\equiv\mathrm{\mathbf{E}}^{\mathrm{tot}}-\mathrm{\mathbf{E}}^{\mathrm{inc}}$ and from the linearity of Eqs. (3,5), we obtain the following scattering problem:

[TABLE]

Note that the support of the effective sources $(\boldsymbol{\varepsilon}_{r,2D}-\boldsymbol{\varepsilon}_{r,3D})$ in this scattering problem has to be bounded to ensure a proper outgoing wave condition [26] to the scattered field $\mathrm{\mathbf{E}}^{\mathrm{d}}$, which is the case in our examples. Finally, 3D cartesian PMLs are used to bound the computational domain [17, 27] as shown in grey lines in Fig. 3. Compared to a total field approach with a port condition [11], it is stressed that the electromagnetic sources of our equivalent radiation problem are located within the discontinuities. The PMLs of elongated structures are naturally built to damp fields radiating from the center of the computational box more efficiently than the total field radiating from a port located at one extremity of the elongated box, the resulting total field being more grazing than the scattered field when entering the PMLs.

2.C Energy balance

The Poynting vectors associated with the incident, diffracted and total fields are classically defined by respectively $\mathrm{\mathbf{S}}^{\mathrm{inc}}=\mathrm{\mathbf{E}}^{\mathrm{inc}}\times\overline{{\mathrm{\mathbf{H}}^{\mathrm{inc}}}}/2$, $\mathrm{\mathbf{S}}^{\mathrm{d}}=\mathrm{\mathbf{E}}^{\mathrm{d}}\times\overline{\mathrm{\mathbf{H}}^{\mathrm{d}}}/2$ and $\mathrm{\mathbf{S}}^{\mathrm{tot}}=\mathrm{\mathbf{E}}^{\mathrm{tot}}\times\overline{\mathrm{\mathbf{H}}^{\mathrm{tot}}}/2$, where the horizontal bar means complex conjugation. Then, the incoming, transmitted, reflected and absorbed powers can be defined as respectively :

[TABLE]

where $\Gamma_{in}$ and $\Gamma_{out}$ are the transverse plane surfaces before and after the obstacles depicted in transparent grey color in Fig. 3, $\mathbf{n}_{\Gamma}$ is the unit vector normal to $\Gamma_{in}$ and $\Gamma_{out}$ and $\Omega_{d}$ is the support of the diffractive obstacles or of the localized region where the waveguide opto-geometrical parameters are modified. Finally one can define transmission (T), reflection (R) and absorption (A) coefficients as :

[TABLE]

For clarity, Fig. 3 illustrates the quantities at stake in the energy balance. This numerical set up is obtained for an incident field set to $\mathbf{E}_{1,2D}$ (cf Fig. 2(a)) with four ellipsoidal lossy patches placed above the same waveguide as in Section 2.A. This particular configuration will be discussed in detail in Section 4. The red/blue colormap represents the norm of the total electric field $|\mathrm{\mathbf{E}}^{\mathrm{tot}}|$ involved in the computation of the Joule losses, i.e. $P_{abs}$ in (7). The purple/yellow colormap represents three cuts of the $z$ component of the total Poynting vector $\mathrm{\mathbf{S}}^{\mathrm{tot}}$ on three selected plane surfaces. The first cut is taken at $z=z_{min}$ (see left side of the figure), another one at $z=z_{max}$ (right side of the figure) and the last one along the symmetry plane $x=0$ of the structure. In this last cut, it is clear that the perturbation induced by the objects affects the $z$ component of the total Poynting vector since in absence of scattering objects above the ridge, this map would be constant along $z$.

In this example, the transmission $T$ reaches 0.688, the reflection $R=0.007$ and the absorption $A=0.224$ ($T+R+A=0.904$). Note that $R$, $T$ and $A$ are defined here in order to match commonly measured quantities does but do not add up to unity. It is nonetheless expected since their sum does not correspond to a full Poynting balance. It can be easily completed by adding the flux contributions from the surfaces parallel to the $zOx$ and $zOy$ planes, which represent the extra radiative leakage induced by the modified segment ($9.6\%$ in the present case).

As shown in Sec. 2.E, the scattering process can be further precised by expansion of the diffracted and total fields outside the modified segment on the modes of the 2D invariant structure.

2.D Discretization and convergence

In the 2D eigenvalue problem of Sec. 2.A, the longitudinal component of the electric field, which is a continuous scalar field, is discretized using classical $P_{2}$ nodal elements having one Degree Of Freedom (DOF) per node and one DOF per edge. The transverse components are discretized using edge elements of the second order (2 degrees of freedom per edge). For eigenvalue problems, the GetDP software relies on the high performance library SLEPc [28] which implements advanced Krylov subspace methods for computing a small amount of eigenvalues of the large sparse matrices.

The 3D scattering problem uses high order Webb hierarchical edge elements [29, 30, 11] with 26 DOFs per tetrahedron (3 DOFs per edge, 2 DOFs per face). The direct problem described in section 2.B is solved using the direct solver MUMPS [31] interfaced in GetDP.

The convergence of the absorption $A$ as a function of the mesh refinement is shown in Fig. 4(a). The mesh size is parametrized by $n$ (in abscissa) and decreases as $\lambda_{0}/(n\,\mathrm{Re}\{\sqrt{\varepsilon_{r}}\})$ as $n$ increases. In other words, $n$ represents the average number of tetrahedrons per wavelength inside a dielectric material of relative permittivity $\varepsilon_{r}$. Note that in metals, the relevant physical length to consider for a proper spatial sampling of the field would be the skin depth rather that the wavelength. For $n$=1, that is roughly one tetrahedron per wavelength, the computational box in Fig. 3 leads to 32000 DOFs solved in $3\text{\,}\mathrm{s}$ on a laptop equipped with 4 cores 16 Gb of RAM memory. In this case, the local values of the field are poorly approximated, but the order of magnitude of integral quantities such as the absorption is relevant, as can be noticed on the left side of the convergence plot shown in Fig. 4(a). For $n$=4, the number of DOFs is 650000, the model still runs on the same laptop within $4\text{\,}\mathrm{min}$. For $n$=7, the number of DOFs becomes about 3 millions and a workstation equipped with 24 cores and 256 Gb RAM memory was used for a runtime of $30\text{\,}\mathrm{min}$. Five significant digits are then obtained on energy-related quantities such as the absorption.

2.E Modal expansion of the scattered field

The modes of the 2D invariant structure satisfy the following bi-orthogonality condition equivalent to the one given in Ref. [32, 19] which provides the normalization of each leaky modes:

[TABLE]

$S$ is an infinite cross-section of the open waveguide. In the case of leaky modes, it is suggested in [32] to perform a complex change of space variable as one moves far away from the waveguide to damp the exponential growth of the leaky mode. In our finite element approach that includes the PMLs which are an analytical continuation of the space variables, this integration simply corresponds to integrate over a full cross-section of computational domain including the PML regions. Therefore it is stressed that the cross-sections considered hereafter include the PMLs. Hence, away from the obstacles, it is possible to expand the scattered field as $\mathrm{\mathbf{E}}^{\mathrm{d}}=\sum_{k}r_{k}\mathbf{E}_{k,2D}$ and the total field as $\mathrm{\mathbf{E}}^{\mathrm{tot}}=\sum_{k}t_{k}\mathbf{E}_{k,2D}$ where the reflection and transmission coefficients are simply given as:

[TABLE]

where $S_{in}$ can be any transverse section before the obstacles ($z<z_{min}$) and $S_{out}$ can be any transverse section after the obstacles ($z>z_{out}$). Note that this formalism can be extended to the computation of the scattering matrix of the waveguide, by considering sequentially all the modes of the invariant structure computed in Sec. 2.A as incident fields. In this example, the modulus of the mean value of the off-diagonal coefficients of the $12\times 12$ bi-orthogonality matrix (see (9)) is less than $10^{-3}$ smaller than the modulus of the mean value of the diagonal coefficients $A_{k}$.

The total field within a leaky guiding structure can be expanded [19] as discrete sum over bounded modes plus an integral over the continuous spectrum. When using finite size PMLs, the continuous spectrum becomes discrete and the integral contribution turns into a discrete sum. Besides, the power orthogonality between the modes holding for self-adjoint eigenvalue problems such as the perfect metallic waveguide now fails in the our non-Hermitian case. This can be seen in the bi-orthogonality relations (see (9)) that involve the magnetic field rather than its complex conjugate. However, we are in a weakly leaky regime [32] were the weakly leaky modes decay rapidly as the radial distance $r=\sqrt{x^{2}+y^{2}}$ increases until experiencing an exponential blow at even larger radial distances. In this regime, it is interesting to see that the power exchange between the most highly confined modes (associated with eigenvalues with small imaginary part compared to their real part) can be neglected. In short, $\sum_{k}|t_{k}|^{2}\rightarrow T$ and $\sum_{k}|r_{k}|^{2}\rightarrow R$ hold to a very good approximation. In the particular configuration described in Fig. 3, one obtains the values $\sum_{k=1}^{12}|t_{k}|^{2}=0.685$ (with a major contribution from $|t_{1}|^{2}=0.632$, to be compared to $T$=0.688 obtained in Sec. 2.C) and $\sum_{k=1}^{12}|r_{k}|^{2}=0.0007$. The comparison between $T$ and $\sum_{k=1}^{M}|t_{k}|^{2}$ as a function of the truncation order $M$ is shown in Fig. 4(b) for a fine mesh with $n=7$.

3 Modes of the infinitely periodic 3D structure

3.A Variational formulation of the spectral problem

In this section, we are now interested in a 3D spectral problem with one direction of periodicity defined by

[TABLE]

where $\boldsymbol{\varepsilon}_{r,\#}(\mathbf{x},\omega_{0})$ and $\boldsymbol{\mu}_{r,\#}(\mathbf{x})$ are respectively the permittivity and permeability tensor fields at a fixed frequency $\omega_{0}$ exhibiting a 1D $d$-periodicity along $Oz$. Bloch’s theorem states that, without loss of generality, one can look for solutions for the electric field $\mathbf{E}$ under the form [33, 17] :

[TABLE]

where $\mathbf{E}_{\#}$ is a $d$-periodic function in $z$ and $\gamma$ is the Bloch variable lying in the first reduced Brillouin zone [0,$\pi/d$].

One can choose to set $\gamma$ to a real value lying in the first Brillouin zone and to look for ($\omega_{\gamma,i}$,$\mathbf{E}_{\gamma,i}$) eigenvalues and eigenvectors, by imposing Bloch conditions on the z-transverse surfaces of the cell and making the use of (11). An alternative option amounts to set $\omega_{0}$ to a real value, inject the ansatz in (12) into (11) and look for eigenvectors under the form of the periodic part $\mathrm{\mathbf{E}}_{\#}$ of the Bloch wave along with corresponding eigenvalue $\gamma$. In this latter case, two equations are to be fulfilled:

[TABLE]

It is stressed that the invariant 2D problem described in Sec. 2.A is a particular case of this 3D problem, an invariant structure along $z$ being trivially periodic in $z$ with arbitrary period. Unsurprisingly, expanding the $\mathbf{curl}$ term in (13a) in order to get rid of the $e^{i\gamma z}$ dependency, leads to an expression similar to the $z$-invariant counterpart of the problem (see (4)):

[TABLE]

In a variational way, after classically integrating by part two curl operators, it holds that for any $\mathbf{W}\in H_{\#}^{1}(\Omega,\mathbf{curl}\,)$ :

[TABLE]

Note that the two boundary terms recombine into $-\int_{\partial\Omega}[\mathbf{n}_{|\partial\Omega}\times(\boldsymbol{\mu}_{r,\#}^{-1}(\mathbf{curl}\,\mathrm{\mathbf{E}}_{\#}+i\gamma\hat{\mathbf{z}}\times\mathrm{\mathbf{E}}_{\#}))]\cdot\overline{\mathbf{W}}\,\mathrm{d}S\propto\int_{\partial\Omega}[\mathbf{n}_{|\partial\Omega}\times\mathbf{H}]\cdot\overline{\mathbf{W}}\,\mathrm{d}S$ so that setting a Dirichlet or Neumann natural condition for $\mathrm{\mathbf{E}}_{\#}$ on non-periodic faces of the domain (i.e. the PML bounds) actually corresponds to a Dirichlet or Neumann natural condition for $\mathbf{H}$.

The divergence condition in (13b) has to be handled carefully. Indeed, we are looking for divergence free solutions such that $\mathrm{div}\,(\boldsymbol{\varepsilon}_{r,\#}\,\mathbf{E})=0$, that is:

[TABLE]

Consequently, $\boldsymbol{\varepsilon}_{r,\#}\,\mathrm{\mathbf{E}}_{\#}$ is not divergence-free and, from the variational point of view, the following holds for any $\varphi\in H_{\#}^{1}(\Omega)$:

[TABLE]

where the boundary term arising from the integration by part vanishes due to periodicity and homogeneous conditions at the back of the PMLs.

Finally, the proper way to ensure the divergence condition [34] in a weak sense is to use $\varphi$ as a Lagrange multiplier. We are now in position to reformulate the eigenvalue problem at stake in this section. We are looking for non trivial pairs $\gamma_{k},(\mathrm{\mathbf{E}}_{\#,k},\varphi_{k})\in\mathbb{C}\times(H^{1}_{\#}(\Omega,\mathrm{curl})\times H_{\#}^{1}(\Omega))$ such that:

[TABLE]

3.B Discretization

The periodic vector unknown $\mathrm{\mathbf{E}}_{\#}$ is discretized using high order Webb hierarchical edge elements [29, 30] with 26 DOFs per tetrahedron (3 DOFs per edge, 2 DOFs per face). The scalar field $\varphi$ mapping the divergence is discretized using Lagrange $P_{3}$ elements, with 20 DOFs per tetrahedron (4 nodal DOFs, two DOFs per edge, one DOF per face). Periodic boundary conditions are imposed along the $z$ direction for both $\mathrm{\mathbf{E}}_{\#}$ and $\varphi$.

3.C Numerical validation

It is apropos to validate this 3D model numerically using an extruded 2D domain. The eigenvalue resulting from three finite element problems are shown in Fig. 5. Two of them (orange circles and blue crosses) are variants of the 2D problem in (4). The problem is indeed quadratic and can be solved as is using the SLEPc library (orange circles) which implements its own internal numerical linearization. But as detailed in Ref. [17], it is possible to linearize the 2D problem by simply using for unknown (e_{x},e_{y},i{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta}\,e_{z}) instead of $(e_{x},e_{y},e_{z})$. The resulting sparse systems resulting from the two methods are different and it is worth noting that the 30 eigenvalues are identical up to numerical precision. Now the third eigenproblem (red pluses in Fig. 5) corresponds to Eqs. (18a,18b) applied to the 3D problem obtained by simple extrusion along $z$ of the previous 2D problem by a period $d$=$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ along the $z$ direction. The value of $d$ can be arbitrarily chosen because of the translational invariance. The small value of the period along $z$ corresponds to a large first Brillouin zone so that the eigenvalues computed do not belong to a folded dispersion branch and the periodic part of the Bloch vector field is constant along $z$ as depicted in Figs. 5(a-c). Up to the $\pi$/d folding of the dispersion curves expected from the application to the Bloch theorem, these 2D and 3D invariant problems are spectrally equivalent and the eigenvalues are indeed retrieved with excellent accuracy. The discrepancies obtained for Im$\{\gamma_{k}/k_{0}\}>0.2$ can be simply explained by the fact that the 3D mesh used in the simulation is comparatively coarser than the 2D mesh. For the modes with smallest attenuation labelled (a), (b) and (c), the relative error between the 2D and 3D eigenvalues is lower that $10^{-4}$ in modulus. Note that convergence tests have been performed and this discrepancy decreases with the mesh size at the expected rate given the type and order of the chosen finite elements spaces. The real parts of $\mathbf{E}_{\#,k}$ are shown in the insets (a-c) of Fig. 5 (note that they cannot be directly compared to the yellow arrows of Fig. 2(a-c) since an eigenvector is defined up to an arbitrary complex number).

Computing 30 eigenvalues of the periodic 3D structure takes $80\text{\,}\mathrm{s}$ on a laptop with low order elements (twelve DOFs per tetrahedron for the vector unknown, 10 DOFs per tetrahedron for the scalar one, about 100000 DOFs in total) with the mesh paramater $n$=4. The same computation with $n$=8 and higher order edge elements (26 DOFs per tetrahedron for the vector unknown, twenty DOFs per tetrahedron for the scalar one, about 800000 DOFs in total) takes $30\text{\,}\mathrm{min}$ on the 24 cores workstation.

We are finally in position to compare the results derived from the modes of the infinitely structured waveguide (3D modal problem defined in Sec. 3) to the transmission properties of scattering problems with a finite number of periods (3D direct problems defined in Sec. 2).

4 Discussion

The direct 3D approach detailed in Sec. 2 and the 3D modal one based on $z$-periodicity presented in the previous section are now compared. In both cases, the unit cell of the waveguide guide contains one lossy ellipsoidal patch defined by its relative permittivity $\varepsilon_{r,d}=9+2i$, height $h=$1\text{,}\mathrm{\SIUnitSymbolMicro m}, transverse and longitudinal radii $r_{t}=$5\text{\,}\mathrm{\SIUnitSymbolMicro m} and $r_{l}=$1.5\text{,}\mathrm{\SIUnitSymbolMicro m}, and period $d=$4\text{\,}\mathrm{\SIUnitSymbolMicro m} (see ➊ in Fig. 1). In the finite size problem, the number of patches is $N$ and the incident field is the fundamental mode labelled 1 in Fig. 2(a).

Figure 6 shows the transmission $T_{N}$ (black curve, cf (8)) as a function of the number $N$ unit cells (or scatterers). The orange curve represents the coefficient $|t_{1}|^{2}$ (cf (10)), which is a good approximation of the fraction of incident energy carried into mode 1 and remaining in this channel after crossing the modified waveguide segment containing the scatterers. This last curve lies below the green curve that represents the sum of the amplitudes transmitted into all 6 channels (or modes) of the $z$-invariant waveguide shown in Fig. 2. This numerical set up corresponds exactly to the configuration depicted in Fig. 3 with $N=4$.

Finally, one can correlate these results to the infinitely periodic structure and superimpose the last red curve that represents $T_{1}\,e^{-2\gamma^{\prime\prime}(N-1)d}$ (where $\gamma^{\prime\prime}=\mathrm{Im}\{\gamma\}$), the spatial damping of power associated with the mode with smallest propagation losses found using the modal approach detailed in Sec. 3. Note that the normalization factor $T_{1}$ (transmission obtained for the direct 3D problem for one single obstacle) accounts for the fact the damping of the 3D mode computed with one obstacle does not make any sens in absence of obstacle (i.e. for $N=0$) and represents the input impedance of the structured waveguide. The norm of this 3D mode is represented in purple/yellow colors at the top right of the Fig. 6. The direction of the eigenfield is not represented for clarity, but it is globally polarized along $x$. This 3D mode is very similar to mode 1 in Fig. 2(a), which is the one injected in the direct problem to obtain the three other black, orange and green curves. The consistency between these four quantities is remarkable.

5 Conclusion

In conclusion, we present in this paper a general finite element frame for the study of discontinuous waveguides, from isolated discontinuities to fully periodic ones. A first method, adapted to a finite set of discontinuities allows to compute, given the modes of the invariant structure, the field scattered by the local discontinuity, all relevant energy related quantities, and the projection of the scattered field on the modes of the invariant structure (that is the elements of the transition scattering matrix).

When the modified region extends to infinity with periodic discontinuities, the relevant quantity is the dispersion relation of the so formed structured waveguide. An adequate weak treatment of divergence condition allows to determine these modes with accuracy.

The two numerical models presented in this paper show great interest for the design of structured waveguides. Note that the methodology adopted for the direct problem is very general and can readily be applied to a large variety of guiding structure and geometry of objects located in the modified waveguide segment.

The two methods are in fact complementary. To give a concrete example in sensing applications, the adjunction of well chosen periodic scatterers [35] above a waveguide allows for instance to strengthen the interaction of the light flowing in the superstrate near the ridge where e.g. the molecules to detect lie. In this case, it is enough to study the modes profiles of the infinitely periodic 3D structure as detailed in Sec. 3. It is indeed much faster than optimizing the scatterers properties using the direct problem introduced in Sec. 2 applied to a very long finite chain of scatterers. However, once the properties of the scatterers optimal for the targeted application, the practical device consists indeed of a finite chain. Then, the direct problem introduced in Sec. 2 is the ideal tool to study and optimize the coupling of an incident mode into the modified segment.

Finally, open source models allowing to retrieve most of the results of this paper are provided [25]. They can be tuned to handle different geometries and material properties.

This work will later be extended to the case where the input and output invariant structures mismatch using a coupled mode-FE approach [8] to compute the relevant incident fields.

Acknowledgement

This research was supported by ANR Louise project, grant ANR-15-CE04-000164 of the French Agence Nationale de la Recherche. The authors would like to thank the developers of MUMPS [31], PETSc [36], SLEPc [28], Gmsh [23] and GetDP [24] for maintaining and making their respective libraries freely available. Finally, the authors acknowledge Sonia Fliss (INRIA POEMS) for fruitful discussions.

Disclosures

The authors declare no conflicts of interest.