Defining cylindrical space optical resonators through supported mode properties: inverse numerical process

Robert C. Gauthier

doi:10.1364/OE.26.010091

1. Introduction

The optical resonator, configured in various geometries and structured in numerous materials, is one of the key component utilized in integrated and waveguide devices. The usual approach to initial device theoretical design with optical resonators is to select the material properties and geometrical configuration and then utilize numerical solvers to determine the properties of the supported resonator states. One or more of the states are then selected and if required the resonator properties are adjusted to fine tune the states properties and device response. This approach usually requires several iterations to yield a suitable resonator and is a well-established and acceptable design procedure. All current optical design software tools utilize this forward approach as it is linked to manufacturing capabilities and based on parameters sets for available materials.

In this presentation the reverse design procedure is examined. The goal is to propose the optical state and compute the material properties and geometry of the resonator that will support the desired state. The theoretical development presented here is based on the Fourier-Bessel (FB) numerical mode solver for structures that display a high degree of cylindrical symmetry [1]. The presented approach is well suited for the reverse design of numerous axially symmetric structures such as optical and photonic crystal fibers, ring resonators and whispering gallery mode configurations. In Cartesian coordinates the inverse numerical technique would be equivalent to computing the material profile required using the equivalent expressions of the plane wave technique [2]. In spherical coordinates, the Fourier-Bessel-Legendre numerical mode solver has been presented [3] from which the inverse numerical process expressions can be obtained. In a general coordinate system, the curl operators would need to be properly defined as detailed in [4]. Orthogonality of the basis functions utilized in the series expansions is highly desirable but is not an absolute requirement as the volume integrals can be tabulated as done in the FB numerical process [1]. Section 2 indicates the theoretical approach for determining the structure given a particular resonator state. The presentation is general and suitable for all geometries that may be confined in the cylindrical domain.

2. Theoretical development

The theoretical expressions are developed for optical structures that are characterized by relative permittivity and relative permeability tensors and contain variations in the electric and magnetic media properties within the cylindrical domain of the resonator. The orientation of the principle axis for the anisotropy is unrestricted with respect to the resonator. Coordinate axis rotation permits the perpendicular alignment of the resonator plane with the z-axis of a cylindrical coordinate system. Maxwell’s starting curl expressions in a charge and current free environment with time dependence $e^{- j ω t}$ , complex angular frequency $ω = ω_{r} + j ω_{i}$ , are:

\vec{\nabla} \times \vec{E} = - \frac{\partial \vec{B}}{\partial t} = j ω \vec{B}

\vec{\nabla} \times \vec{H} = \frac{\partial \vec{D}}{\partial t} = - j ω \vec{D}

For anisotropic media the constitutive expressions are:

\vec{D} = \overset{\leftrightarrow}{ε} \vec{E} = ε_{o} \overset{\leftrightarrow}{ε_{r}} \vec{E} = ε_{o} [\begin{matrix} ε_{r r} & ε_{r φ} & ε_{r z} \\ ε_{φ r} & ε_{φ φ} & ε_{φ z} \\ ε_{z r} & ε_{z φ} & ε_{z z} \end{matrix}] \vec{E}

\vec{B} = \overset{\leftrightarrow}{μ} \vec{H} = μ_{o} \overset{\leftrightarrow}{μ_{r}} \vec{H} = μ_{o} [\begin{matrix} μ_{r r} & μ_{r φ} & μ_{r z} \\ μ_{φ r} & μ_{φ φ} & μ_{φ z} \\ μ_{z r} & μ_{z φ} & μ_{z z} \end{matrix}] \vec{H}

Using the definition of the free space impedance

Z_{0} = {(μ_{o} / ε_{o})}^{1 / 2}

, and setting

\vec{ℑ} = j Z_{o} \vec{H}

, Faraday’s and Ampere’s law can be written in identical symbolic form as:

\vec{\nabla} \times \vec{E} = \overset{\leftrightarrow}{μ_{r}} \frac{ω}{c} \vec{ℑ}

\vec{\nabla} \times \vec{ℑ} = \overset{\leftrightarrow}{ε_{r}} \frac{ω}{c} \vec{E}

The cylindrical coordinate curl of a vector field, RHS in (3) for a general vector field,

\nabla \times \vec{A}

,

\vec{A} = A_{r} \hat{r} + A_{φ} \hat{φ} + A_{z} \hat{z}

, produces the following three terms:

(\frac{1}{r} \frac{\partial A_{z}}{\partial φ} - \frac{\partial A {}_{φ}}{\partial z}) = A_{\nabla}_{\times} (r) = A_{\nabla}_{\times} (r z) + A_{\nabla}_{\times} (r φ) for \hat{r}

(\frac{\partial A_{r}}{\partial z} - \frac{\partial A {}_{z}}{\partial r}) = A_{\nabla \times} (φ) = A_{\nabla \times} (φ r) + A_{\nabla \times} (φ z) for \hat{φ}

