Dielectric profile segmentation used to reduce computation effort of the Fourier-Bessel mode solving technique

Robert C. Gauthier; Seyed Hamed Jafari

doi:10.1364/OE.23.014288

1. Introduction

The optical properties of localized states attached to a central defect in a photonic crystal and to the central region of a quasi-crystal can be numerically modelled using several well-known techniques [1–3]. Such information is required in the designing stages of numerous optical devices and enhancement of existing structure configuration. If the dielectric structure is considered as finite and presenting rotational symmetry about what would be considered its center, the Fourier-Fourier-Bessel (FFB) numerical technique has been shown to accurately provide frequencies (wavelengths) and field profiles of the centrally localized states [4,5]. The FFB technique cast Maxwell’s equations into an eigenvalue matrix formulation using cylindrical space basis functions, a process similar to casting in linear space using plane waves, PWM [6]. The FFB technique returns information on steady states while PWM primarily returns information on propagated fields. One nice feature of using FFB over PWM is that rotational symmetry present can significantly reduce the order of the matrix and permit the tuning of the matrix for a specific mode type (monopole, dipole, …). Given the available computation resources, the reduced order of the FFB matrix can be increased as higher numerical accuracy is desired on dielectric profiles of larger extent and finer features without exceeding computational resources. Computing the localized states for a circularly symmetric dielectric profile using FFB is a three-step process. The inverse dielectric is discretized and represented by a set of expansion coefficients related to the basis functions; the second step involves the computation of the matrix elements; the third involves solving the matrix to extract the eigenvalues (frequencies and wavelengths) and the eigenvectors (field profiles). Any computation algorithm modifications which speed up the numerical process without sacrificing accuracy are greatly appreciated by users of the technique.

A large number of photonic crystal structures (fibers, resonators, sensors) are configured such that in very specific regions in the structure, the relative dielectric is changed (hole filling for instance) without otherwise changing the underlying dielectric spatial distribution. Some examples include hole filling with an analyte in refractive index sensing [7], gas infiltration [8], humidity in polymers [9], porous materials [10], biological layer deposition [11], bottle resonators [12] and thermal effects [13]. A computation flow algorithm is presented in which the FFB simulation technique is significantly speed enhanced when applied to generating the response curves for these particular configurations. The enhancement results from an elimination of the first two steps in the computation process by recognizing that the inverse dielectric profile can be decomposed into separate profiles for each different dielectric constituent present and that a change in a constituent’s relative dielectric constant is achievable through a rescaling process of the original structure. The next section will present the numerical details of the rescaling process with matrix element expressions provided in Appendix A. The application of the technique is presented on a two dielectric constituent structure composed of Bragg rings and on a three constituent dielectric structure consisting of a hexagonal array of air holes in silicon in which an enlarged central core region can be selectively filled with different refractive index media.

2. Theoretical considerations

The Fourier-Fourier-Bessel (FFB) technique has been presented in several recent publications indicating the details on converting Maxwell’s vector wave equation into an eigenvalue problem [4, 5, 14, 15]. The technique requires that the electromagnetic field components and inverse of the relative dielectric constant profile be expanded using a cylindrical space set of basis functions composed from the lowest order Bessel function for the radial direction and exponentials for the angular and azimuthal directions. The eigen-matrix is constructed by making use of the orthogonality property of the basis functions. The FFB technique requires that the individual matrix elements be computed and the numerical value of each matrix element are strong functions of the expansion coefficients of the inverse dielectric series. The matrix building process, depending on the number of basis functions used in the series, can range from a few minutes, for obtaining indicative results, to several hours for well converged results. Time values are quoted for current high-end desktop PC technology. For dielectric simulation situations in which the spatial distribution profile of the dielectric constituents remains unchanged and only the values of the dielectric constituents change, the matrix building process can be greatly simplified by segmenting the inverse relative dielectric, $Ω = 1 / ε_{r}$ , into an offset value, $Ω_{o f f}$ , and a residue, $Ω_{r e s}$ , with $Ω_{r e s} = Ω - Ω_{o f f}$ . The matrix elements of the modified dielectric constituents can be obtained by rescaling the offset matrix elements separated from the rescaled residue matrix elements. Figure 1, which is the central region of the Bragg resonator of the insert to Fig. 3, demonstrates the procedure. The top left image illustrates a simple two level inverse dielectric profile which contains only radial variations. This structure can be viewed as made up of two profiles, one which is constant over the extent of the dielectric and another which is the deviation between the original dielectric profile and the constant offset level. The two dielectric constituents involved here are silicon, $ε_{r, S i} = 12.1104$ and oxide, $ε_{r, o} = 2.4000$ . The lower left image presents the same dielectric constituent distribution but with two different inverse dielectric values. This structure can also be considered as an offset, $Ω'_{o f f}$ , and residue, $Ω'_{r e s}$ . Since the two left hand side profiles have the same dielectric spatial distribution, the offset and residue dielectric profiles are related by:

Ω'_{o f f} (r, ϕ, z) = (\frac{ε_{o f f}}{ε'_{o f f}} = K_{o f f}) Ω_{o f f} (r, ϕ, z)

Ω'_{r e s} (r, ϕ, z) = (\frac{\frac{1}{ε'_{1}} - \frac{1}{ε'_{o f f}}}{\frac{1}{ε_{1}} - \frac{1}{ε_{o f f}}} = K_{r e s}) Ω_{r e s} (r, ϕ, z)

with

Ω_{1} = 1 / ε_{1}

,

Ω'_{1} = 1 / ε'_{1}

,

Ω_{o f f} = 1 / ε_{o f f} = 1 / ε (r = R)

, and

Ω'_{o f f} = 1 / ε'_{o f f} = 1 / ε' (r = R)

. Within the confines of the FFB technique, the inverse of the dielectric is usually expanded into a series of the form:

Ω = \sum_{p n m} κ^{Ω} J_{o} (ρ_{p} \frac{r}{R}) e^{j m ϕ} e^{j G_{n} z}

where

κ^{Ω}

are the expansion coefficients determined through an orthogonality integration procedure. Also,

J_{o} ()

is the lowest order Bessel function with roots indicated by

ρ_{p}

and distinguished using the integer index

p \geq 1

. The exponentials are used to Fourier decompose both the angular and azimuthal coordinate dependences. The angular index,

m

, is an integer and

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

are azimuthal related reciprocal lattice vectors (one direction thus scalar) with

T

the axial height of the dielectric extent and

n

an integer. Segmenting the dielectric into a constant offset and residue pushes all of the dielectric variations of the profile into the residue and the inverse dielectric series can be written as:

Ω = Ω_{o f f} + \sum_{p n m} κ^{Ω_{r e s}} J_{o} (ρ_{p} \frac{r}{R}) e^{j m ϕ} e^{j G_{n} z}

Notice that the series expansion of the original dielectric profile expressed in (2) is forced to zero at

r = R

through the roots of the Bessel function. Selecting the offset value in (3) to coincide with the inverse dielectric value at

r = R

removes the series imposed zero of Ω at

r = R

[16]. In addition, dielectric series representation convergence is achieved with fewer radial basis function terms using (3). Maxwell’s wave equations in charge free, current free region, expressed using the inverse dielectric profile, Ω, and time variation,

e^{- j ω t}

are:

Ω \vec{\nabla} \times \vec{\nabla} \times \vec{E} = {(\frac{ω}{c})}^{2} \vec{E}

\vec{\nabla} \times Ω \vec{\nabla} \times \vec{H} = {(\frac{ω}{c})}^{2} \vec{H}

With the introduction of the dielectric profile segmentation and the possibility of rescaling each of the dielectric profiles as shown in Fig. 1, equation set (4) becomes:

K_{o f f} [Ω_{o f f} (\vec{\nabla} \times \vec{\nabla} \times \vec{E})] + K_{r e s} [Ω_{r e s} (\vec{\nabla} \times \vec{\nabla} \times \vec{E})] = {(\frac{ω}{c})}^{2} \vec{E}

and

K_{o f f} [Ω_{o f f} (\vec{\nabla} \times \vec{\nabla} \times \vec{H})] + K_{r e s} [(\vec{\nabla} \times Ω_{r e s} \vec{\nabla} \times \vec{H})] = {(\frac{ω}{c})}^{2} \vec{H}

The eigen-matrix is developed by expressing the fields as a spectral decomposition using the same basis functions as in (2-3):

(\vec{E}, \vec{H}) = \sum_{p n m} κ^{(\vec{E}, \vec{H})} J_{o} (ρ_{p} \frac{r}{R}) e^{j m ϕ} e^{j (G_{n}) z}

The FFB technique returns steady states localized at the center of the dielectric profile. As a result, the desired states display negligible field amplitude at the radial domain edge

r = R

. The fact that the field series expansion in (7) is forced to zero through the roots of the Bessel functions is not an issue when centrally localized states are required. A similar reasoning is introduced in PWM when defect states are required and a super-cell approach is utilized [17]. Expression (8) summarizes the form of the eigen-matrix for a two constituent dielectric configuration.

Fig. 1 Top – Inverse relative dielectric constant profile for the central region of the Bragg structure (insert to Fig. 3). It can be considered as the summation of a constant offset, $Ω_{o f f}$ , and a residue $Ω_{r e s} = Ω - Ω_{o f f}$ as shown on the top-right. Lower left - Inverse relative dielectric constant profile in which the spatial distribution is the same as the top left trace. The lower-right indicates that the offset and residue values can be obtained by scaling the values of the top-right image using (1).

Download Full Size | PDF

K_{o f f} [\begin{matrix} R_{r, o f f} & Φ_{r, o f f} & Z_{r, o f f} \\ R_{ϕ, o f f} & Φ_{ϕ, o f f} & Z_{ϕ, o f f} \\ R_{z, o f f} & Φ_{z, o f f} & Z_{z, o f f} \end{matrix}] [\begin{matrix} R_{E, H} \\ Φ_{E, H} \\ Z_{E, H} \end{matrix}] + K_{r e s} [\begin{matrix} R_{r, r e s} & Φ_{r, r e s} & Z_{r, r e s} \\ R_{ϕ, r e s} & Φ_{ϕ, r e s} & Z_{ϕ, r e s} \\ R_{z, r e s} & Φ_{z, r e s} & Z_{z, r e s} \end{matrix}] [\begin{matrix} R_{E, H} \\ Φ_{E, H} \\ Z_{E, H} \end{matrix}] = {(\frac{ω}{c})}^{2} [\begin{matrix} R_{E, H} \\ Φ_{E, H} \\ Z_{E, H} \end{matrix}]

The 3 by 3 matrices along with their column vector are produced by the corresponding expression in Eqs. (5) and (6). If the spatial dependence of the dielectric constituents is unchanged, these matrices can be determined for an initial structure and then rescaled for generating the matrix for other dielectric values of the same dielectric layout. Each of the 3 by 3 matrices in (8) represents a matrix block, with order equal to the number of basis functions used, formed during the orthogonality integration step. The upper case letter represents the wave equation vector component and the lower case subscript letter indicates the applied orthogonality field component. The additional factor to the subscript indicates if the matrix block is formed from the offset contribution or residue contribution to the total wave equation. The matrix block generating expressions are provided in Appendix A. The possibility of making use of the offset and residue can also be applied to dielectric systems examined using PWM. For PWM the offset level can be selected as the smallest inverse dielectric constant value and the residue formed with only positive values.

Should the dielectric be composed of more than two constituents and that spatial dielectric distribution is unchanged, scaling through the offset and residue can still be applied. The offset is still taken as the value at the radial boundary. A residue spatial profile is produced for each of the dielectric values present as shown in Fig. 2. The scaling factors are calculated using the expressions given in (1) applicable to each of the residue profiles. The resulting eigenmatrix is composed through the following expression which now contains a summation over the number of separate residue profiles:

Fig. 2 Left –Bragg structure consisting of three different relative dielectric constant materials. Right - The offset level is indicated by the dashed line. The central core region forms one residue spatial profile with zero values in the confining ring region. The second residue profile is formed from the wide segments of the confining rings. The second residue profile has zero values in the central core region. The eigenmatrix can be rescaled to accommodate structures having the same spatial relative dielectric distribution.

Download Full Size | PDF

K_{o f f} [\begin{matrix} R_{r, o f f} & Φ_{r, o f f} & Z_{r, o f f} \\ R_{ϕ, o f f} & Φ_{ϕ, o f f} & Z_{ϕ, o f f} \\ R_{z, o f f} & Φ_{z, o f f} & Z_{z, o f f} \end{matrix}] [\begin{matrix} R_{E, H} \\ Φ_{E, H} \\ Z_{E, H} \end{matrix}] + \sum_{i = 1}^{n} {K^{i}_{r e s} [\begin{matrix} R^{i}_{r, r e s} & Φ^{i}_{r, r e s} & Z^{i}_{r, r e s} \\ R^{i}_{ϕ, r e s} & Φ^{i}_{ϕ, r e s} & Z^{i}_{ϕ, r e s} \\ R^{i}_{z, r e s} & Φ^{i}_{z, r e s} & Z^{i}_{z, r e s} \end{matrix}]} [\begin{matrix} R_{E, H} \\ Φ_{E, H} \\ Z_{E, H} \end{matrix}] = {(\frac{ω}{c})}^{2} [\begin{matrix} R_{E, H} \\ Φ_{E, H} \\ Z_{E, H} \end{matrix}]

The matrix generating expressions are still provided by those indicated in Appendix A. It is interesting to indicate that the technique of rescaling the dielectric can also be applied when thermal effects takes place which result in a change of refractive index through the thermo-optic effect and isotropic thermal expansion [13]. The rescaling applies to the change in the relative dielectric constants and can be used in Eqs. (8) or (9). The thermal expansion (or contraction) can be accounted for by invoking the scalability of Maxwell’s equations [18]. Thus the eigenvalues returned from the diagonal matrix would need to be rescaled by the ratio of the expansion coefficient and temperature change.

4. Numerical computation

The rescaling technique is presented for two well-known circularly symmetric dielectric configurations, the Bragg rings and hexagonal photonic crystal with central defect region. Families of curves are produces for the wavelength variation in the localized steady states versus constituent relative dielectric constant variation.

4.1. Bragg Rings

The Bragg structure is shown is Fig. 3-insert and has the following parameters used to define the reference offset and residue values for the structure. The background, edge and white portions of the rings (ring width of 0.111 µm) have a relative dielectric constant of silicon, 12.1104, giving $Ω_{o f f} = 0.0826$ . The black regions in the insert correspond to the low relative dielectric constant constituent (2.400) and define the central region, radius 0.615 µm, and the 0.250 µm wide black sections of the Bragg rings. The reference residue value computed from the dielectric constants is $Ω_{r e s} = Ω - Ω_{o f f} = 0.3341$ . The starting set of Bragg ring parameterswere chosen such that they are a quarter wavelength in width for a free-space wavelength of 1.55 µm. Figure 3 shows a plot of the residue rescale factor computed for the low relative dielectric varied from 1.000 to 12.1104. This factor is used in matrix expression (8) to rescale the elements of the residue matrix. At a constituent value of 2.400 the rescale factor is 1 as expected for the original structure parameters. For the computations presented the background value is unchanged thus the offset rescale factor is 1.00. Figure 4 shows a segment of the eigen-state wavelengths returned versus constituent dielectric constant rescaled based on the values of Fig. 3. All eigenvalues are computed using the magnetic field eigen-matrix and the eigenvalues are plotted for eigenvectors dominated by the $H_{z}$ field component. The constituent index was varied in 0.01 increments resulting in 1111 computations of the eigenvalues returned from (8) for a matrix tuned to solve for the monopoles only. The Bessel index, p, ranged from 1 to 200 and ensures sufficiently converged results. The structure has no variations in φ and z, thus the indices q and n are set to 0. The double edge arrow placed in the dominant bandgap region is positioned when the dielectric constituent has the design value of 2.400. The centrally localized monopole state shows up as a bandgap line and is indicated by the single arrow. The open regions at lower wavelengths correspond to the second and third bandgap regions for the structure. This figure also displays the closing of the bandgaps as the dielectric contrast is reduced. Figure 5 is a compilation of the centrally localized states returned when monopoles, dipoles, hexapoles and quadrupoles tuned matrices are solved for the low relative dielectric range of 1 to 12.1 in 0.01 increments. Selected mode field profiles are plotted for each type of mode confirming the centrally localized nature of the states. Eigenvalues are plotted for eigenvectors dominated by the $H_{z}$ field component. During the development of the enhanced solver, results were compared directly to results obtained using the original solver. Numerical results in the comparisons are effectively identical. Also note that the enhanced solver programs are separate from the original solver programs. This ensure that the integrity of original mode solvers which have been tested and confirmed accurate with other techniques such as FDTD and PWM.

Fig. 3 Plot of the residue scaling factor versus relative dielectric constant $ε'_{1}$ .Residue rescale factor = 1 corresponds to original dielectric structure for which the matrix elements in (8) are determined. For this computation example $ε'_{o f f} = ε_{o f f}$ . Insert shows the Bragg structure with white being silicon and black the oxide.

Download Full Size | PDF

Fig. 4 Plot of the eigen-wavelengths in the 0.5 to 5 µm range versus, upper scale - matrixfile index, lower scale altered relative dielectric constant, for the 1111 eigen-matrices solved for monopole states with eigenvectors dominated by the $H_{z}$ field component. Double arrow used to locate the primary bandgap and single arrow points to a state localized to the central region of the dielectric structure. Two other bandgaps at lower wavelengths (higher frequencies) are observed in the figure. One of the bandgaps displays that a higher order localized monopole state is present.

Download Full Size | PDF

Fig. 5 Monopole, dipole, quadrupole and hexapole states of the primary bandgap dominated by the $H_{z}$ field component in the eigenvector plotted versus file index and altered relative dielectric constant. Selected mode profile plotted displaying the rotational symmetry and localized nature of the bandgap states.

Download Full Size | PDF

4.2 Hexagonal array

The hexagonal array structure is shown in the insert of Fig. 6. The background material in white is silicon and serves to define the offset inverse dielectric value (same as for the Bragg ring structure). The low dielectric value regions shown in black are initially set to air. The central air core region has a radius of 0.8 µm. The hexagonal array of holes has a lattice constant of 0.630 µm and hole radius of 0.270 µm. The band structure of the hexagonal array was computed using the plane wave method (PWM) and shows a bandgap between wavelengths 1.360 and 2.290 µm when the holes are air filled. The pinching off of the bandgap as the hole-background dielectric contrast is reduced is plotted in the main trace of Fig. 6, also computed using PWM. The lower axis corresponds to the hole region relative dielectric constant.

Fig. 6 Primary band gap range fo the silicon background hexagonal array plotted versus relative dielectric constnat within the holes of the array. Insert - Hexagonal array dielectric profile with central large defect (white – Silicon,black - air) shown in insert. Additional line in bandgap indicates the location of the monopole states (see Fig. 7 for enlarged state plot).

Download Full Size | PDF

Full infiltration is examined first where all low relative dielectric constant values are rescaled (center and holes of the hexagonal array). This corresponds to a single residue dielectric profile and matrix expression (8) is applicable when rescaling the dielectric values. This is followed by a partial infiltration where the relative dielectric constant is changed only in the central region. The holes in the hexagonal array remain at a value of 1.0 and the central core region’s value is changed making this a two residue dielectric profile configuration governed by matrix Eq. (9).

4.3 Full infiltration of the array

Full infiltration consists of changing the dielectric constant in the central region and in the holes of the hexagonal array. The effect of the hexagonal array’s bandgap pinching off for decreasing dielectric contrast results in a structure that losses the monopole state with increasing infiltration dielectric value. For each infiltration value the eigen-matrix was obtained through the rescaling of Eq. (8) with 70 Bessel terms $(p = 1 \to 70)$ and 61 angular $q = - 180 \to 180$ in steps of 6 due to 6-fold rotational symmetry of the structure. No azimuthal terms are required in the basis functions as the z direction is uniform. The order of the square matrix produced is 12810 which require approximately 3 hours to create the first time and a few minutes to create when rescaled (excluding file load and save times which depend greatly on the file format employed). The monopoles state wavelengths dominated by the $H_{z}$ field component are plotted in Fig. 7 as a function of the infiltration relative dielectric constant. This plot is produced with data points indicating greater than 75% field confinement to the central disk region. The missing data point at $ε_{r} = 1.5$ results from the monopole state mixing with an edge state and reduces the field confinement to below the threshold of 75%. Selected mode profiles are plotted showing a well confined mode when air filled and radially extended mode profiles close to cut-off. The behavior of higher order modes (dipole, quadrupole, …) were not examined for the hexagonal array structure.

Fig. 7 Monopole states in the 1.2 and 2.2 µm wavelength range with 75% or greater field confinement to the central disk region of the hexagonal structure. States with eigenvector dominated by the $H_{z}$ field component plotted. Lower axis represents the file index and altered relative dielectric constant varied from 1.0 to 2.4 in 0.1 increments (15 data files). Due to the pinching off nature of the bandgap as the dielectric contrast decreases the localized monopole is not supported for infiltration relative dielectric constant above 2.3.

Download Full Size | PDF

4.4 Core infiltration

The situation were only the central core inner region is selectively filled “increased” results in the hole regions of the hexagonal array remaining as air and maintains the hexagonal array’s bandgap intact. The plot in Fig. 8-left shows the $H_{z}$ field component dominated states obtained for all eigen-wavelengths between 1.2 and 2.2 µm which correspond to the monopole states which overlaps the band gap region of the hexagonal array. The plot was obtained by increasing the dielectric value from 1.0 to 12.1 in 0.1 increments. This gives 111 eigenvalue data files, each obtained through the rescaling in (9) and computed using the same basis function space as used for full infiltration calculations. The centrally localized states can be isolated from the full eigenvalue spectrum by reconstructing the field profiles and retaining those eigenvalues that have a field with 75% confinement to the central disk. The plot of Fig. 8-right shows the result of applying the field profile based filter. The localized states produce a sequence of mode lines with the number of states increasing with center “core” dielectric value. Selected monopole mode profiles with a circle around the eigen-wavelength are plotted in this figure. The outer circle represents the radial extent of the dielectric profile while the inner circle passes through the second inner most set of holes in the hexagonal array (see grey circles of insert of Fig. 6).

Fig. 8 Left – All monopole states in the 1.2 and 2.2 µm wavelength range returned from the eigen-solver dominated by the $H_{z}$ field component. States with 75% or greater field confinement to the central region are plotted on the right. The field confinement filter identifies and extracts localized states from all states provided by the solver. States identified with an addition O are plotted in top down column format for the intensity of the $H_{z}$ field component. They display the various mode symmetries composing the monopole family.

Download Full Size | PDF

The following discussion is included in order to clarify the classification of the states plotted in Fig. 8 as monopole states. Figure 9-top shows a collection of the dielectric series expansion coefficients plotted as a function of the azimuthal index for integers between −33 and + 33. The hexagonal dielectric structure of Fig. 6 has rotational symmetry of 6 thus non-zero expansion coefficients exist for $q_{Ω}$ indices which are an integer multiple of the rotational symmetry of the dielectric (6). The base rotational symmetry of a monopole state is zero as for a pure monopole there is no variation in the field profile with in-plane angle. Therefore through the matrix element generating expressions the non-zero elements in the matrices correspond to locations were the field expansion coefficients are related to basis functions with $q_{f}^{*} - q_{f} = q_{Ω}$ . Thus the $q_{f}$ possible non-zero expansion coefficients for the monopole follow the sequence $q_{f} = {..., - 30, - 24, - 18, - 12, - 6, 0, 6, 12, 18, 24, 30, ...}$ . This sequence determines the family grouping for the monopole states. States in the family with expansion coefficient spectrum dominated by $q_{f} = 0$ define the usual monopoles with high field amplitude at the center of the pattern and may show field variation with angle due to the presence of non-zero field expansion coefficients at higher $q_{f}$ values. Such a field profile has its intensity plotted in Fig. 9-left along with its series expansion coefficients. The expansion coefficients are non-zero for integer multiple of the dielectric symmetry and dominated by the $q_{f} = 0$ values. States in the family dominated by expansion coefficients with $q_{f} = \pm 6$ show a strong 6-fold angular symmetry. Such a field profile has its intensity plotted in Fig. 9-right along with its series expansion coefficients. The expansion coefficients are non-zero for integer multiples of the dielectric symmetry and dominated by the $q_{f} = \pm 6$ values. States in the family dominated by expansion coefficients with $q_{f} = \pm 12$ will show a strong 12-fold angular symmetry. The FFB technique returns all states for the mode family not just the lowest and most highly confined states. For a dielectric with rotational symmetry of 6 as presented here, mode families exist for monopoles (0), dipoles (1), quadrupoles (2) and hexapoles (3). The largest mode family order is given by half of the rotational symmetry of the dielectric when even and half of the (rotational symmetry - 1) for odd rotationally symmetry dielectric profiles. Mode families with (1), (2), (3) were not examined here.

Fig. 9 Top – Real part of the expansion coefficients for the 6-fold rotationally symmetric hexagonal array. Only the coefficients in the angular index range from −33 to + 33 are shown. The 6-fold rotational symmetry present forces all expansion coefficients with angular index which are not an integer multiple of the rotational symmetry to zero. Left – Positive side and absolute value of the expansion coefficients of a strong monopole (center region relative dielectric constant of value 11, wavelength 1.461 µm). Right – Positive side and absolute value expansion space for the hexapole state at wavelength 1.323 µm. Both these states show the 6-fold symmetry of the dielectric. The center figure is dominated by the zero order angular index and defines the usual monopole while the left state is dominated by the angular order of 6 of a hexapole. Inserts show field amplitude profiles.

Download Full Size | PDF

5. Conclusion

We have shown that the splitting of a dielectric profile into a constant offset and residue can significantly speed up the matrix building process of the Fourier-Fourier-Bessel mode solver. The splitting technique can be applied when the spatial distribution of the dielectric constituents remains the same but the dielectric values are changed. The new matrices can be obtained through a simple rescaling multiplication factor. Such a technique may find application in numerical analysis of index of refraction sensors, gas infiltration, humidity in polymers, porous materials, biological layer deposition, bottle resonators and thermal effects. In FFB the rotational symmetry present and the possibility of tuning the matrix to a particular mode type is first applied to reduce the order of the matrix, combined with the speed enhancements of the rescaling property provides the researcher with quick access to numerical data.

Appendix A:

The matrix elements within each of the 3 by 3 blocks in Eqs. (8) and (9) can be computed using the expressions and the integrals listed below. A fair degree of bookkeeping is required in placing the matrix elements into the array as the matrix element indices must match the sequence of the indices in the field’s column matrix. The term in $‖ ‖$ indicate the condition for non-zero matrix elements relative to the azimuthal index q, ( $q + q^{Ω} = q^{*}$ ) and axial index n, ( $n + n^{Ω} = n^{*}$ ). The volume of the computation domain is $V = π R^{2} T$ . Equations are written for the residue matrix corresponding to magnetic field wave equation. The offset matrix element generating expressions are obtained by setting dielectric derivative related terms to zero. This is easily accomplished by setting $(p_{Ω}, q_{Ω}, G_{n_{Ω}}) = (0, 0, 0)$ and converting the integral expressions to integrals involving only the field contributing Bessel. The integrals convert from the product of three Bessel functions to the product of two Bessel functions. All integrals can be tabulated as they are all normalized to unit radial extent and used in all applications of FFB technique. The FFB technique can also be used to compute modes for axially propagated fields. The axial propagation constant, $k_{z}$ is added to $G_{n}$ in (7) and all occurrences of $G_{n}$ related to the field components are changed to $(G_{n_{r, φ, z}} + k_{z})$ .

R_{r, r e s} = \frac{1}{V {[J_{1} (ρ_{p_{r}^{*}})]}^{2}} \sum_{Ω} κ^{Ω} ([q_{r} q_{Ω} + q_{r}^{2}]) T_{- 1} (0_{p_{r}}, 0_{p_{Ω}}, 0_{p}_{r^{*}}) + [G_{n_{r}} G_{n_{Ω}} + G_{n_{r}}^{2} R^{2} T_{1} (0_{p_{r}}, 0_{p_{Ω}}, 0_{p}_{r^{*}})] {‖ 4 π^{2} ‖}_{q n}

φ_{r, r e s} = \frac{1}{V {[J_{1} (ρ_{p_{r}^{*}})]}^{2}} {\sum_{Ω} κ^{Ω} j (q_{φ} + q_{Ω}) [T_{- 1} (0_{p φ}, 0_{p_{Ω}}, 0_{p}_{r^{*}}) + [ρ_{p φ} T_{o} (0_{p φ}, 0_{p_{Ω}}, 0_{p}_{r^{*}})]] ‖ 4 π^{2} ‖}_{q n}

Z_{r, r e s} = \frac{1}{V {[J_{1} (ρ_{p_{r}^{*}})]}^{2}} \sum_{Ω} κ^{Ω} (- j R ρ_{p_{z}}) (G_{n_{Ω}} + G_{n_{z}}) T_{1} (1_{p_{ϕ}}, 0_{p_{Ω}}, 0_{p}_{r^{*}}) {‖ 4 π^{2} ‖}_{q n}

R_{φ, r e s} = \frac{1}{V {[J_{1} (ρ_{p_{φ}^{*}})]}^{2}} {\sum_{Ω} κ^{Ω} (- j) ([q_{r} T_{- 1} (0_{p_{r}}, 0_{p_{Ω}}, 0_{p}_{φ^{*}}) + q_{r} ρ_{p_{r}} T_{1} (1_{p_{r}}, 0_{p_{Ω}}, 0_{p}_{φ^{*}}) - q_{r} ρ_{p φ} T_{o} (0_{p_{r}}, 1_{p_{Ω}}, 0_{p}_{φ^{*}})]) ‖ 4 π^{2} ‖}_{q n}

ϕ_{ϕ, r e s} = \frac{1}{V {[J_{1} (ρ_{p_{ϕ}^{*}})]}^{2}} \sum_{Ω} κ^{Ω} R^{2} [(G_{n_{ϕ}} G_{n_{Ω}} + G_{n_{φ}}^{2} + \frac{ρ_{p_{ϕ}}^{2}}{R^{2}}) T_{1} (0_{p_{φ}}, 0_{p_{Ω}}, 0_{p}_{φ^{*}}) + \frac{1}{R^{2}} T_{- 1} (0_{p_{φ}}, 0_{p_{Ω}}, 0_{p}_{φ^{*}}) + \frac{ρ_{p_{ϕ}}}{R^{2}} T_{o} (0_{p_{φ}}, 1_{p_{Ω}}, 0_{p}_{φ^{*}}) - \frac{ρ_{p_{φ}} ρ_{p_{Ω}}}{R^{2}} T_{o} (1_{p_{φ}}, 1_{p_{Ω}}, 0_{p}_{φ^{*}})] {‖ 4 π^{2} ‖}_{q n}

Z_{φ, r e s} = \frac{1}{V {[J_{1} (ρ_{p_{φ}^{*}})]}^{2}} {\sum_{Ω} κ^{Ω} R q_{z} (- G_{n_{Ω}} - G_{n_{z}}) T_{o} (0_{p_{z}}, 0_{p_{Ω}}, 0_{p}_{φ^{*}}) ‖ 4 π^{2} ‖}_{q n}

R_{z, r e s} = \frac{1}{V {[J_{1} (ρ_{p_{z}^{*}})]}^{2}} \sum_{Ω} κ^{Ω} (j R) [G_{n_{r}} T_{0} (0_{p_{r}}, 0_{p_{Ω}}, 0_{p}_{z^{*}}) - G_{n_{r}} ρ_{p_{r}} T_{1} (1_{p_{r}}, 0_{p_{Ω}}, 0_{p}_{z^{*}}) - G_{n_{r}} ρ_{p_{Ω}} T_{1} (0_{p_{r}}, 1_{p_{Ω}}, 0_{p}_{z^{*}})] {‖ 4 π^{2} ‖}_{q n}

φ_{z, r e s} = \frac{1}{V {[J_{1} (ρ_{p_{z}^{*}})]}^{2}} \sum_{Ω} κ^{Ω} (- R G_{n_{ϕ}}) (q_{φ} + q_{Ω}) [T_{0} (0_{p_{ϕ}}, 0_{p_{Ω}}, 0_{p}_{z^{*}})] {‖ 4 π^{2} ‖}_{q n}

Z_{z, r e s} = \frac{1}{V {[J_{1} (ρ_{p_{z}^{*}})]}^{2}} {\sum_{Ω} κ^{Ω} [ρ_{p_{z}}^{2} T_{1} (0_{p_{z}}, 0_{p_{Ω}}, 0_{p}_{z^{*}}) + (q_{z} q_{Ω} + q_{z}^{2}) T_{- 1} (0_{p_{z}}, 0_{p_{Ω}}, 0_{p}_{z^{*}}) - ρ_{p_{z}} ρ_{p_{z^{*}}} T_{1} (1_{p_{z}}, 1_{p_{Ω}}, 0_{p}_{z^{*}})] ‖ 4 π^{2} ‖}_{q n}

T_{i} (A, B, C) = \int_{0}^{1} J_{A} (ρ_{A} τ) J_{B} (ρ_{B} τ) J_{C} (ρ_{C} τ) τ^{i} d τ

Acknowledgement

The authors wish to thank NSERC for supporting this line of photonics research.

References and links

1. R. C. Gauthier and K. Mnaymneh, “FDTD analysis of 12-fold photonic quasi-crystal central pattern localized states,” Opt. Commun. 264(1), 78–88 (2006). [CrossRef]

2. Bahrampour, A. Seifalinezhad, A. Iadicicco, and A. R. Bahrampour, “Zero birefringence 8-fold photonic quasicrystal (QC) fiber,” in Proceedings of IEEE Conference on Photonics Third Mediterranean (IEEE, 2014), pp. 1–3.

3. S. Guo and S. Albin, “Simple plane wave implementation for photonic crystal calculations,” Opt. Express 11(2), 167–175 (2003). [CrossRef] [PubMed]

4. R. C. Gauthier and M. Alzahrani, “Cylindrical space Fourier-Bessel mode solver for Maxwell’s wave equation,” Adv. Mater. 2(3), 32–35 (2013).

5. S. R. Newman and R. C. Gauthier, “Fourier-Bessel analysis of localized states and photonic bandgaps in 12-fold photonic quasi-crystals,” J. Opt. Soc. Am. A 29(11), 2344–2349 (2012). [CrossRef] [PubMed]

6. J. C. Slater, Quantum Theory of Molecules and Solids McGraw-Hill, 1963, Chap. 7.

7. G. E. Town, W. Yuan, R. McCosker, and O. Bang, “Microstructured optical fiber refractive index sensor,” Opt. Lett. 35(6), 856–858 (2010). [CrossRef] [PubMed]

8. I. Shavrin, S. Novotny, A. Shevchenko, and H. Ludvigsen, “Gas refractometry using a hollow-core photonic bandgap fiber in a Mach-Zehnder-type interferometer,” Appl. Phys. Lett. 100(5), 051106 (2012). [CrossRef]

9. W. Yuan, L. Khan, D. J. Webb, K. Kalli, H. K. Rasmussen, A. Stefani, and O. Bang, “Humidity insensitive TOPAS polymer fiber Bragg grating sensor,” Opt. Express 19(20), 19731–19739 (2011). [CrossRef] [PubMed]

10. H. Bao, K. Nielsen, H. K. Rasmussen, P. U. Jepsen, and O. Bang, “Fabrication and characterization of porous-core honeycomb bandgap THz fibers,” Opt. Express 20(28), 29507–29517 (2012). [PubMed]

11. G. Emiliyanov, P. E. Høiby, L. H. Pedersen, and O. Bang, “Selective serial multi-antibody biosensing with TOPAS microstructured polymer optical fibers,” Sensors (Basel) 13(3), 3242–3251 (2013). [CrossRef] [PubMed]

12. G. Senthil Murugan, M. N. Petrovich, Y. Jung, J. S. Wilkinson, and M. N. Zervas, “Hollow-bottle optical microresonators,” Opt. Express 19(21), 20773–20784 (2011). [CrossRef] [PubMed]

13. M. Luo, Y. G. Liu, Z. Wang, T. Han, J. Guo, and W. Huang, “Microfluidic assistant beat-frequency interferometer based on a single-hole-infiltrated dual-mode microstructured optical fiber,” Opt. Express 22(21), 25224–25232 (2014). [PubMed]

14. S. R. Newman and R. C. Gauthier, “Fourier-Bessel expansions of localized light in disorder dielectric lattices,” Proc. SPIE 8269, 82690T (2012). [CrossRef]

15. R. C. Gauthier, M. A. Alzahrani, and S. H. Jafari, “Theoretical examination of the slot channel waveguide configured in a cylindrically symmetric dielectric ring profile,” Opt. Commun. 329, 154–162 (2014). [CrossRef]

16. F. Bowman, Introduction to Bessel Functions Dover Publications, New York, 1958, pp. 17–18.

17. K. Busch, C. T. Chan, and C. M. Soukoulis, “Techniques for band structures and defect states in photonic crystals,” Photon. Band Gap Mater. 315, 465–485 (1996). [CrossRef]

18. M. Diem, T. Koschny, and C. M. Soukoulis, “Wide-angle perfect absorber/thermal emitter in the terahertz regime,” Phys. Rev. B 79(3), 033101 (2009). [CrossRef]

Dielectric profile segmentation used to reduce computation effort of the Fourier-Bessel mode solving technique

Abstract

1. Introduction

2. Theoretical considerations

4. Numerical computation

4.1. Bragg Rings

4.2 Hexagonal array

4.3 Full infiltration of the array

4.4 Core infiltration

5. Conclusion

Appendix A:

Acknowledgement

References and links

Cited By

Figures (9)

Equations (21)

Optics Express