WMMS

Vectorial or semi-vectorial mode analysis of dielectric integrated optical waveguides with two-dimensional light confinement. Multilayer step-index waveguides with rectangular, piecewise constant refractive index profiles are addressed. For a waveguide configuration specified in terms of refractive indices, layer thicknesses, slice widths, and the vacuum wavelength, the online script calculates the propagation constants / effective indices of guided modes and allows to inspect the corresponding optical field patterns. The quasi-analytical solver relies on expansions of the modal fields into locally exact solutions of the mode equations, applied piecewise on the rectangular regions that constitute the waveguide cross section (wave matching method, WMM), as described in Refs. [1, 2, 3, 4]. Further implementational details can also be found in the documentation that accompanies the sources of the WMM solver.

Input

This concerns lossless dielectric optical waveguides with a rectangular, piecewise constant refractive index profile. For a waveguide geometry with N_l interior layers and N_s interior slices, the input mask receives the layer thicknesses t₁, ..., t_Nl, the slice widths w₁, ..., w_Ns, and a matrix of refractive index values n_0,0, ..., n_Nl+1,Ns+1. All dimensions are given in units of micrometers. The figure illustrates an example for a waveguide with N_l = 3 interior layers, and N_s = 1 interior slice:

The entire structure is meant to be constant along the longitudinal z-direction; the light propagates along that axis. Note the orientation of the transverse coordinate axes x and y. In case a waveguide structure is mirror symmetric with respect to a central x-z-plane, this is detected automatically and taken into account for the calculations.

Further, a value for the vacuum wavelength needs to be provided, and a choice of the target polarization, or of the type of analysis, is required. Here the alternatives QTE, QTM, and VEC indicate the semi-vectorial computation of quasi-TE or quasi-TM modes, or the fully vectorial modal analysis, respectively:

• QTE, semivectorial modal analysis. The vertical component E_x of the electric field is assumed to be zero. The horizontal electric component E_y serves as the principal component for the representation of the QTE mode profiles.
• QTM, semivectorial modal analysis. The vertical component H_x of the magnetic field is assumed to be zero. The horizontal magnetic component H_y serves as the principal component for the representation of the QTM mode profiles.
• VEC, fully vectorial simulations. Both transverse directions are treated alike (exception: only vertical symmetry planes are explicitly recognized), there are no assumptions about the polarization of the output modes. Note that simulations of these type require much higher numerical effort than the semivectorial calculations, and that problems due to almost degenerate modes are more likely to be encountered.

Pre-analysis

Given the parameters in the primary input mask, the problem is pre-analyzed with the help of the OMS solver. Limits for potential effective indices of any guided modes supported by the 2-D waveguide are established in the following way.

1. The largest refractive index on one of the four doubly-half-infinite corner rectangles of the cross section decomposition establishes a rough lower limit. A rough upper limit is given by the largest refractive index of one of the interior bounded rectangular regions.
2. The layers stack formed by the regions of the half-infinite left- and rightmost slices that constitute the 2-D cross section are analyzed for guided slab modes, of either polarization. The highest effective index of some slab mode that is guided in these stacks forms an increased lower limit to the effective indices of guided modes of the 2-D waveguide.
3. The layer stacks formed by the interior slices of the 2-D waveguide are subjected to a 1-D slab mode analysis, also for both TE and TM polarization. The highest effective index of any guided slab mode found for those structures, or alternatively the highest refractive index of the half-infinite upper- or lowermost regions (limits to the effective indices of the local continuum of radiation modes), establishes an upper limit for the effective indices of guided modes of the 2-D waveguide.
4. Steps 2 and 3 are repeated for the layer stacks formed by the exterior and interior layers of the 2-D waveguide. This procedure results in the interval for effective mode indices N_eff as shown by the solver.

Next, the symmetry properties of the structure are considered. In case there is a vertical x-z-mirror plane at a central y position, the solver offers analysis options related to modes that are symmetric (SYM) or antisymmetric (ASY) with respect to that plane. Here the interpretation depends on the polarization / the type of analysis.