\frac{1}{r} (\frac{\partial [r A_{φ}]}{\partial r} - \frac{\partial A {}_{r}}{\partial φ}) = A_{\nabla \times} (z) = A_{\nabla \times} (z r) + A_{\nabla \times} (z φ) for \hat{z}

In order to perform the required derivatives, each field component for

\vec{A}

can be series expanded using Fourier-Bessel basis functions [1, 4]:

A_{f} = \sum_{q p n} κ_{f} C_{q p n} J_{o} (ρ_{p} \frac{r}{R}) e^{j q φ} e^{j G_{n} z} e^{j k_{z} z}

In the series,

κ_{f}

are the expansion coefficients of field component

f \Rightarrow (r, φ, z)

,

ρ_{p}

are the

p

indexed zeroes of the lowest order Bessel function,

G_{n} = n (2 π / T)

is an

n

indexed reciprocal axial scalar with structure axial repetition length

T

,

q

the azimuthal rotational symmetry index,

C_{q p n}

the basis function normalization constant and

R

is the radial extent of the computation domain. Note that only the lowest order Bessel functions are required and that the radial and azimuthal coordinate directions are decoupled. For a full discussion on the properties and numerical simplifications they provide consult the following [5]. The series form of (7) contains an axial propagation factor,

e^{j k_{z} z}

, with

k_{z}

the axial propagation constant. In the study of optical resonators supporting localized states, the axial propagation constant is set to zero when no out of plane field propagation is involved. For axially propagated fields, such as field propagated in optical fibers,

k_{z} \neq 0

and can be quantified in units of

2 π / T

. Using the expansion series form for each component, (7), Eqs. (4)-(6) are:

A_{\nabla}_{\times} (r) = \sum_{q p n} {(\frac{j q}{r}) κ_{z} + (- j (G_{n} + k_{z})) κ_{φ}} C_{q p n} J_{o} (ρ_{p} \frac{r}{R}) e^{j q φ} e^{j G_{n} z} e^{j k_{z} z}

A_{\nabla \times} (φ) = \sum_{q p n} {(j (G_{n} + k_{z}) J_{o} (ρ_{p} \frac{r}{R}) κ_{r}) + (\frac{ρ_{p}}{R} J_{1} (ρ_{p} \frac{r}{R}) κ_{z})} C_{q p n} e^{j q φ} e^{j G_{n} z} e^{j k_{z} z}

A_{\nabla \times} (z) = \sum_{q p n} {(\frac{κ_{φ}}{r} J_{o} (ρ_{p} \frac{r}{R})) + (\frac{- κ_{φ} ρ_{p}}{R} J_{1} (ρ_{p} \frac{r}{R})) + (\frac{- j κ_{r} q}{r} J_{o} (ρ_{p} \frac{r}{R}))} C_{q p n} e^{j q φ} e^{j G_{n} z} e^{j k_{z} z}

These expressions can be organized into a matrix form and equation set (3) written as:

\nabla \times \vec{A} = {[\begin{matrix} {(4)}_{A} \\ {(5)}_{A} \\ {(6)}_{A} \end{matrix}] = \frac{ω}{c} [\begin{matrix} [B_{r r} + B_{r φ} + B_{r z}] \\ [B_{φ r} + B_{φ φ} + B_{φ z}] \\ [B_{z r} + B_{z φ} + B_{z z}] \end{matrix}] [\begin{matrix} ζ_{r} \\ ζ_{φ} \\ ζ_{z} \end{matrix}]} = \frac{ω}{c} \overset{\leftrightarrow}{B} \vec{ζ}

When (11) represents Faraday’s law (or Ampere’s law), the following change of variables are performed:

[\begin{matrix} \vec{A} \\ \overset{\leftrightarrow}{B} \\ \vec{ξ} \end{matrix}] \Rightarrow [\begin{matrix} \vec{E} \\ \overset{\leftrightarrow}{μ_{r}} \\ \vec{ℑ} \end{matrix}] or ([\begin{matrix} \vec{A} \\ \overset{\leftrightarrow}{B} \\ \vec{ξ} \end{matrix}] \Rightarrow [\begin{matrix} \vec{ℑ} \\ \overset{\leftrightarrow}{ε_{r}} \\ \vec{E} \end{matrix}])

Faraday’s and Ampere’s laws can be combined into a single matrix expression using a 6 component field vector as:

[\begin{matrix} {(4)}_{ℑ} \\ {(5)}_{ℑ} \\ {(6)}_{ℑ} \\ {(4)}_{E} \\ {(5)}_{E} \\ {(6)}_{E} \end{matrix}] = \frac{ω}{c} [\begin{matrix} ε_{r r} & ε_{r φ} & ε_{r z} & 0 & 0 & 0 \\ ε_{φ r} & ε_{φ φ} & ε_{φ z} & 0 & 0 & 0 \\ ε_{z r} & ε_{z φ} & ε_{z z} & 0 & 0 & 0 \\ 0 & 0 & 0 & μ_{r r} & μ_{r φ} & μ_{r z} \\ 0 & 0 & 0 & μ_{φ r} & μ_{φ φ} & μ_{φ z} \\ 0 & 0 & 0 & μ_{z r} & μ_{z φ} & μ_{z z} \end{matrix}] [\begin{matrix} E_{r} \\ E_{φ} \\ E_{z} \\ ℑ_{r} \\ ℑ_{φ} \\ ℑ_{z} \end{matrix}]

In its most general form, Eq. (13) contains 6 equations and 18 unknowns. As a result the most general application of (13) will yield a material profile which is not unique. The non-unique nature of (13) can be displayed by examining the matrix operator form of (13). Equations (3.A) and (3.B) in operator form are:

[\nabla_{\times}] | E = \frac{ω}{c} \overset{\leftrightarrow}{μ_{r}} | ℑ a n d [\nabla_{\times}] | ℑ = \frac{ω}{c} \overset{\leftrightarrow}{ε_{r}} | E

Where

[\nabla_{\times}]

is the curl matrix operator generated from the derivatives in (4)-(6) and

| 〉

represents the three components of the field. When combined to provide a single field solution:

[\nabla_{\times}] | E = {(\frac{ω}{c})}^{2} \overset{\leftrightarrow}{μ_{r}} {[\nabla_{\times}]}^{- 1} \overset{\leftrightarrow}{ε_{r}} | E

The matrix multiplication,

\overset{\leftrightarrow}{μ_{r}} {[\nabla_{\times}]}^{- 1} \overset{\leftrightarrow}{ε_{r}}

, will result in the product of the relative permittivity and relative permeability. The material properties required to support a state can be shared between the electric and magnetic materials provided the product of the two is preserved. Thus a unique solution is not possible. This does provide the designer with the ability to share the material properties required for a state between electric permittivity and magnetic permeability. The numerical examples presented below will employ non-magnetic media and as such all material properties required to support a state are located within the electric permittivity.

The solution to (13) is to yield the material properties which are the elements of the 6 X 6 matrix on the right hand side. Therefore the spacial variations of the 6 field component vector on the RHS is proposed and from these the curls can be computed using the series expansions (7) to generate the 6 component vector on the left hand side. The angular frequency, $ω$ , is chosen for the proposed state, can be complex, enabling the selection of the state’s wavelength and quality factor. Since Faraday’s and Ampere’s laws, when written in curl form are point functions, these laws must be preserved at each point in the computational domain. Therefore, if the field components are specified over a grid discretized cylindrical computation domain, the electric and magnetic material properties can be determined at each of these grid points. The collection of the returned values provides the material space suitable to support the initially proposed mode. Note that once the material space dependence has been determined over the grid space, the material properties can also be expanded in a Fourier-Bessel series and used in further numerical analysis. Determining material and geometry of the resonator using (13) when the state field profile and angular frequency are provided is termed the reverse FB numerical process. Using (13) to compute the states angular frequency and field profile when the material and geometry of the resonator are provided is termed the forwards FB numerical process. The forward and reverse numerical processes, although theoretically related, utilize completely independent numerical engines as they solve for different unknowns. The material space series expansion for each element of the tensors are:

[\begin{matrix} ε_{i j} \\ μ_{i j} \end{matrix}] = \sum_{q p n} [\begin{matrix} κ_{ε_{i j}} \\ κ_{μ_{i j}} \end{matrix}] C_{q p n} J_{o} (ρ_{p} \frac{r}{R}) e^{j q φ} e^{j G_{n} z}

In the most general case for the application of Eq. (13), each row of the matrix system forms a separate equation. The three first row expressions are:

\begin{matrix} \frac{c}{ω} {(4)}_{ℑ} = ε_{r r} E_{r} + ε_{r φ} E_{φ} + ε_{r z} E_{z} \\ \frac{c}{ω} {(5)}_{ℑ} = ε_{φ r} E_{r} + ε_{φ φ} E_{φ} + ε_{φ z} E_{z} \\ \frac{c}{ω} {(6)}_{ℑ} = ε_{z r} E_{r} + ε_{z φ} E_{φ} + ε_{z z} E_{z} \end{matrix}

Any combination of

(ε_{r r}, ε_{r φ}, ε_{r z})

