QuEPS

The scalar polarized 2-D Helmholtz equations are addressed, on a rectangular computational domain with transparent boundary conditions that permit guided wave in- and outflux. Given a problem specification in terms of interface positions, a matrix of refractive index values, polarization and wavelength parameters, and guided wave input, the script determines modal output amplitudes (elements of the scattering matrices), and the power levels associated with guided and nonguided directional outgoing waves (transmittances / reflectances, power balance). Facilities for detailed inspection of the optical electromagnetic field are provided, including animations of the harmonic oscillations in time, with options for exporting figures and data.

Input

Cartesian coordinates x and z span the 2-D plane of interest. For convenience, we refer to these directions as vertical (x) and horizontal (z). The terms "top", "right", "bottom", and "left" are used to identify the boundaries of the computational domain. For a rectangular circuit with N_l inner layers and N_s inner slices, the input mask receives the vacuum wavelength λ, a specification of the polarization, a matrix of refractive index values n_l,s, values t_l for the thicknesses of the inner layers, and values w_s for the widths of the inner slices. All dimensions are meant in micrometers. At least one inner layer and one inner slice is required. The figure illustrates the relevant geometry, where blue elements relate to the boundaries of the computational domain:

The optical electromagnetic field is assumed to be constant along the y-direction. One then distinguishes waves that are polarized perpendicular to the x-z-plane (TE polarization), and waves that are polarized parallel to that plane (TM polarization).

In certain cases the default-fill-functionality of the script can ease the data entry. Select the appropriate number of inner layers and slices, and Clear the input mask. Enter refractive indices and width for one slice. Then Fill will repeat the values of refractive indices and width row-wise over all slices. Likewise, if only a single refractive index value, width, or thickness is provided, the Fill procedure enters these values into all other fields.

Guided wave excitation

Initialize starts a pre-analysis of the structure. Guided modes supported by the slab waveguides that cross the boundaries of the computational domain are identified. Initial amplitudes can be specified for each of these modes.

• Input- and output amplitudes relate to power-normalized modes, and are specified for propagation coordinates at the positions of the respective boundaries of the computational domain (cf. the slab mode solver OMS for a detailed analysis of the access modes).
• The guided input and output modes are mutually power orthogonal, and power orthogonal to all other, non-guided incoming and outgoing fields. The guided incoming and outgoing optical power can thus be calculated as the sum of the squared modal amplitudes. The solver shows the total incoming guided power P_in associated with the given amplitudes.
• In case only a single modal amplitude is set to 1, the complex amplitudes associated with the outgoing modes represent directly the elements of the scattering matrix of the circuit, where the external boundaries act as port "planes". The squared output amplitudes can then be interpreted as guided wave reflectance and transmittances.

Computational parameters

The solver works on a cross-shaped computational domain, with an inner rectangular computational window (x, z) ∈ [x_b, x_t] × [z_l, z_r] (cf. the figure above). That window is specified by the distances Δx_t, Δz_r, Δx_b, Δz_l between the outermost interface positions, and the boundary locations.

The electromagnetic field is expanded into the sets of normal modes (1-D slab modes with Dirichlet boundary conditions) supported by the refractive index profiles associated with the layers and slices that make up the circuit. Inside the inner computational window, quadridirectional expansions into modes that propagate in the ±x- and ±z-directions apply. Bidirectional expansions in the half-infinite exterior regions ensure that the boundaries of the inner window become transparent for outgoing (scattered) waves from the interior. The numbers M_x, M_z of terms taken into account for these expansions are specified through the number of terms spend for each (vacuum) wavelength-unit in the respective x-, z-extension of the inner window.

Remarks:

• The computational effort required for the simulations depends on the number of expansion terms, and, for given spectral density, on the size of the computational domain (where, due to the semi-analytical method, the precise functional dependence differs from configuration to configuration).
• While the four boundaries of the inner computational window are transparent for outgoing waves, the four corner points act as point scatterers (and also the boundaries of the "arms" of the cross shaped domain that begin at these corners act as perfect reflectors). In case of configurations with strong radiative losses, the corners should thus be positioned (in a follow-up simulation, by readjusting / increasing Δx_t, Δz_r, Δx_b, Δz_l) away from any pronounced scattered waves.
• The boundary distances should be chosen such that the evanescent tails of the incoming and outgoing guided modes are sufficiently covered by the respective port extension.
• While the default spectral density of 10/λ can be sufficient for a configuration with low or moderate refractive index contrast, higher values are typically required for structures with stronger contrast. An insufficient spectral density typically shows up by a violation of the power balance (the total output power, summed over all four boundaries, deviates from P_in), accompanied by nonphysical discontinuities of the electromagnetic fields.
• In singular instances it may happen that the solver produces nonphysical results, even with actually (i.e. for a slightly differing configuration) appropriate parameters. Presumably this is caused by an accidental degeneracy of modes, of different slices/layers, in the internal expansions. You might wish to try a simulation with slightly modified computational window in that case.
• Please refer to the original proposal [1] of the QUEP method for details of the algorithm.

