anisotropic waveguide

Guided modes of an anisotropic waveguide

Lithium niobate, liquid crystals and many other optical materials are anisotropic: the refractive index that light sees depends on its polarization. Build a waveguide out of such a material and the anisotropy splits modes that would otherwise travel at the same speed, and makes their effective index depend on how the crystal is cut. This notebook computes the guided modes of a square waveguide with an anisotropic core, using the finite element package gyptis, for four orientations of the crystal (one of them tilted out of the cross section), and checks the results against a classic benchmark (Koshiba, Hayata and Suzuki, J. Lightwave Technol. 4, 121 (1986)).

Along the way we meet the main ingredients of a modal simulation: the effective index, the guiding condition, vector edge elements, an eigenvalue problem for the propagation constant, and the polarization of the modes.

Guided modes and effective index

A waveguide is invariant along its axis $z$. We look for fields that keep their shape as they propagate,

with the $e^{j\omega t}$ time convention used by gyptis. The propagation constant $\beta$ tells how fast the phase advances along the guide. Comparing it with the vacuum wavenumber $k_0 = 2\pi/\lambda$ defines the effective index

the mode travels as a plane wave would in a uniform medium of index $n_\text{eff}$. Substituting this field into Maxwell’s equations (non-magnetic media, $\mu = 1$) gives

where $\nabla_t$ acts on $x$ and $y$ only. For a given frequency, this is an eigenvalue problem for $\beta$. Its solutions are a few discrete guided modes, plus a continuum of radiation modes that leak into the cladding. A mode is guided when

above $n_\text{clad}$ the field decays exponentially away from the core, so it carries its power along the guide without loss; and no mode can travel slower than light in the densest material it sees.

A uniaxial core

The core is a uniaxial crystal. One direction, the optic axis $\hat{\mathbf{c}}$, is special: an electric field along $\hat{\mathbf{c}}$ sees the extraordinary permittivity $\varepsilon_e$, a field perpendicular to it the ordinary permittivity $\varepsilon_o$. As a tensor,

A plane wave travelling perpendicular to $\hat{\mathbf{c}}$ therefore has two indices: $n_o = \sqrt{\varepsilon_o}$ when polarized perpendicular to the axis, $n_e = \sqrt{\varepsilon_e}$ when polarized along it. This is birefringence.

We take the values of the benchmark: a square core of side $t$ with $\varepsilon_o = 2.31$ and $\varepsilon_e = 2.19$, buried in an isotropic cladding with $\varepsilon = 2.05$. Only the orientation of the optic axis changes:

• along $z$, the axis of the guide;
• along $x$, parallel to one side of the core;
• diagonal, along $(\hat{\mathbf{x}} + \hat{\mathbf{y}})/\sqrt{2}$ in the cross section;
• tilted, along $(\hat{\mathbf{y}} + \hat{\mathbf{z}})/\sqrt{2}$, halfway between the cross section and the axis of the guide.

All lengths are in units of $t$, so the frequency enters through the dimensionless product $k_0 t$.

The four tensors share the same eigenvalues, only their principal axes turn. With the axis along $z$, the cross section is isotropic: both transverse components see $\varepsilon_o$. With the axis along $x$, the tensor is diagonal but $E_x$ and $E_y$ see different permittivities. With the diagonal axis, an off-diagonal term $\varepsilon_{xy}$ couples $E_x$ and $E_y$. With the tilted axis, the off-diagonal term $\varepsilon_{yz}$ couples a transverse component, $E_y$, to the longitudinal one, $E_z$: this case needs more care.

Finite element formulation

We split the field into its transverse part $\mathbf{E}_t = (E_x, E_y)$ and its longitudinal part $E_z$. They call for different finite elements.

• $\mathbf{E}_t$ uses edge (Nédélec) elements, whose degrees of freedom are tangential components along the triangle edges. They enforce continuity of the tangential field across element boundaries and nothing more, which is exactly what Maxwell’s equations demand at a material interface, where the normal component jumps with $\varepsilon$. Nodal elements, which force all components to be continuous, produce spurious non-physical modes in this problem.
• $E_z$ is continuous everywhere and uses ordinary nodal (Lagrange) elements.