may be chosen to satisfy expression (17-first row) given the field components at the grid point. In general, the selection of the permittivity contributions should be such that adjacent grid points are usually connected through smooth variations in structure parameters, with discontinuities restricted to media interface crossings (for instance core to cladding interfaces). There are several forms of (13) which do not require an ensemble set solution and each relative permittivity element can be obtained individually. The symmetry properties of the sought after resonator configuration may be utilized to provide relationships between the tensor elements and reduce the number of unknowns in the relative permittivity tensor of (2). For isotropic or a material with a diagonal representation in cylindrical space the three equations of (17) decouple in the electric field dependence and may be solved separately. These conditions are presented as numerical examples of the reverse material property determination process. Transposing the theoretical steps presented here to other coordinate system representations is easily performed once a suitable basis expansion set is selected.

3. Computations

Several computation examples of the applicability of the reverse material formulation of the FB numerical mode solver are provided for structures that are uniform infinite in the axial direction. The structures are also selected to be non-magnetic and the relative permeability tensor is diagonal with free space values of 1.000. The lower three rows of (13) or (3.A) provides expressions relating the scaled magnetic field and the curl of the electric field. For a non-magnetic material, or one which is magnetic and uniform over the computation domain, it is sufficient to specify the electric field components of the user defined state. If the electric field profiles are Fourier-Bessel series represented through (7), the expansion coefficients of each field component are determined through orthogonal-basis-function-integration over the computation domain. Equations (4)-(6) can then yield the scaled magnetic field components over the computation domain. Here $(κ_{r}, κ_{φ}, κ_{z})$ are the expansion coefficients for $(E_{r}, E_{φ}, E_{z})$ and $(A_{\nabla}_{\times} (r), A_{\nabla \times} (φ), A_{\nabla \times} (z))$ are related to the scaled magnetic field components after multiplying by, $c / ω$ , with the complex angular frequency selected by the state designer.

3.1 In-plane Bragg rings isotropic resonator in air

An isotropic Bragg ring structure with its axis aligned along the z-axis of the cylindrical coordinate system is uniform-infinite in z and shows no variations with respect to the $φ$ coordinate. The structure displays only a radial dependence in the relative permittivity and the series representation of the tensor elements reduces to a single summation over the lowest order Bessel functions. The relative permittivity tensor is also diagonal with equal elements. The central cylinder, $ε_{r} = 2.400$ and radius of 0.615 µm, plays the role of the core region with the surrounding alternating layers, $ε_{r} = 5.225$ width of 0.169 µm for the high relative permittivity and $ε_{r} = 2.400$ width of 0.2501 µm, forming the cladding region. A surface plot of the structure is shown in Fig. 1 along with a plot of the relative permittivity versus radial direction. The state space can be separated based on polarization, with TM composed of field components $(ℑ_{r}, ℑ_{φ}, E_{z})$ and TE composed of the field components $(E_{r}, E_{φ}, ℑ_{z})$ . The structure is chosen to be non-magnetic and through Eq. (3.A) only one electric field component is required to compute the corresponding magnetic fields for TM states. In general two electric field components are needed for the computation of the single magnetic field for TE states as Eq. (3.B) cannot be used since the relative permittivity profile in the reverse FB application is the unknown. To verify the reverse FB technique, the states of the structure are computed using the forward FB numerical computation technique with $k_{z} = 0.$ The computation radial edge was bordered by a lossy PML layer ensuring zero field values at computation domain perimeter. The proposed structure is discretized in cylindrical space using 2400 radial points (1-D since no $φ$ or z dependence) and the material space series expanded using Eq. (16) with 1000 lowest order Bessel terms $(p_{ε} = {1 \to 1000}, q_{ε} = 0, n_{ε} = 0)$ . The monopole state space was computed using 1000 Bessel terms in the field component series, $(p_{f} = {1 \to 1000}, q_{f} = 0, n_{f} = 0)$ and the state space obtained for modes with wavelengths between 1 and 2 µm are shown plotted in Fig. 2. Two of the states were selected as the input states for application of the reverse FB numerical process. The TM state with wavelength 1.5175 µm ( $ω = 1.2412 \times 10^{15} + j 2.8506 \times 10^{11}$ ) (rad/s) has the axial electric field profile plotted in grey scale. The TE state at wavelength 1.7922 µm ( $ω = 1.0509 \times 10^{15} + j 2.4825 \times 10^{12}$ ) (rad/s) has the azimuthal field component plotted in grey scale. The selected states are highly localized to the central region of the structure.

Fig. 1 Left – Top view of the Bragg Rings resonator structure. Right – Radial plot of the relative permittivity for the Bragg ring structure.

Abstract

1. Introduction

2. Theoretical development

3. Computations

3.1 In-plane Bragg rings isotropic resonator in air

3.2 Uniform infinite, isotropic Bragg rings, axial propagation

3.3 Hexagonal Array

4. Reverse engineering

4.1 Single radial peaked Gaussian

4.2 Double radial identical peaked Gaussian

5. Comments on the technique

6. Summary

Acknowledgment

References and links

Cited By

Figures (16)

Equations (23)

Optics Express