6.5.2 Limitations

There are two sources of limitations of the differential method, namely, computational requirements that are getting too large and failure of the stabilized march technique due to strongly damped solution modes resulting in numerical instabilities. The first can be overcome by an increase of the available computer resources, e.g., with the next computer generation it will be possible to solve problems of double size. The second, however, is a principle limitation and thus more serious. We will discuss both impacts from a qualitative point of view.

Computational Requirements. In (6.77) is shown that the computational costs and storage demands grow with the second and the third orders of the number N_ODE of the Fourier coefficients considered, respectively. Since N_ODE is determined by the two truncation frequencies N_x and N_y of the Fourier expansion of the lateral field components (cf. (6.29)), they should be chosen as small as possible. The physical impacts of the truncation cause a low-pass filtering of the field as well as of the geometry and have already been discussed in Section 6.2.5. An analytical investigation of a lateral stepped permittivity including error-bounds is presented in Appendix D. Hence, there is a fundamental tradeoff between computational requirements and accuracy. To be more precise we have to answer the question of how many Fourier modes are required to represent the EM field with a precision that is high enough. For that we reconsider the sampling of the lateral wavevectors introduced in (6.13)

$$k_{x,n}=n\frac{2\pi}{a} \quad\qquad\text{and}\quad\qquad k_{y,m}=m\frac{2\pi}{b}.$$ (6.52)

The following considerations are performed for a homogeneous medium only. However, the results are useful for a qualitative study of a general inhomogeneous medium. Denoting the permittivity of the medium by $\varepsilon$ = $\varepsilon_{r}^{}$ + j$\varepsilon_{i}^{}$ the vertical wavevector component writes to

$$k_{z,nm} = \sqrt{k_0^2\varepsilon- k_{x,n}^2 - k_{y,m}^2}$$

$$= \frac{2\pi}{\lambda} \sqrt{(\varepsilon_r-(n\lambda/a)^2-(m\lambda/b)^2) +j\varepsilon_i},$$ (6.53)

whereby k₀ = 2$\pi$/$\lambda$ is inserted for the wavenumber. Since the vertical dependence E_nm(z) of the EM field is an exponential one, i.e., we have plane waves due to the homogeneity and thus E_nm(z) $\propto$ exp(jk_{z, nm}z), it is described by a damped oscillating function like

$$\mathbf{E}_{nm}(z) \propto e^{j\operatorname{Re}[k_{z,nm}]z} e^{-\operatorname{Im}[k_{z,nm}]z}.$$ (6.54)

The oscillation period $\lambda_{z,nm}^{}$ and the penetration depth $\delta_{z,nm}^{}$ are simply given by

$$\lambda_{z,nm} = \frac{2\pi}{\operatorname{Re}[k_{z,nm}]}\quad\qquad\text{and}\quad\qquad\delta_{z,nm} = \frac{1}{\operatorname{Im}[k_{z,nm}]},$$ (6.55)

whereby the real and imaginary parts of k_{z, nm} are calculated from (6.79) to

 

$$\begin{aligned}\operatorname{Re}[k_{z,nm}] &= \sqrt{\frac{+\varep... ...varepsilon_{r,nm}+\sqrt{\varepsilon_{r,nm}^2+\varepsilon_i^2}}{2}}\end{aligned}$$

with $\varepsilon_{r,nm}^{}$ defined as the real part of (cf. (6.79))

 

$$\varepsilon_{nm}= \underbrace{\varepsilon_r-(n\lambda/a)^2-(m\lambda/b)^2}_{\textstyle\varepsilon_{r,nm}} + j\varepsilon_i.$$ (6.56)

We are now interested in the minimal values of both, the oscillation period $\lambda_{z,nm}^{}$ and the penetration depth $\delta_{z,nm}^{}$, depending on the index pair (n, m). This means that we have to determine (n, m) for which

 

$$\begin{alignedat}{2}\lambda_{z,nm} &\longrightarrow\max \qquad&\o... ...angle}{=}\qquad\operatorname{Im}[k_{z,nm}] &\longrightarrow\max.\end{alignedat}$$