Write the permittivity tensor in blocks: $\varepsilon_{tt}$, the $2\times2$ transverse block, $\varepsilon_{zz}$, and the coupling column $\boldsymbol{\varepsilon}_{tz} = (\varepsilon_{xz}, \varepsilon_{yz})$. With the scaled unknowns $\mathbf{e}_t = \beta\,\mathbf{E}_t$ and $e_z = -j E_z$ (Lee, Sun and Cendes, IEEE Trans. Microw. Theory Tech. 39, 1262 (1991)), the weak form is a polynomial in $\beta$: for all test functions $(\mathbf{v}_t, v_z)$,

with

Each form becomes a sparse matrix. When the coupling column vanishes, so does $c$, and the problem is a generalized eigenvalue problem linear in $\beta^2$. That is the case when the optic axis lies in the cross section or along the guide. A tilted axis keeps the term linear in $\beta$, and the problem is a quadratic eigenvalue problem, which SLEPc solves with a polynomial eigensolver. gyptis inspects the tensors and picks the right one. In both cases, an iterative (Krylov–Schur) eigensolver with a shift-and-invert transform finds the eigenvalues closest to a chosen target; we aim just above $n_o k_0$, the largest $\beta$ a guided mode can have.

Geometry and mesh

gyptis builds the cross section as a layered box surrounded by perfectly matched layers (PMLs). We give both layers the cladding permittivity, so the core sits in a uniform cladding. Guided modes decay away from the core and barely reach the PMLs; radiation modes run into them and get absorbed, which turns their eigenvalues complex. That is how we will tell them apart.

The fundamental modes vary slowly across the core, so a moderate mesh suffices: triangles of size $t/10$ in the core and $t/3$ in the cladding. We check this choice at the end.

Solving for the modes

guided_modes sets up the problem for one orientation of the optic axis and one frequency (with quadratic=None, gyptis decides whether the problem is linear or quadratic), solves it, keeps the modes with a real effective index above $n_\text{clad}$, and sorts them from the most to the least confined. polarization measures which way a mode’s transverse field points: the fraction of $\int|\mathbf{E}_t|^2$ carried by the component along an in-plane unit vector $\hat{\mathbf{u}}$, integrated over the physical domain (PMLs excluded). It assembles the overlaps $M_{ij} = \mathrm{Re}\int E_i E_j^*\,\mathrm{d}\Omega$ once, so that the fraction along any $\hat{\mathbf{u}}$ is $\hat{\mathbf{u}}^{\mathsf{T}} M \hat{\mathbf{u}} / \mathrm{tr}\,M$.

Optic axis along the guide

We start at $k_0 t = 10$ with the optic axis along $z$. The cross section is then isotropic for transverse fields, and the square core is symmetric under a quarter turn, which maps a mode polarized along $x$ onto one polarized along $y$. The two fundamental modes, traditionally called $E^x_{11}$ and $E^y_{11}$ after their dominant field component, must share the same effective index.

The fundamental pair is degenerate up to a splitting of order $10^{-5}$, which comes from the mesh: it is unstructured and does not share the exact symmetry of the square. Any combination of two degenerate modes is also a mode, so in principle the pair has no preferred polarization; here the tiny asymmetry of the mesh happens to select modes polarized almost purely along $x$ and $y$. The next four modes form the second-order family, split into close pairs in the same way.

Optic axis in the cross section

Now we turn the optic axis into the cross section, first along $x$, then along the diagonal.

