6.5.3 Comparison with the Waveguide Method

Since the waveguide model is the second frequency domain method close similarities exist to the differential method. In fact we will show that the waveguide model equals the differential method if a zeroth-order vertical discretization scheme is used for the solution of the BVP. Hence the differential method comprises the waveguide model and is thus superior to it.

For that we reconsider (6.32). By comprising the four lateral vectors e_x(z), e_y(z), h_x(z), and h_y(z) carrying the relevant Fourier coefficients (cf. (6.30)) into two new vectors

$$\mathbf{e}(z)= \Big(\mathbf{e}_x(z) \;\; \mathbf{e}_y(z)\Big)^{\m... ...\quad\mathbf{h}(z)= \Big(\mathbf{h}_x(z) \;\; \mathbf{h}_y(z)\Big)^{\mathrm{T}}$$ (6.59)

we can rewrite the ODE system (6.32) to

 

$$\mathbf{h}^\prime(z)= \underline{\mathbf{R}}(z)\mathbf{e}(z)\qqua... ...xt{and}\quad\qquad\mathbf{e}^\prime(z)= \underline{\mathbf{G}}(z)\mathbf{h}(z),$$ (6.60)

whereby the two matrices $\underline{\mathbf{R}}$(z) and $\underline{\mathbf{G}}$(z) are given by

$$\underline{\mathbf{R}}(z)= \begin{pmatrix}\underline{\mathbf{R}}^... ... \underline{\mathbf{G}}^{yx}(z) & \underline{\mathbf{G}}^{yy}(z) \end{pmatrix}.$$ (6.61)

In the waveguide method the medium is divided into thin layers with vertically constant refractive index [8,190,69,73,186,187,188]. The two matrices are assumed to be constant within such a layer, i.e., for z $\in$ [z_l, z_{l + 1}]

$$\underline{\mathbf{R}}(z) = \underline{\mathbf{R}}_l \qquad\quad\text{and}\qquad\quad\underline{\mathbf{G}}(z) = \underline{\mathbf{G}}_l.$$ (6.62)

One vector of unknowns, e.g., h(z), can thus simply be eliminated from the ODE system (6.89). By differentiating the left equation of (6.89) and inserting the right equation of (6.89) we obtain^g

$$\mathbf{e}^{\prime\prime}_l(z)= \underline{\mathbf{G}}_l\mathbf{h... ...e\prime}_l(z)= \underline{\mathbf{G}}_l\underline{\mathbf{R}}_l\mathbf{e}_l(z).$$ (6.63)

The system matrix of the reduced second-order ODE system has a constant coefficient matrix $\underline{\mathbf{G}}_{l}^{}$$\underline{\mathbf{R}}_{l}^{}$. We set it equal to the square $\underline{\mathbf{K}}_{l}^{2}$ of a matrix $\underline{\mathbf{K}}_{l}^{}$, i.e., the N_ODE/2-dimensional second-order ODE system

 

$$\mathbf{e}^{\prime\prime}_l(z)= \underline{\mathbf{K}}^2_l\mathbf... ...quad\underline{\mathbf{K}}^2_l=\underline{\mathbf{G}}_l\underline{\mathbf{R}}_l$$ (6.64)

has to be solved within one layer z $\in$ [z_l, z_{l + 1}]. This ODE system corresponds to the Helmholtz equation (4.2), whereby the wavenumber k is replaced by the matrix $\underline{\mathbf{K}}_{l}^{}$. Note that the system matrix has to be normal to write it as the square $\underline{\mathbf{K}}^{2}_{l}$ of another matrix since only then an eigenvalue decomposition exists, and a matrix $\underline{\mathbf{K}}_{l}^{}$ can be found so that $\underline{\mathbf{K}}^{2}_{l}$ = $\underline{\mathbf{G}}_{l}^{}$$\underline{\mathbf{R}}_{l}^{}$.^h The existence proof of the normal property of $\underline{\mathbf{G}}_{l}^{}$$\underline{\mathbf{R}}_{l}^{}$ as well as the calculation of the square root $\underline{\mathbf{K}}_{l}^{}$ = $\sqrt{\underline{\mathbf{G}}_l\underline{\mathbf{R}}_l}$ can both be performed with the help of (6.22). Since we are only interested in a theoretical comparison the actual numerical evaluation is of minor interest. The formal solution of (6.93) is given by [203, pp. 937-947]

 

$$\mathbf{e}_l(z) = \mathbf{c}_l^+ e^{+\underline{\mathbf{K}}_lz} + \mathbf{c}_l^- e^{-\underline{\mathbf{K}}_lz},$$ (6.65)

whereby the matrix exponential exp( $\pm$ $\underline{\mathbf{K}}_{l}^{}$z) can be calculated with the methods described in [202, pp. 572-577].

The solution (6.94) of the ODE system (6.93) corresponds to the solution of the eigenvalue problem occurring in the waveguide model [187, eq. (33)]. Since the matrix exponential can also be evaluated by an eigenvalue decomposition,ⁱ both methods--waveguide and zeroth-order differential--involve the same numerical operations and thus require the same computational costs. All other steps like matching the solutions of adjacent layers to determine the coefficient vectors c_l^$\pm$ of (6.94) and thereby propagating the solution through the whole simulation interval are equal. This proves theoretically that the differential method comprises the waveguide model as a sub-class, i.e., as the zeroth-order approximation.

As a direction for future work it would be interesting to compare both methods numerically. For a fair and correct benchmark a zeroth-order discretization scheme should be implemented for the differential method. This can be done most efficiently by determining the solution (6.94) of (6.93) since otherwise the higher order discretization is compared with the zeroth-order discretization of the waveguide model. Our implementation is not optimized for that purpose and the benchmark could potentially yield misleading results.

Footnotes

... obtain^g

Also for vertically dependent matrices the vector h(z) can be eliminated from (6.89) yielding $$\mathbf{e}^{\prime\prime}_l(z) - \underline{\mathbf{G}}_l^\... ...hbf{G}}_l\underline{\mathbf{R}}_l\mathbf{e}_l(z) = \mathbf{0}.$$ Future research has to be done to investigate whether a direct solution of the above second-order ODE is advantageous to our implementation which solves the corresponding first-order ODE.