• QTE: The property denotes even (SYM) or odd (ASY) symmetry of the principal field E_y.
• QTM: The property denotes even (SYM) or odd (ASY) symmetry of the principal component H_y.
• VEC: SYM enforces even symmetry for the mode components H_x, E_y, and H_z, and odd symmetry for E_x, H_y, and E_z; ASY enforces even symmetry for the mode components E_x, H_y, and E_z, and odd symmetry for H_x, E_y, and H_z.

Otherwise, if no mirror symmetry y↔-y is detected, the script continues with mode analysis for fields of no particular symmetry (NOS). Then the number of unknowns is twice as large as for a simulation with prescribed symmetry, requiring a much higher numerical effort.

In a final step, the 2-D waveguide structure is pre-analyzed by the EIMS solver. The script then offers a selection of targets for the further WMM analysis procedure. For the polarization / analysis type pol as selected initially, these options differ with respect to the symmetry alternatives, if applicable, and with respect to the interval of effective indices that is being searched for potential guided modes.

• pol NOS: The full effective index range is investigated with symmetry option NOS.
• pol SYM&ASY: The full effective index range is investigated, in a first pass with symmetry option SYM, then with symmetry option ASY.
• pol SYM: The full effective index range is investigated with symmetry option SYM.
• pol ASY: The full effective index range is investigated with symmetry option ASY.
• EIMS mode id: An interval around the EIMS effective index value for the respective mode is investigated, with the symmetry setting as appropriate for that mode.

The WMM progress window (see the remarks in the next section) shows the approximate modal effective indices as obtained from the EIMS analysis, if any, as thin white dash-dotted lines.

Computational parameters

While details on the algorithm can be found elsewhere, a few remarks are necessary to make reasonable choices on the values of the parameters that affect the mode analysis process.

Separately for each rectangular region with constant refractive index, the mode profiles are expanded into a number of factorizing harmonics or exponentials that solve the basic wave equation. This spectral discretization depends on the trial value for the propagation constant. Each trial function included in the set is characterized by a transverse wave vector. The approximation of the mode profile improves with the width of the interval for these wave vectors and with the density of the selected points in the wavenumber space, i.e. with the number of trial functions. At the same time, the numerical effort increases significantly with the spectral discretization density.

To isolate valid propagation constants, the mode solver searches the interval of prospective effective mode indices for the minima of an error function, the least squares mismatch of the electromagnetic fields on the dielectric boundaries. This error function is shown in the progress window that opens while the actual computations are carried out. The error functions appear as series of sections of more or less narrow parabolas, each of which can be associated with a particular guided mode.

The minimum searching procedure begins with an initial survey of this function. Marching with constant stepsize backwards from the upper end of the specified interval of potential mode wavenumbers, the error function is evaluated for a usually coarse preliminary spectral discretization. Once a minimum in the error function is trapped (i.e. a smaller ordinate appears between two larger ones), the mode solver refines the spectral discretization in several steps and repeatedly fixes the location of the minimum. The process stops when a prescribed final spectral discretization density is reached; this discretization serves to compute the final approximation for the propagation constant. The initial survey then continues until the entire interval is investigated.

The following list of parameters contains quantities both for controlling the spectral discretization and for influencing the behavior of the search algorithm:

• Δα, N_α, α_max: Spectral discretization parameters. Δα is a kind of normalized stepsize for the selection of transverse wavenumbers. The lower this value, the higher is the density of the spectral discretization. Reasonable values are between 0.1 and 0.001. N_α restricts the number of trial functions (for each type of trial functions, rectangle, and basic component …). Reasonable values are between 8 and 30-50. α_max limits the wavenumber interval allowed for combined exponential and harmonic trial functions. A larger value means that more trial functions with steep exponential decays are included. Reasonable values are between 2.0 and 4.0.
  Two sets of these parameters need to be specified, where the first set concerns the initial survey of the error function, aiming at a rough localization of potential guided modes. Proper choices depend on the waveguide structure under investigation. For a simple single mode waveguide, a coarse initial discretization Δα = 0.05, N_α = 10, α_max = 2.0 can be sufficient and efficient. Note that usually a certain spectral density is required for the minima to appear in the error function. Hence, to reveal the minima for all modes of a multimode waveguide, a finer initial spectral discretization Δα = 0.01, N_α = 25, α_max = 2.5 might be required.
  The second set specifies the spectral discretization for finally fixing the propagation constant and for computing the mode profile. For certain benchmark structures, we found that a parameter set of Δα = 0.01, N_α = 30, α_max = 3.0 is sufficient to achieve state of the art accuracy in the propagation constant. These are the default values.