The degeneracy is lifted: the two fundamental modes now differ in effective index by a few $10^{-2}$, the same order as the bulk birefringence $n_o - n_e$, though smaller, because part of each mode travels in the isotropic cladding. The two in-plane orientations give the same pair to within $2\times10^{-5}$. The plots below show why: in both cases one mode is polarized perpendicular to the optic axis and sees $\varepsilon_o$, the other along it and sees $\varepsilon_e$. The orientation of the axis relative to the sides of the core hardly matters, and symmetry explains it. Treat the anisotropy as a small perturbation of an isotropic core: it shifts the effective index of a mode polarized along $\hat{\mathbf{u}}$ in proportion to the overlap of the mode with the core, weighted by $\hat{\mathbf{u}}^{\mathsf{T}}\hat{\mathbf{c}}$. The quarter-turn symmetry of the square makes that overlap the same for every in-plane direction, so to first order the core cannot tell a side from a diagonal. Only the vector corrections near the corners of the core, a weaker effect, separate the two cases.

Optic axis tilted out of the cross section

Before solving, consider plane waves in the bulk tilted crystal. Polarized along $x$, perpendicular to the optic axis, a wave sees $\varepsilon_o$ whatever its direction in the $yz$ plane. Polarized in the $yz$ plane, it is more subtle. A wave travelling exactly along $z$ must have a transverse displacement, $D_z = \varepsilon_{zy}E_y + \varepsilon_{zz}E_z = 0$, so its electric field gains a longitudinal component $E_z = -(\varepsilon_{zy}/\varepsilon_{zz})\,E_y$ and it sees

the familiar extraordinary index at $45°$ from the optic axis. But its electric field is then no longer perpendicular to $z$, and its energy does not flow along $z$: it walks off. A guided mode, on the contrary, must carry its energy along the guide. The plane wave that does so has a wavevector tilted in the $yz$ plane, with a phase varying as $e^{-j(py + \beta z)}$, $p = -\beta\,\varepsilon_{yz}/\varepsilon_{yy}$, and its propagation constant is $\beta = k_0\sqrt{\varepsilon_{yy}}$, the largest of all directions. So the bulk reference for the guide is $\varepsilon_{yy}$, not $\varepsilon_y$. The coupling term $c$ of the weak form is what lets a mode tilt its phase fronts and build its longitudinal field. To see its effect, we also solve the same problem with quadratic=False, which keeps the linear eigenproblem and drops $c$.

The fundamental pair splits again, by about half as much as with the axis in the cross section, because only half of the anisotropy shows up in the cross section: $\varepsilon_{yy} = 2.25$ sits halfway between $\varepsilon_e$ and $\varepsilon_o$. The guide keeps the same fraction of the bulk birefringence as before, $0.0168/0.0199 \approx 0.84$ against $0.0331/0.0400 \approx 0.83$. Dropping the coupling shifts the effective indices by only $5\times10^{-5}$ and $2\times10^{-4}$: $\varepsilon_{yz}$ enters at second order, and here it is small ($0.06$). With a strongly birefringent crystal, the coupling would no longer be a small correction.

To draw the transverse field, we project $E_x$ and $E_y$ onto nodal elements, remove the arbitrary global phase of the eigenvector, and plot $|\mathbf{E}_t|$ with arrows for its direction. The optic axis is marked in green: a double arrow for its projection on the cross section (shorter when the axis is tilted), a dot for its component along the guide.

With the optic axis in the cross section, the mode with the higher index is polarized perpendicular to the axis (ordinary-like, it sees $\varepsilon_o = 2.31$), the one with the lower index along it (extraordinary-like, $\varepsilon_e = 2.19$). For the diagonal axis, the modes are polarized along the diagonals of the square: the crystal, not the shape of the core, sets their polarization. With the tilted axis the same rule holds with the projection of the axis on the cross section, here $\hat{\mathbf{y}}$: the mode polarized along $x$ is ordinary and keeps the higher index, the one polarized along $y$ is extraordinary, with $\varepsilon_{yy}$ as its bulk reference. Both modes are almost linearly polarized, but not exactly: near the corners of the core the field bends, and the small $E_z$ component (not shown) gives the modes a hybrid character that no scalar model captures.

Dispersion curves

Waveguides are characterized by how the effective index changes with frequency. We sweep $k_0 t$ from 3.5 to 15 for the four orientations and record each guided mode together with its polarization: the fraction of the transverse field along the projection of the optic axis on the cross section ($x$ when the axis points along $z$). The sweep takes a few minutes.