Convergence with respect to all computational parameters, up to a reasonable degree, should be ensured before attempting any interpretation of results. You might wish to try the structural data and the values for computational parameters summarized in the square resonator example to gain some inital experience with the solver and the user interface.

Output

Upon completion of a simulation, the solver lists the complex amplitudes of guided modes of different order (if any), that cross the boundary lines in the four directions, and the levels of outgoing guided and total optical power, per boundary line and for the entire domain. Depending on the specific excitation setting, the values can be directly interpreted as elements of the scattering matrix of the circuit, or as reflectance and transmittance levels, respectively.

The solver then offers a panel with controls for detailed inspection of the optical field solution. Referring to the coordinate system as introduced above, this concerns electric fields E and magnetic fields H that depend on the spatial coordinates x, z and on time t, with angular frequency ω, as

All fields are constant along the y-direction.

The frequency-domain electric and magnetic fields E, H depend on the in plane coordinates x and z only. For TE polarization, these fields are of the form E = (0, E_y, 0), H = (H_x, 0, H_z). Likewise, TM-polarized fields are of the form E = (E_x, 0, E_z), H = (0, H_y, 0). The plot functionality refers to the complex components of these field functions. Further, the x- and z-components S_x, S_z of the time-averaged Poynting vector and the time averaged energy density w can be inspected. Note that the y-component of the Poynting vector vanishes for the present y-constant fields. Single components ψ ∈ {E_x, . . . , H_z} of the complex-valued frequency domain electromagnetic field relate to time-varying physical fields Ψ(x, z, t) = Re ψ(x, z) exp(i ω t). The global phase of the solution can be adjusted; the animations show the respective component at recurring points in time, equally distributed over the period of 2π/ω.

Being solutions of linear differential equations, the field profiles scale linearly with the incident waves. No units are shown for the electric or magnetic field components. Still the given values correspond to a normalization to unit incident power per unit length (y-direction), of the incoming modes (W/µm). 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).

Plot controls are offered that permit to enlarge (+) and to reduce (-) the size of the figure, to Export the curve data, to export the figure in SVG format, and to Detach the figure into a separate browser window. The boundary positions of the mapped rectangle can be adjusted; defaults are averages between the outermost interfaces and the computational window boundaries. The plot rectangle can extend beyond the inner computational window, e.g. in order to visualize corner effects. 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.

Example

The following list of parameters reproduces the field solution that corresponds to the resonant state of a square 2-D resonator with perpendicular bus waveguides, as reported in section 3.4 of Ref. [1] (different: position of the coordinate origin; size of the computational window, precise spectral density, precise plot region). Quantities are as introduced in the preceding paragraphs of this page.

• TE,  λ = 1.55µm,
• N_l = 3,  N_s = 3,

• n4,0 = 1.0,	n4,1 = 1.0,	n4,2 = 1.0,	n4,3 = 3.4,	n4,4 = 1.0,
  n3,0 = 1.0,	n3,1 = 3.4,	n3,2 = 1.0,	n3,3 = 3.4,	n3,4 = 1.0,
  n2,0 = 1.0,	n2,1 = 1.0,	n2,2 = 1.0,	n2,3 = 1.0,	n2,4 = 1.0,
  n1,0 = 3.4,	n1,1 = 3.4,	n1,2 = 3.4,	n1,3 = 3.4,	n1,4 = 3.4,
  n0,0 = 1.0,	n0,1 = 1.0,	n0,2 = 1.0,	n0,3 = 1.0,	n0,4 = 1.0,

• t₃ = 1.786µm,  t₂ = 0.355µm,  t₁ = 0.1µm,
• w₁ = 1.786µm,  w₂ = 0.385µm,  w₃ = 0.1µm,
• input (TE₀):   a_left = 1.0,  a_top = a_right = 0.0,
• Δx_t = Δz_r = Δx_b = Δz_l = 4.0µm,
• # spectral terms x: 20/λ,  # spectral terms z: 20/λ.

The solver predicts a guided reflectance of 22% (output to the left port), and transmittances of 22% (right port) and 46% (top port).

Reference

[1]  M. Hammer
Quadridirectional eigenmode expansion scheme for 2-D modeling of wave propagation in integrated optics
Optics Communications 235 (4-6), 285-303 (2004)