....^h

The eigenvalue decomposition of a normal matrix $\underline{\mathbf{A}}$ is given by $\underline{\mathbf{A}}$ = $\underline{\mathbf{U}}$ $\underline{\mathbf{\Lambda}}$ $\underline{\mathbf{U}}^{\mathrm{H}}_{}$, whereby $\underline{\mathbf{U}}$ is a unitary matrix comprising the eigenvectors $\underline{\mathbf{A}}$u_j = $\lambda_{j}^{}$u_j and $\underline{\mathbf{\Lambda}}$ is a diagonal matrix with the eigenvalues $\lambda_{j}^{}$ as diagonal entries. The square root is defined as $$\underline{\mathbf{A}}^{\frac{1}{2}} = \underline{\mathbf{U... ...}\,\underline{\mathbf{U}}^{\mathrm{H}}=\underline{\mathbf{A}},$$ whereby the unitary property $\underline{\mathbf{U}}^{\mathrm{H}}_{}$$\underline{\mathbf{U}}$ = $\underline{\mathbf{I}}$ of $\underline{\mathbf{U}}$ has been used and the square root $\underline{\mathbf{\Lambda}}^{\frac{1}{2}}_{}$ of the diagonal matrix $\underline{\mathbf{\Lambda}}$ is simply given by taking the square root $\sqrt{\lambda_j}$ of the diagonal entries $\lambda_{j}^{}$. Numerically more efficient and stable methods exist to calculate the square root [202, pp. 572-577]. Note that any matrix polynomial $\underline{\mathbf{A}}^{n}_{}$ can be calculated similarly to the square root, and thus a broad class of functions--including transcendental ones like the exponential--can be calculated from their Taylor series expansions.

... decomposition,ⁱ

See footnote  f on previous page.