## Abstract

Faraday’s and Ampere’s laws are converted to matrix operator form and rearranged such that the unknown relative permittivity and relative permeability tensors can be determined. The material and geometry of cylindrically symmetric optical resonator structures are determined through the electric and magnetic field component profiles and complex angular frequency of a proposed localized state. This differs from the usual utilization of the electromagnetic wave equations, solving for states given the material properties and geometry. Thus the technique presented here is an inverse numerical process. The theoretical expressions are provided based on a Fourier-Bessel numerical approach which is highly suitable for cylindrical geometry resonators. Without loss of the generality of the technique, examples of resonant structure determination are presented for non-magnetic and diagonal relative permittivity tensor. Axial field propagation is included to demonstrate the design capabilities related to optical fiber and photonic crystal fiber structures.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 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\omega t}$, complex angular frequency $\omega ={\omega}_{r}+j{\omega}_{i}$, are:

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, $\omega $, 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:

**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 $\left({\kappa}_{r},{\kappa}_{\phi},{\kappa}_{z}\right)$ are the expansion coefficients for $\left({E}_{r},{E}_{\phi},{E}_{z}\right)$ and $\left({A}_{\nabla}{}_{\times}\left(r\right),{A}_{\nabla \times}\left(\phi \right),{A}_{\nabla \times}\left(z\right)\right)$ are related to the scaled magnetic field components after multiplying by, $c/\omega $, 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 $\phi $ 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, ${\epsilon}_{r}=2.400$ and radius of 0.615 µm, plays the role of the core region with the surrounding alternating layers, ${\epsilon}_{r}=5.225$ width of 0.169 µm for the high relative permittivity and ${\epsilon}_{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$\left({\Im}_{r},{\Im}_{\phi},{E}_{z}\right)$ and TE composed of the field components $\left({E}_{r},{E}_{\phi},{\Im}_{z}\right)$. 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 $\phi $ or *z* dependence) and the material space series expanded using Eq. (16) with 1000 lowest order Bessel terms $\left({p}_{\epsilon}=\left\{1\to 1000\right\},{q}_{\epsilon}=0,{n}_{\epsilon}=0\right)$. The monopole state space was computed using 1000 Bessel terms in the field component series, $\left({p}_{f}=\left\{1\to 1000\right\},{q}_{f}=0,{n}_{f}=0\right)$ 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 ($\omega =1.2412\times {10}^{15}+j2.8506\times {10}^{11}$) (rad/s) has the axial electric field profile plotted in grey scale. The TE state at wavelength 1.7922 µm ($\omega =1.0509\times {10}^{15}+j2.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.

For the TM polarized state the third row of (13) yields the following expression which must be satisfied at every grid point in the computation domain:

Equation (**19**) and the relative permittivity values obtained are valid for all grid points were the denominator is non-zero. The radial relative permittivity profile for${\epsilon}_{zz}$ computed from the selected state is shown in Fig. 3 as the solid line. The radial dependence of the axial electric field, ${E}_{z}$ dashed line, is plotted as well. The Bragg rings relative permittivity profile is well reproduced using the reverse FB numerical process when the ${E}_{z}$ field profile is utilized as the starting input along with the angular frequency. The relative permittivity profile shows the following features. Slight discontinuities are observed where the field component values pass through zero due to the denominator in (**19**) being zero. The computation accuracy is greater for large amplitude field values as displayed in the central region of the radial domain. In the outer radial region, the field values are lower making (**19**) less accurate and subjected to larger numerical fluctuations. Never the less a reasonable interpolation to a permittivity value of 1.0 can be determined. The region at the extreme outer radial edge displays very large variations as the field is very small and is a manifestation of the presence of PML border imposed on the Bragg rings structure. Note that the numerical process applied here consisted of proposing the structure and then computing the supported fields. Then one of the fields is selected and the reverse FB numerical process is applied. The extremely good match of the relative permittivity profile obtained through the reverse FB application compared to the initial structure is an indication to the highly converged nature of the computation engines utilized.

For monopole TE polarized states the radial field component is zero since$\left({q}_{f}=0,{n}_{f}=0,{k}_{z}=0\right)$. Figure 4 shows the relative permittivity profile and radial ${E}_{\phi}$ field profile plotted to the same radial domain. The numerical process provides the ${\epsilon}_{\phi \phi}$ field profile through the second row equation of (13):

**20**) are at different locations with respect to the zeros of (

**19**) providing addition information on the production of the relative permittivity profile. The TE state is also stronger in the outer radial region providing a more accurate determination of the relative permittivity in that region. The presences of the PML layer still has an effect though of less radial influence on relative permittivity determination.