Note that minimal penetration depth means maximal damping. From (6.82) we see that Re[k_{z, nm}] is a monotone increasing function whereas Im[k_{z, nm}] is a monotone decreasing function with respect to $\varepsilon_{r,nm}^{}$, i.e.,

$$\frac{\partial \operatorname{Re}[k_{z,nm}]}{\partial \varepsilon_... ...ad\frac{\partial \operatorname{Im}[k_{z,nm}]}{\partial \varepsilon_{r,nm}} < 0.$$ (6.57)

The extreme values are thus given by (cf. (6.84))

 

$$\begin{alignedat}{4} \lambda_{z,nm} &\longrightarrow\max \quad&\o... ...rt N_y\vert} &= \varepsilon_r-(N_x\lambda/a)^2-(N_y\lambda/b)^2,\end{alignedat}$$

whereby N_x and N_y are the truncation frequencies of the Fourier expansions of the EM field. The mode (0, 0) propagating in vertical direction exhibits the highest oscillation and the least damping, whereas the four most oblique ones ( $\pm$ N_x, $\pm$ N_y) are least oscillating but strongest damped--a result that is physically reasonable. Hence the oscillating nature of the EM field is reflected by the ODE system independently of the truncation frequencies. The truncation frequencies have to be chosen in a way that also damped modes are resembled with sufficient precision. But the higher N_x and/or N_y are, the stronger the damping becomes, which causes stability problems of the algorithm. Before we consider this problem another consequence from (6.86) shall be discussed. As can be seen the ratio between wavelength $\lambda$ and the lateral dimensions a and b of the simulation domain determines $\varepsilon_{r,\vert N_x\vert\vert N_y\vert}^{}$ and thus also the imaginary part Im[k_{z,| Nx|| Ny|}] of the vertical wavevector component (cf. (6.82)). These ratios are the critical quantities for the choice of appropriate truncation frequencies, i.e.,

$$N_x \propto \frac{a}{\lambda} {}\quad\qquad{\text{and}}\quad\qquad{}N_y \propto \frac{b}{\lambda}.$$ (6.58)

The rank of the resulting ODE system and thus the computational requirements are determined by it. For example, in the simulation results presented in Chapter 8 this ratio was approximately four, i.e., for 248 nm lithography the lateral extensions were 1 $\mu$m. This shows that the differential method is only suited for single feature simulation, a situation that applies to all rigorous EM field solvers.

Numerical Instabilities. Numerical instabilities are caused by damped solution modes. The preceding discussion shows that there are two reasons for the existence of such modes, namely, a strongly absorptive medium and high truncation frequencies. The first is self-evident, whereas the second results from the fact that for larger N_x and/or N_y the real part $\varepsilon_{r,nm}^{}$ of the permittivity $\varepsilon_{nm}^{}$ becomes negative (cf. (6.83)). These modes are then called evanescent since in a loss-less medium with $\varepsilon_{i}^{}$ = 0 they propagate in lateral direction and are damped vertically, i.e., in the orthogonal direction to the propagation direction. In a lossy or absorbing medium with $\varepsilon_{i}^{}$ > 0 the real part Re[k_{z, nm}] of the vertical component of the wavevector almost vanishes whereas the negative part Im[k_{z, nm}] becomes dominating. In the limiting case | N_x| $\longrightarrow$ $\infty$ and | N_y| $\longrightarrow$ $\infty$ we have Re[k_{z, nm}] $\longrightarrow$ 0 and/or Im[k_{z, nm}] $\longrightarrow$ $\infty$. This situation is graphically illustrated in Figure 6.6. For a proper function of the stabilized march algorithm the truncation frequencies have to be chosen carefully so that all propagating modes are accounted for, but at most only a few evanescent ones are considered.

The discussion shows that the truncation of the EM field expansion is crucial for the operation of the differential method--not just for the computational requirements but also for the convergence of the algorithm. The same situation applies to the waveguide model, since it it is actually a special case of the differential method as will be shown in the next section.

  

Figure 6.6: The truncation frequencies N_x and N_y in the differenital method have to be chosen so that all propagating modes are reflected by the solution, whereas only a small number of evanescent modes can be considered since they might cause numerical instabilities.