• N_s: Number of equidistant steps to initially scan the effective index interval for minima. Reasonable values range between 10 and 100. Note that the stepsize must be sufficient to bracket all minima that are expected. A value of 10 should be sufficient for an effective index interval with only one guided mode in it, i.e. for a single mode waveguide. Choose higher values for a multimode structure with an error function that comprises a higher number of parabolic sections.
• N_ref, e_exp, S_ref: Parameters influencing the search refinement. N_ref is the number of refinement steps. Values between 1 and 10 are reasonable, 5 is the default. e_exp serves as an exponent for the refinement of the search stepsize. Reasonable values lie between 1.0 and 10.0. A value of 1 means a linear decrease of the stepsize from the initial value, as given by N_s and the effective index interval, to the final value given by the tolerance for the propagation constant. e_exp = 4.0 is the default. S_ref defines the amount of spectral refinement in dependence of the number of the refinement step. Values are to be between 0.0 and 1.0. S_ref = 0.0 means that the three spectral discretization parameters change linearly from the initial to the final values. With S_ref = 1.0 this is a hyperbolic behavior: the spectral discretization remains almost constant for the first steps but is refined steeply in the last few steps. S_ref = 0.5 is the default value.
• δβ: Relative tolerance for the approximation of the propagation constant. Do not choose this value below 1.0e-8 (minimum search for an almost square function with double precision). 1.0e-7 is the default.
• M_shift: Parameter to enforce numerical positivity of the least squares error. The default is 1.0e-10. If the program should stop unexpectedly with a message like Error: Mode analysis failed (the browser console might then shows an error from the original WMM routine like Matrix not positive), first try to alter slightly the initial and final spectral discretization parameters. If that does not help, try setting the M_shift parameter to a moderately (!) higher value (1.0e-9, 1.0e-8).
• p_ext: Weight factor for separating closely spaced modes, here for focusing the solver on actually guided modes, while disregarding close-by leaky-mode-like solutions with higher external field levels. The weight factor is applied to the contributions to the mode normalization on all partly infinite outer rectangles. A value of 1.0 effects uniform weighting; smaller values prefer modes with high field levels in the interior rectangular regions. In case the solver fails to identify a guided mode very close to the lower limit of potential effective indices, you might try to lower the value of this parameter.

Further, in case of fully vectorial calculations, the online solver implements the variant of the WMM algorithm as specified by the parameters vform=HXHY, vnorm=NRMMH, ccomp=CCALL; see the respective remarks in the WMM documentation.

If the solver seems to have missed a guided mode, one should consider the shape of the least squares error function, as displayed in the solver progress window. If no minimum appears at all in the region of the prospective propagation constant, one could try to refine the initial spectral discretization. If a minimum is visible which was not trapped by the mode solver, an increase in the number of initial survey steps N_s and / or a slight alteration of the stepsize for the initial survey Δα may help.

Internally, the solver uses a simple inverse iteration procedure to determine smallest eigenvalues of the matrices that relate to the least-squares error of the field at the dielectric interfaces. That numerical procedure relies on a random initial vector. Hence, repeated runs of the solver for identical problem settings can well lead to different results, what concerns the last digits of the effective index values shown. Obviously, no physical significance should be attached to these.

Output

A table shows, for each guided mode that has been identified by the solver, the propagation constant β (in µm^-1) and the effective index N_eff = β/k, where k = 2π/λ is the vacuum wavenumber associated with the specified vacuum wavelength λ. The mode identifiers indicate the type of analysis (QTE, QTM, VEC) that led to the solution in question, further an indicator (s, a) for the symmetry of the mode (cf. the respective previous remarks), and an index that counts modes with this polarization/analysis type and symmetry in the set of modes as shown, ordered according to effective index (highest first).