Koshiba and co-workers plotted $(\beta/k_0)^2 = n_\text{eff}^2$ against $k_0 t$. Their curves (Figs. 3, 4, 7 and 8 of the paper, digitized for the gyptis example this notebook is based on) appear as grey dots, our modes as coloured circles.

The computed modes fall on the reference curves of all four cases. Every curve starts near the cladding line $(\beta/k_0)^2 = \varepsilon_\text{clad}$, where the mode spreads far into the cladding, and rises towards the core permittivity as the frequency grows and the mode concentrates in the core. Our curves follow the reference almost down to the cladding line; closer still, the mode reaches the PMLs, which give its effective index an imaginary part, and our filter discards it.

The colours show the polarization. With the optic axis along $z$, the fundamental pair is degenerate: its two members (polarized along $x$ and $y$) overlap on the plot. In the three other cases, each family of modes splits into a branch polarized across the projection of the axis (blue, upper) and a branch polarized along it (red, lower), each heading towards its own bulk permittivity: $\varepsilon_o$, and $\varepsilon_e$ (axis in the cross section) or $\varepsilon_{yy}$ (tilted axis).

Modal birefringence

The difference between the effective indices of the two fundamental modes, the modal birefringence, is what sets the length over which a polarization state beats in the guide. We extract it from the sweep: at each frequency, the most confined mode polarized mostly across the axis and the most confined one polarized mostly along it.

The two in-plane orientations give the same curve, as the symmetry argument above predicts, and the tilted axis gives about half of it at every frequency. At low frequency the modes live mostly in the isotropic cladding and the birefringence is small, about a third of the bulk value at $k_0 t \approx 5$. As the frequency grows, the modes concentrate in the core and the birefringence climbs towards the bulk value $n_o - n_e \approx 0.040$, which it reaches only for a core much larger than the wavelength. The extraordinary-like mode lags behind because it also sees a smaller core–cladding contrast, $\varepsilon_e - \varepsilon_\text{clad}$ instead of $\varepsilon_o - \varepsilon_\text{clad}$, so it spreads further into the cladding. The size of the core, relative to the wavelength, is therefore a knob on the birefringence of the guide.

How accurate is the mesh?

We refine the mesh, keeping the ratio between cladding and core element sizes, and follow the two fundamental modes at $k_0 t = 10$ for the optic axis along $x$. As a reference we use quadratic elements on the finest mesh.

The error on the effective index falls roughly as $h^2$, as expected for lowest-order elements on a smooth field. With ten elements across the core, as in the sweep, it is of order $10^{-4}$, well below anything visible on the dispersion curves.

Things to try

Every cell above is editable. A few ideas:

• Stronger anisotropy: lithium niobate at telecom wavelengths has $n_o \approx 2.21$, $n_e \approx 2.14$. Change eps_o, eps_e and the cladding (silica, $\varepsilon \approx 2.09$), and watch the modes become much more confined.
• Any in-plane angle: set the optic axis to [np.cos(a), np.sin(a), 0] and check that the polarization of the fundamental modes follows the angle a.
• Tilt angle: set the optic axis to [0, np.sin(a), np.cos(a)] and follow the birefringence from zero (axis along the guide) to its in-plane value. Compare with quadratic=False: dropping the coupling matters most at a = π/4, since the error grows with $\varepsilon_{yz}^2 \propto \sin^2 a \cos^2 a$.
• Stronger tilt effect: with a strongly birefringent crystal (say $\varepsilon_o = 2$, $\varepsilon_e = 4$), the coupling is no longer a small correction.
• Rectangular core: make the core wider than it is thick (in make_geometry). Shape birefringence then adds to, or cancels, the birefringence of the crystal; find the aspect ratio where the two fundamental modes become degenerate again.
• Near cutoff: enlarge L and the box to follow the curves closer to the cladding line, and see how far the field reaches into the cladding there.