#### 3.2 Uniform infinite, isotropic Bragg rings, axial propagation

For a uniform infinite isotropic axial waveguide structure the states separate into TE and TM polarized sets. When axial field propagation is present, ${k}_{z}\ne 0$, Eq. (8) re-acquires the azimuthal field dependence and Eq. (9) re-acquires a radial field dependence. Equation (10) does not acquire any ${k}_{z}$dependence other than in the complex exponential factor. For the computation of the relative dielectric profile using TM polarized states, Eq. (19) {third row of (13)} is utilized. Examination of this expressions indicates that the left and right hand side expressions have the common axial propagation factor, ${e}^{j{k}_{z}z}$, which cancel. Thus the relative permittivity profile computed using (19) for ${k}_{z}=0$ is the same for TM states when${k}_{z}\ne 0$. A solution space for state with wavelengths between 1 and 2 µm is shown in Fig. 5 plotted as a function of ${k}_{z}$in units of $2\pi /T$, $T=1$ µm for computation purposes for monopole azimuthal family. TM states are characterized by the horizontal, ${k}_{z}$independent, dotted lines. For TE polarized states and monopole with $q=0$, when ${k}_{z}\ne 0$, Eq. (8) acquires an azimuthal electric field dependence indicating that a radial component to the magnetic field is present for the axially propagated monopole states. When ${k}_{z}\ne 0$ Eq. (9) acquires a radial electric field dependence and an azimuthal magnetic field component is present. Note that the complex propagation constant dependent exponential still cancel for TE polarized states, but due to the derivatives taken with respect to the z axis coordinate, a ${k}_{z}$dependence remains. The dependence of the TE polarized states on the axial propagation constant is displayed in Fig. 5 by the downward curving state dotted lines with increasing ${k}_{z}$. Note in Fig. 5 that the data can be presented in a different format by using the modes effective index, ${n}_{eff}={k}_{z}/k$ versus $k=2\pi /\lambda $. In this representation the effective index of both polarized states are a function of the propagation constant. The determination of the relative permittivity profile using the TE states can in general use either (13) first or second row. However, the first row of (13) cannot be used as (4) = 0 for monopoles. Equation (13) second row utilizes (5) which is ${k}_{z}$ dependent. Figure 6 shows the computed ${\epsilon}_{\phi \phi}$profile for the TE state characterized with ${k}_{z}=0.5\left(2\pi /T\right)$. The Bragg ring structure is re-obtained indicating that the reverse FB numerical process presented here provides a valid technique for the determination of axial waveguide configurations when an axially propagated mode and angular frequency are utilized as inputs selected by the designer.

#### 3.3 Hexagonal Array

A structure that contains material variations in the radial and azimuthal directions, hexagonal array with central hole missing, was utilized as a starting structure. The monopole states for the structure were computed using the forward FB numerical process. The structure on the left of Fig. 7 has air holes in black with a silicon, ${\epsilon}_{r}$ = 12.1104 in white. The holes have a radius of 0.24 µm and lattice constant of 1 µm. The radius of the structure is 2.5 µm. The structure was discretized on a 300 point per micron grid and series expanded using 375 Bessel indices,${p}_{\epsilon}$, and 360 azimuthal indices ${q}_{\epsilon}$. The axial direction is uniform and ${n}_{\epsilon}$ = 0. The ${E}_{z}$field profile for a TM localized state is shown in the central image of Fig. 7. The wavelength of the state is $\lambda $ = 1.625267 µm ($\omega =1.15898\times {10}^{15}$rad/s) which contains only real values as the complex PML bordering region is not present for this structure. The states were computed using 140 Bessel and 150 azimuthal values for the field expansion indices. Given the state the inverse FB process was utilized to determine the ${\epsilon}_{zz}$material profile required to support the state using (**19**). The image on the right of Fig. 7 displays the computed structure shown on a two level gray scale. The hexagonal array is clearly reproduced. As for the Bragg rings structure examined above, the discrepancy in the reproduced profile compared to the original profile occur when the field values pass through zero values. An additional effect is observed due to the lower number of basis functions utilized in determining the states. The state computation is effectively a 2-D computation for the hexagonal array and thus the matrix utilized with the indices selected has order 20580, which is the limit of the memory used in our PC for computations. Greater resolution and better match to the original structure is possible using a larger matrix, but this requires additional solution time. Another way to improve the material space resolution is to compute the structure using two or more states and then selecting peicewise the material profile were field values are not passing through a zero.