Arrow symbols hint at the polarization character of the modes, i.e. the predominant orientation of the electric field, either along the horizontal y-axis (↔), with a major electric profile component E_y, or along the vertical x-axis ( ↕ ), with a major electric profile component E_x. The polarization characters are determined by comparing ratios of integrals of the squared transverse magnetic components of the mode profiles. By construction, QTE modes are predominantly horizontally polarized, while QTM modes are predominantly vertically polarized.

For a structure that supports more than one guided mode, the program calculates the coupling lengths or half beat lengths L_c = π/|β₀ - β₁| that correspond to the interference pattern of each pair of modes with different propagation constants β₀ and β₁.

Integrals of the mode profiles can serve for purposes of normalization, for an assessment of the confinement, or for the evaluation of perturbational expressions, where standard ratios concern the influence of material loss or gain, or phase shifts caused by small changes in refractive index. The script evaluates respective integrals of the squared absolute values |E|², |H|² of the vectorial electric and magnetic fields, and integrals of the longitudinal component S_z of the Poynting vector. Absolute and relative expressions (confinement factors, Γ) are evaluated piecewise for each rectangular part of the cross section with constant refractive index, and for the layers and slices the cross section is composed of.

Field inspection

Mode profile plots show the field (real or imaginary, as discussed above), its absolute value, or the squared field. Thin white lines indicate the dielectric interfaces. After selecting "Plot", the extent of the "color axis" is being adjusted such that it covers the maximum values, determined separately for the electric field strength, magnetic field strength, Poynting vector, and the energy density, over all modes (and all their field components) that have been identified by the solver. This is to make plots comparable. Select the button labeled "↕" to adjust the range of the colormap to the functions that are actually displayed. Note that, due to the approximations involved, some of the profile components of modes of QTE or QTM type are zero by construction.

Click in the figure for a precise evaluation of field levels. A click outside the actual axis closes the coordinate display. The "▯"-button toggles a colorbar. Select a field level on the colorbar to superimpose the field plot with a pseudo-contour at that level. Also here the contour is removed by a click in the colorbar area outside the axes.

Referring to the coordinate system as introduced above, the profiles of the guided modes supported by the present lossless dielectric waveguides are of the form E(x, y) = (E_x, E_y, i E_z)(x,y), and H(x, y) = (H_x, H_y, i H_z)(x,y), where E and H are the electric and magnetic parts of the 2-D mode profile, respectively, depending on the transverse coordinates x and y. The overall phase of the mode can be adjusted such that the component functions E_x, E_y, E_z, H_x, H_y, and H_z are real (note the imaginary units in the expressions for E and H). Consequently, after selection of "E_z" or "H_z", the plot window shows the imaginary part of the longitudinal electric or magnetic component of the mode profile, while for all other components the real part is displayed.

Being solutions of eigenvalue problems, the mode profiles are determined up to some complex constant only. No units are shown for their electric or magnetic fields. Still the given values correspond to a normalization of the modes to unit power flow: The integral over the x-y-plane of the longitudinal component S_z of the Poynting vector evaluates to a power of 1 W. Correspondingly, all electric fields are given in units of V/µm, magnetic fields are measured in A/µm, the components of the Poynting vector S have units of V·A/µm² = W/µm², and the electromagnetic energy density w is measured in W·fs/µm³. In this context the vacuum permittivity and permeability, respectively, are ε₀ = 8.85·10^-3 A·fs/(V·µm) and µ₀ = 1.25·10³ V·fs/(A·µm).

References

[1]  M. Lohmeyer
Wave-matching-method for mode analysis of dielectric waveguides
Optical and Quantum Electronics 29 (9), 907-922 (1997)      

[2]  M. Lohmeyer
Vectorial wave-matching-method mode analysis of integrated optical waveguides
Optical and Quantum Electronics 30 (5/6), 385-396 (1998)      

[3]  M. Lohmeyer, N. Bahlmann, O. Zhuromskyy, P. Hertel
Wave-matching-simulations of integrated optical coupler structures
Proceedings of SPIE, Vol. 3620, Integrated Optics Devices III, 68-78 (1999)      

[4]  M. Lohmeyer
Guided waves in rectangular integrated magnetooptic devices
Cuvillier Verlag, Göttingen, 1999