## 4. Reverse engineering

Having verified that the inverse Fourier-Bessel technique can reproduce the original material profile based on a selected state, the attention now focuses on providing engineered states and computing the material profile required to support the state. The states are once again restricted to being uniform and infinite in the axial direction. This simplifies the computation process at no loss to the applicability of the technique to address axially variant mode profiles. Two examples are provided. The first consists of a single Gaussian peaked function in the radial direction and displaying an azimuthal rotational symmetry of 10. The second consists of a double identical Gaussian peaked mode profile. Both states have the wavelength of 1.55 µm and are lossless.

#### 4.1 Single radial peaked Gaussian

The modal field ${E}_{z}$ profile is shown plotted in Fig. 8 Left in the mapped $\left(r,\phi \right)\to \left(x,y\right)$ representation and on the right in grey scale polar form. The Gaussian profile is peaked at *r* = 2 µm away from the center and has a σ = 0.25 µm. The azimuthal rotational order is 10 and follows a sinusoidal variation with angle. Through the use of Faraday’s law the radial and azimuthal scaled magnetic field components are computed and shown plotted in Figs. 9 and 10. Equation (**19**) can be used to determine the axial dependence of the relative permittivity tensor element ${\epsilon}_{zz}$. Due to the modal properties selected the relative permittivity shows variations only the radial direction and the radial dependence is shown plotted in Fig. 11 (solid line) along with the radial profile of the proposed state (${E}_{z}$dashed line). Within this plot a zone labelled “Region of Interest” is identified. It corresponds to a section of the radial range for which the proposed modal field is non-zero. Within this region, the relative permittivity plotted is required to reproduce the state. Outside of this region, the proposed modal profile values used in the denominator of (19) should be zero and the corresponding values for the relative permittivity are undefined. However, the modal profile used in the computation process is series expanded using the Fourier Bessel basis functions and due to the finite number of basis functions employed the modal field profile reconstructed outside the Region of Interest is small and non-zero. Thus the field values employed in (19) outside of the Region of Interest are non-zero and finite relative permittivity values are obtained. Increasing the number of basis functions used in the field expansion series has a negative effect on the fluctuations of the relative permittivity values obtained outside the region of interest as the field component is driven closer to zero with increasing basis function number. Knowing that the field value is zero outside the Region of Interest, the relative permittivity values obtained there are not required for state reproduction and can be set to a background value. Note that the region of interest is sufficiently wide to ensure zero field components at the edge of the region. This is similar to selecting the super-cell size in the plane wave technique when defect states in photonic crystals are examined numerically [6]. To verify the suitability of the numerical technique, the radial profile suggested in the Region of Interest of Fig. 11 is placed in a computation domain with relative permittivity of 1.00 elsewhere. The masked relative permittivity profile is shown in Fig. 12. This material profile was used in the forward FB numerical solver and the supported states computed. Figure 13 shows the state obtained at the initial state wavelength of 1.55 µm for the initial proposed ${E}_{z}$field components plotted left side and the recomputed ${E}_{z}$field plotted on the right side. The reproduction of the state indicates that the numerical technique is suitable for the determination of the material profile required to support the state. It also indicates that regions external to the “Region of Interest” can be masked out to a background level provided the proposed state is effectively zero at the region edges. The computation process also indicates that the numerical process is sufficiently converged as the reverse and forward numerical techniques utilize different computation engines and reproduces state profiles accurately.

#### 4.2 Double radial identical peaked Gaussian

A slightly more complicated radial dependence to the ${E}_{z}$ field profile is introduced by adding a second identical Gaussian peaked at *r* = 4 µm as shown in Fig. 14 left in mapped coordinates and on the top right in polar representation. It consist of azimuthal mode order 10 with sinusoidal angular variation. The proposed state is used to compute the radial and azimuthal scaled magnetic field profiles (not shown) through Faraday’s law. These are then utilized in (**19**) to obtain the axial dependence of the relative permittivity as plotted in the solid line of Fig. 15. The dashed line in the figure is a plot of the ${E}_{z}$ field profile. Two non-overlapping Regions of Interest are identified now corresponding to the double peaked Gaussian mode proposed and are the relative permittivity profiles required to support the desired state. Outside these regions the relative permittivity can be masked to the background level as the proposed field values are effectively zero. The two Regions of Interest relative permittivity profiles are similar in shape but the radial outward relative permittivity has larger values as this portion of the relative permittivity must support the same azimuthal order on a greater circumference than the inner modal profile. Figure 16 shows the masked out relative permittivity profile inserted in to the forward FB numerical solver in order to verify the existence of the state. Figure 14 lower right shows a plot of the recomputed state plotted below the proposed state confirming the suitability of the reverse FB technique.

## 5. Comments on the technique

The numerical technique presented here involved only axially uniform and infinite structures making the TE and TM polarized fields separable. In this way the axial relative permittivity could be computed separately from the in plane relative permittivity values. If the proposed states are hybrid states, then the TE and TM field components are no longer separable. If the structure is non-magnetic or magnetic and isotropic, then Faraday’s law can be used to define the scaled magnetic field components. The matrix form of (17) can then be used to obtain the relative permittivity values at each grid point. Note that in this case the three equations are coupled if the anisotropy in cylindrical coordinates is not a diagonal matrix. It is important to note that a relative permittivity tensor for the material present in the structure expressed in $\left(x,y,z\right)$ coordinates and diagonal may not be represented by a diagonal relative permittivity tensor in cylindrical coordinates. Please consult the following when considering anisotropy in general cylindrical resonator configurations [7]. The numerical process will be more involved but the set of equations can be solved for the unknowns provided there are an equal number of equations as unknowns. As is common practice, symmetry properties can be used to reduce the number of unknowns. Computations were performed on a desk-top PC running windows 10 with 4 cores and 32 gigabytes of memory. Typical computations of the relative permittivity profiles require less than one hour to complete.

If the proposed structure is to contain variations in both the relative permittivity and relative permeability, then all 6 field components must be specified as Faraday’s law can no longer be used since the permeability tensor is unknown. The matrix form of (13) is to be solved for the 18 unknowns which requires an equal number of equations. Symmetry properties of the states and imposed symmetry properties of the material distribution can reduce the computation unknowns. In general the computation process is more involved but the overall technique and expressions are those presented here for simpler examples.

As indicated in the single and double Gaussian peaked proposed state section. The relative permittivity profiles obtained in the “Region of Interest” are graded with positive and negative values. Manufacturing such a profile would be difficult and most likely require new material processing techniques or homogenization at the sub wavelength scale. For information on the Maxwell-Garnett effective medium theory please consult [8]. The proposed reverse numerical technique does provide the material profile and geometry required to support the proposed state at the design angular frequency. This information may be used to guide resonator design such that a state close to the desired state may be obtained using available manufacturing techniques.

## 6. Summary

A numerical technique is presented using the Fourier-Bessel spectral approach in cylindrical space for which the material properties of a resonant structure can be determined through the specification of the desired states angular frequency and field profiles. Effectively the process is the inverse of the usual practice of specifying structure and hoping for the desired state. The technique is presented for non-magnetic configurations which results in a unique determination of the relative permittivity in the regions of non-zero state field values. Theoretical expressions provided are general and may be extended for higher demanding field configurations.

## Acknowledgment

The author wishes to acknowledge NSERC for supporting this and other interesting research endeavors.

## References and links

**1. **R. C. Gauthier, “Reformulation of the Fourier-Bessel steady state mode solver,” Opt. Commun. **375**, 63–71 (2016). [CrossRef]

**2. **M. Plihal and A. A. Maradudin, “Photonic band structure of two-dimensional systems: The triangular lattice,” Phys. Rev. B Condens. Matter **44**(16), 8565–8571 (1991). [CrossRef] [PubMed]

**3. **M. A. Alzahrani and R. C. Gauthier, “Spherical space Bessel-Legendre-Fourier localized modes solver for electromagnetic waves,” Opt. Express **23**(20), 25717–25737 (2015). [CrossRef] [PubMed]

**4. **G. Arfken, *Mathematical Methods for Physicists* (Academic Press, 1970)

**5. **R. C. Gauthier and S. H. Jafari, “Dielectric profile segmentation used to reduce computation effort of the Fourier-Bessel mode solving technique,” Opt. Express **23**(11), 14288–14300 (2015). [CrossRef] [PubMed]

**6. **H. Benisty, “Modal analysis of optical guides with two-dimensional photonic band-gap boundaries,” J. Appl. Phys. **79**(10), 7483–7492 (1996). [CrossRef]

**7. **R. C. Gauthier, “Anisotropic resonator analysis using the Fourier-Bessel Mode solver,” Opt. Commun. **410**, 317–327 (2018). [CrossRef]

**8. **T. C. Choy, *Effective Medium Theory: Principles and Applications* (Oxford University Press, 2016).