1. Introduction

Manufacturing imperfections and tight tolerances in photonic crystal (PC) structures present great challenge on the road of transferring this promising technology into the domain of commercial applications. Much work has been done to study an impact of imperfections on the performance of PCs [1–17]. It was established quite generally that small randomness in PC geometry and/or material constants leads to an overall reduction in a band gap size, as well as an increased back scattering and radiation loss, while stronger randomness can lead to appearance of localized impurity states. In the majority of theoretical studies various simplified models of randomness are assumed. Such models are typically chosen for simplicity of parameterization of a particular type of randomness, or because some modelling methods could only handle certain types of geometries. In 1D PC multilayers [1, 2] one typically considers disorder in the thickness and value of a dielectric constant of individual layers. In 2D planar PCs and microstructured fibers [3–13] one frequently considers random displacement of features from an underlying ideal lattice, disorder in a feature size (radius of a hole, for example), disorder in the refractive index, distortion of feature shapes (ellipticity), as well as roughness of walls [14, 15] which is sometimes modelled by protrusions of some average characteristic hight and width. In 3D PCs derived from lithographical techniques and opals [16, 17] additional imperfections are stacking faults, and surface roughness. In all these calculations disorder parameters are scanned from small to large and conclusions are drawn about their relative impacts. Propagation parameters can be sensitive functions of disorder parameters. For example, power in the back scattered modes from wall roughness in a planar TIR waveguide scales quadratically with roughness hight (perturbation theory) and is a very sensitive function of a roughness correlation length. Thus, for a rigorous comparison of theoretical estimates with experimental observations one has to be precise about the types and statistical importance of realistic imperfections.

The goal of this paper is to understand which are the statistical parameters of importance when describing disorder in PC lattices, and then to characterize such parameters quantitatively by analyzing high resolution experimental images of 2D planar slab PCs. We find that three intuitive sets of parameters are necessary to create a comprehensive statistical model of PC imperfections. First set of parameters describe coarse properties of features such as radius, ellipticity and other low angular momenta components in a feature shape. Such coarse variations of a shape can be either deliberately designed or result from an imperfect manufacturing process. Another set of parameters describe roughness of feature edges on a nanometer scale, that is wall roughness, which is ultimately determined by the random physical processes of electron scattering in a resist, resist development and etching. A final set of parameters describes deviations of feature centers from ideal periodic lattice.

When interpreting Scanning Electron Microscopy (SEM) images one has to always keep in mind that SEM image is a convolution of an imperfect fabrication (see Appendix) and an imperfect SEM capture. Reconstruction of actual dielectric profiles from SEM images can be a very non-trivial task far beyond the scope of this paper. In the following we apply our statistical analysis to SEM images assuming that they represent true dielectric profiles. Developed statistical formulation is, however, general and can be applied to any images.

The paper in organized as follows. We first characterize coarse variations and wall roughness in individual features from which the periodic lattice is constructed. Next, we characterize an imperfect lattice formed by individual features. Finally, we evaluate errors in deduced statistical parameters due to finite image resolution. We demonstrate our approach by analyzing various high resolution images of 2D planar PCs manufactured by e-beam lithography and detailed in the following publications: InP/InGaAsP/InP [22, 23], Air/InP membranes [24, 25], Air/Si membranes [26, 27], and SiO ₂/Nb ₂ O ₅/SiO ₂ [28].

2. Statistical description of feature edges

In this section we introduce a parameterization model to characterize feature shapes. It is a general property of statistical fitting for the error to decrease with the number of parameters included in the model. The challenge is to define as small as possible set of parameters and yet to capture physically significant variations in a feature shape. Here we present a “common-sense” criterion to decide on a minimal number of parameters required to decompose the shape into a regular curve plus edge roughness noise. In what follows we adapt statistical methods developed for characterization of rheology of complex surfaces [18, 19, 20] to describe fabrication imperfections in planar PCs.

2.1. Fitting a single feature edge

First, object recognition algorithm [21] (see discussion in section 4) is used to extract circular features and their edges Fig. 1(a). We start with images of highest resolution having few features an example of which is Fig. 1, [22] with resolution of 0.46nm. Edge of each feature is then fitted to extract coordinates of its center, radius, ellipticity and higher order Fourier components in the edge shape Fig. 1(b). Particularly, we define an edge objective function

where r_edge (θ_i ) is a distance from a feature center (X ₀,Y ₀) to an edge point (X_i,Y_i ),

Feature radius in (2) is an m=0 term, implicit m=1 terms correspond to the feature center coordinates (X ₀,Y ₀), while m=2 expansion coefficients in (2) define feature ellipticity as ${\delta }_{el}=\sqrt{{\left(A_2^M\right)}^2+{\left(B_2^M\right)}^2}$ with an angle of a major axis defined by cos(2θ_el )=$B_{\mathit{2}}^M$/δ_el . For a given M there are (1+2M) fit parameters. We fit these parameters by minimizing an objective function (1) (finding zeros of its derivatives) using standard multidimensional Newton method, with a root mean square (RMS) of a fit error defined as $\sigma \left(M\right)=\sqrt{Q_{\mathit{edge}}^M}$. In what follows M≥1, where for M=1 only the feature radius and coordinates of its center are fitted.

When more Fourier components are used in a fit, fit error σ(M) monotonically decreases, and the change in the individual fit coefficients becomes smaller than sub-nanometer image resolution for even modest values of M<10. Thus, for example, for a single hole in Fig. 1

Note that, in general, assumption that feature edge r _edge (θ_i ) can be fitted with a single valued analytical curve $r_{\mathit{\text{fit}}}^M$ (θ) might not be true on a small enough scale (in a particular case of Fig. 1(c) this scale is below 2nm) where rough feature edge is fractal-like. For all the analyzed pictures we find that analytical form (2) is applicable on a larger than a nanometer scale.

If there are n=[1,N_f ], N_f >1 features in the image each containing $N_{\mathit{\text{edge}}}^n$ edge points, their shapes are first fitted individually. Then, all the relevant parameters and correlation functions are averaged over the features. For example, if $R_0^n$, ${\sigma }_n^2$ (M), [C,Γ,S${]}_n^M$ (λ) are the radius, variance and correlation functions of a feature n, then their averaged counterparts are defined as

2.2. Coarse parameters defining feature edges

The goal of this section is to define coarse parameters that characterize feature edges “globally” such as radius, ellipticity, quadruple contribution, etc. and to establish their relative relevance. Particularly, we consider statistics of deviations of feature edges from the corresponding smooth fits when only a small number of low angular components are included in a fit. We define random variable describing the fit error as

where n=[1,N_f ], i=[1,$N_{\mathit{\text{edge}}}^n$ ]. Note that <${\delta }_r^M$ >=0 as it is proportional to the derivative of (1) with respect to a feature radius. Variance of ${\delta }_r^M$ is given by (4). In what follows, we find that for large enough values of M distribution of the remaining edge roughness ${\delta }_r^M$ can always be fitted by a Gaussian probability density distribution

 

Fig. 1. (a) Image of a hole together with a detected edge. (b) Shape of a rugged edge is fitted with Fourier series in θ. Smooth curve is an M=1 circle fit. (c) On a scale <2nm hole edge can not be represented by a single valued analytical function $r_{f\mathit{\text{it}}}^M$ (θ). (d) Edge roughness is self-similar on very different scales suggesting fractal description.

Download Full Size | PDF

In Figs. 2(a,b) we analyze an image with resolution 0.86nm of a PC lattice [22] where direct e-beam writing was used and circular features were coded as polygons with 12 sides. In Fig. 2(a) probability density distribution (PDD) of ${\delta }_r^M$ is presented as a function of the number of angular momenta components M in a fit. We observe that image data (solid curves) and Gaussian distribution (dotted curves) with mean 0 and variance σ²(M) defined by (4) match very well for all M≥1, indicating that error of a fit is random and normally distributed. As the number of angular components M in a fit increases RMS of ${\delta }_r^M$ becomes smaller Fig.2(b). For M>1 we observe a slow decrease of fit error RMS with the number of angular momenta components, suggesting that there is no simple coarse description of a feature shape (such as ellipticity), and that higher order angular components contribute substantially. Simplest statistical model of a feature edge for this image can be defined in terms of an average radius of a circle R=124.68±1.76nm and a RMS deviation of an edge from such a circle σ(1)=3.09nm. One can also specify an averaged over features ellipticity and any number of other higher order angular components in an effort to provide a more complete statistical model of a feature edge. Thus, for example, if ellipticity is included in a fit, statistical model of coarse parameters defining feature edges will be specified by an average radius of an underlying circle R=124.68±1.76nm, feature ellipticity δ_el =2.64±1.17nm, direction of an ellipse major axis θ_el =18±36 ^o , and a RMS deviation of an edge from an elliptical fit σ(2)=2.32nm. Note that averaged over features ellipticity has a large deviation from its mean, as well as an ill-defined direction of an ellipse major axis, thus signifying that in this system ellipticity is not a clearly identifiable property of the features, but rather a part of an edge roughness. Frequently, when ellipticity and other low order angular components in a feature shape are not deliberate (as in Fig.2(a,b)), their contributions are more natural to account for in terms of statistical properties of an edge roughness (see next section), rather than through individual coarse parameters.

In Figs. 2(c,d) we analyze an image with resolution 2.45nm of a PC lattice [26] where a pattern of elliptical features of graded radii was written to form a central vertical waveguide. In Fig. 2(c) we present probability density distribution (PDD) of ${\delta }_r^M$ as a function of the number of angular momenta components M included in a fit.We observe that image data (solid curves) and Gaussian distribution (dotted curves) with mean 0 and variance σ²(M) defined by (4) match well for M≥2 indicating that after including ellipticity, the error of a fit is mostly random and normally distributed. From Fig. 2(d) it is clear that for this structure ellipticity and to a lesser extent quadruple contributions are important when describing feature shapes. That is, from Fig. 2(d) we observe that inclusion of ellipticity M=2 reduces fit error by a factor of 4, inclusion of M=3 does not change fit error substantially, while inclusion of a quadruple M=4 reduces the error further by almost a factor of 2. As the number of angular components in a fit increases beyond M=4 we observe a slow decrease (power law as established later) of RMS of fit error. Statistical model of coarse parameters defining feature edges can, thus, be specified by an average radius of an underlying circle R=178.46±7.82nm, feature ellipticity δ_el =22.11±1.87nm, direction of an ellipse major axis θ_el =87.3±1.5o, and a RMS deviation of an edge from an elliptical fit σ(2)=3.42nm. Note that averaged over features ellipticity has a very small deviation from its mean, and a well-defined direction of an ellipse major axis, thus signifying that in this system ellipticity is an intrinsic property of the features.

 

Fig. 2. (a,c)Probability density distribution of fit error for different number of angular momenta components M in a fit. (b) RMS of fit error decreases slowly as the number of angular momenta components M in a fit increases, suggesting that there is no simple coarse description of a feature shape. (c) RMS of fit error decreases dramatically when ellipticity M=2 of a feature is included in a fit, suggesting ellipticity as a dominant coarse parameter.

Download Full Size | PDF

In general, relevance of coarse parameters can be judged from dependence of a fit error on the number of included angular components M. Typically, we observe that coarse parameters corresponding to the first several angular components M≤10 are of major importance, and their inclusion leads to a considerable reduction in the fit error (Fig. 2(d), 1≤M≤4). Inclusion of higher order moments leads to a slow decrease of the fit error (Fig. 2(b), M>1 ; Fig. 2(d), M>4) signifying an onset of noise-like edge roughness. After differentiating coarse variations from noise-like edge roughness we can now complete statistical description of feature edges by specifying the parameters of edge roughness.

2.3. Characterization of edge roughness

Model of disorder or roughness requires a model for the statistics of the correlation functions. While simple analytical forms of the correlation functions, such as exponential or Gaussian, are frequently used due to their integrability, we aim to derive the proper statistical forms directly from analysis of the images. Introduction of fractal dimensions allows us to develop a “family” of possible statistical distribution functions to describe roughness of features in PC lattices.

The most straightforward way of describing feature edge roughness is by considering statistical properties of a fit error function, which for a feature n is defined as

where n=[1,N_f ], i=[1,$N_{\mathit{\text{edge}}}^n$ ]. Values of ${\delta }_n^M$ (θ) for θ≠θ_i are interpolated. In our studies we used “nearest neighbor”,“linear”, and “cubic” interpolation schemes with almost identical final values of statistical parameters. While interpolation of roughness, in general, is not trivial, in our case on a several pixel scale feature edges are continuous and relatively smooth curves as they correspond to physical boundaries, thus allowing us to perform interpolation in a consistent way. In what follows, for a feature n we consider interpolated values ${\delta }_n^M$ (θ) on a uniform mesh θ=[0,2π-2π/$N_{\mathit{\text{edge}}}^n$ ] with $N_{\mathit{\text{edge}}}^n$ points, which allows a straightforward use of FFT transforms.

We now consider spectral and fractal properties of roughness, which present an alternative description to the angular momenta parameter approach of the preceding section. Fractal curves are scale-invariant structures, having a similar shape independent of the scale of observation. In practice, fractal stability should cover at least two decades in order to be unambiguously identified [20]. In the case of 2D PCs, even for the images with the highest sub-nanometer resolution, self-similar behavior of edge roughness extended only for a maximum of two decades in spatial wavelength. Nevertheless, fractal methodology seems to be useful for our purposes as fractal exponents inferred from spectral and fractal analysis are consistent with each other.

To introduce fractal dimension we consider Lipschitz function f (θ) having the property

where exponent H is called Lipschitz-Holder or Hurst exponent. When H=0, f (θ) is discontinuous, while if H=1, f (θ) is differentiable. For 0<H<1, function f (θ) is continuous but not differentiable and is known as fractal.

In order to perform a fractal analysis a “height to height” correlation function is introduced. For individual features it is defined as

Assuming that ${\delta }_n^M$ (θ) is a fractal curve with Hurst exponent H, from definition (10) it follows that when λ→0, $C_n^M$ (λ) ∝ λ^2H . By explicit squaring of an integrand in (10) and after minor manipulations we write

where autocorrelation function is defined as

For large enough values of λ that exceed noise correlation length λ>${\lambda }_{\mathit{\text{nc}}}^M$ , ${\mathrm{\Gamma }}_n^M$ (λ)→0, and consequently, $C_n^M$ (λ)→2${\mathrm{\sigma }}_n^2$ (M). In summary, asymptotics of spectral functions are

We find that for all the images of 2D PCs analyzed the following parameterizations (consistent with asymptotics (13)) of functions $C_n^M$ (λ) and ${\mathrm{\Gamma }}_n^M$ (λ) can be used to describe statistics of noise-like edge deviation from a smooth fit

 

Fig. 3. (a) “Height to height” correlation function and (b) auto-correlation function of an edge deviation from smooth fits with M angular components.

Download Full Size | PDF

These parameterizations work especially well when λ≲${\lambda }_{\mathit{\text{nc}}}^M$ , while for λ≲${\lambda }_{\mathit{\text{nc}}}^M$ oscillating features persist due to aliasing effects. Note, from their definitions on a periodic domain [C,Γ](λ)=[C,Γ](-λ)=[C,Γ](2${\pi R}_{\mathit{0}}^n$ -λ), therefore we will only consider 0≤λ≤${\pi R}_{\mathit{0}}^n$ . Moreover, for asymptotics (13) to hold correlation length of fit error has to be ${\lambda }_{\mathit{\text{nc}}}^M$ ≪${\pi R}_{\mathit{0}}^n$ .

In Fig. 3(a) “height to height” correlation function of an edge deviation from a smooth fit with M angular components is presented as a function of a spatial wavelengthλ. From (13) it follows that we can extract Hurst exponent of edge roughness by fitting a straight line to $C_n^M$ (λ) plotted on a log-log scale when λ≲${\lambda }_{\mathit{\text{nc}}}^M$ . For wavelengths larger than correlation wavelength λ≲${\mathrm{\lambda }}_{\mathit{\text{nc}}}^M$ , $C_n^M$ (λ) approaches a constant value. To get a reliable fit of Hurst exponent one typically needs fractal behavior to persist over several decades of λ. In all the high resolution images that we analyzed, fractal behavior persisted over one to two decades, thus making determination of Hurst exponents from scaling of$C_n^M$ (λ) somewhat imprecise. Thus, for example, from Fig. 3(a) Hurst exponent for M=1 curve is H=0.5 when 2nm<λ<20nm interval is considered, while H=0.43 when curve is fitted over the 2nm<λ<90nm interval. The upper value of this interval, however exceeds correlation length ${\lambda }_c^{\mathit{1}}$ =35nm and the fit underestimates the value of a Hurst exponent. As we pointed out earlier, on a smallest scale 0.46nm<λ<2nm our description of an edge as a single valued curve of θ is not valid any longer (see Fig. 1(c)) and this region can not be used in a fit. Note, that Hurst exponent of the remaining roughness is almost constant for the values M=(1,2,4,8) somewhat decreasing from H=0.5 to H=0.45 for largerM’s. Correlation length can be determined from Fig. 3(a) using parameterization (14) from which it follows that $C_n^M$ (λ$M_c^n$ )=2${\sigma }_n^{\mathit{2}}$ (M)(e-1)/e. Thus, for M=1,2,4,8 the values of correlation lengths are ${\lambda }_{\mathit{\text{nc}}}^M$ =35nm,22nm,11.2nm,6.4nm. We notice that ${\lambda }_{\mathit{\text{nc}}}^M$ is a decreasing function of the number M of angular components in a fit.

In Fig. 3(b) auto-correlation function of edge deviation from a smooth fit with M angular components is presented as a function of spatial wavelength λ. Correlation length can be also determined from Fig. 3(b) using parameterization (14) from which it follows that ${\mathrm{\Gamma }}_n^M$ (${\lambda }_{\mathit{\text{nc}}}^M$ )=${\mathrm{\sigma }}_n^{\mathit{2}}$ (M)/e. Both “height to height” and autocorrelation functions give very similar values of correlation lengths. To demonstrate that we plot in dotted lines parameterizations of autocorrelation function (14) with H=0.43,0.5, and correlation lengths deduced from asymptotics of “height to height” correlation function, and observe a good fit. Remaining oscillatory features for λ≫${\lambda }_{\mathit{\text{nc}}}^M$ are due to aliasing effectss.

An alternativeway of extracting Hurst exponents of a fractal data is using spectral techniques. First, we use a small number of angular components M to fit feature center coordinates $X_{\mathit{0}}^n$${\mathit{,}Y}_{\mathit{0}}^n$ and radius $R_{\mathit{0}}^n$ by minimizing edge objective function (1). These parameters converge rapidly as M increases, and we found that M=8 gave a reliable fit for all the analyzed images. Given feature center coordinates and a radius we consider deviation ${\delta }_n^{\mathit{1}}$ (θ) of an edge from a circle, interpolate it onto a uniform grid with the same number of $N_{\mathit{\text{edge}}}^n$ points as in a discretized image of an edge θ=(0:2π/$N_{\mathit{\text{edge}}}^n$ : 2π (1-1/$N_{\mathit{\text{edge}}}^n$ ), and use Fourier representation (2)

where coefficients A_m and B_m can now be efficiently computed using standard FFT. Substituting expansion (15) into (12) we get the following expression of the autocorrelation function

Next, we define power spectral density function $S_n^{\mathit{1}}$(λ_m) of edge deviation from a circle as

where λ_m =2${\pi R}_{\mathit{0}}^n$/m. Alternatively, using parameterization (14) and performing integration (17) we arrive to the following scaling relation

from which one can extract Hurst exponent by plotting S_n (λ_m ) versus λ_m on a log-log scale.

Another spectral method that can be used to find Hurst exponent of a fractal data involves plotting on a log-log scale RMS of edge deviation from a smooth fit σ_n(M) as a function of the number of angular components M in a fit (same plot as Fig.2(b,d) but for all the values of M=(2:$N_{\mathit{\text{edge}}}^n$ )). In principle, to calculate σ_n(M) we have to first solve minimization problem (1) that includes (1+2M) fit parameters, which for even moderate values of M>10 becomes time consuming. A considerably faster way to evaluate σ_n(M) even for large M is to assume that coefficients $A_m^M$${\mathit{,}B}_m^M$ in expansion (2) are independent of M. Then, in the same way as in calculating power spectral density we first find all the expansion coefficients in (15) using FFT. Then, ${\delta }_n^M$ (θ) can be expressed via expansion coefficients (2) as

Taking into account scaling (18) we find that for M large enough the following holds

Finally, when several features are present in the image we extract Hurst exponents from the averaged statistical functions (4).

In solid blue lines in Fig. 4(a,b) we present power spectral density and RMS of fit error for the same PC as in Fig. 1(a). Spectral density in Fig. 4(a) can be well fitted by a straight line over 2 decades starting from the largest length scale justifying our earlier finding that features in this PC can be characterized simply by an average radius and remaining roughness around a circular fit. Because only one feature is considered in this image the data is somewhat noisy and a range of Hurst exponents H_S =0.45-0.6 is possible. RMS of fit error in Fig. 4(b) can be fitted by a straight line starting from the lowest angular momentum M=1 all the way to M=300 with a Hurst exponent H_σ =0.45. In solid red lines we present spectral functions of estimated discretization noise due to finite image resolution (for more discussion see section 4). In solid blue lines in Fig. 4(c,d) we present power spectral density and RMS of fit error for the same PC as in Fig. 2(c). Power spectral density in Fig. 4(c) can be well fitted by a straight line over 1 decade in the interval 30nm≲λ≲200nm giving an estimate of the Hurst exponent H_S =0.3. RMS of fit error in Fig. 4(d) can be fitted by a straight line for angular momenta M > 4 all the way to M=40 with a Hurst exponent H_σ =0.28. As it was established earlier, ellipticity and a quadruple component are important contributions in the shape of the features, which is also clearly visible from Fig. 2(d), where power dependence of a spectral function is clearly observed only for M>4.

 

Fig. 4. (a) Power spectral density (blue). Linear fit is over 2 decades starting from the largest length scale. (b) RMS of a fit error (blue). Linear fit spans the lowest angular momenta starting with M=1. (c) Power spectral density (blue). Linear fit is over 1 decade in the interval 30nm≳λ≳200nm (d) RMS of a fit error (blue). Linear fit is in the range 4<M≲40. In red are the statistical functions of a noise level due to finite resolution of an image.

Download Full Size | PDF

Following discussions of this section we now present in Tables 1,2 several parameterizations of features in PCs considered in Figs. 1(a), 2(c). Parameters R_av, δ_el, θ_el, A_m, B_m -the radius of an underlying circle, ellipticity and other higher order moments of importance, are the coarse parameters of feature shapes averaged over features, whileσ (M), ${\lambda }_c^M$ , HC,S,σ describe statistical properties of the remaining roughness (edge deviation from a coarse fit).

3. Statistical description of feature lattices

In this section, we investigate variations in the feature positions from an ideal periodic lattice.

In Fig. 5(a) we present an image with resolution 5.54nm containing a PC lattice of 437 holes and a waveguide made of two rows of missing holes [22]. At first, coordinates of the hole centers ${\overline{r}}_{\mathit{0}}^n$ =($X_{\mathit{0}}^n$${\mathit{,}Y}_{\mathit{0}}^n$ ) are found by minimizing objective function (2) for various values of M. It was found that statistics of deviations of the hole centers from an underlying perfect lattice is not sensitive to a particular choice of M, and in what follows we choose M=3. Parameters of an underlying perfect lattice are then found by minimizing lattice objective function

Table 1. Parameterization of features in Fig. 1(a), InP/InGaAsP/InP [22].

View Table | View all tables in this article

Table 2. Parameterization of features in Fig. 2(c), Air/Si membranes [26].

View Table | View all tables in this article

 

Fig. 5. (a)PC lattice of holes with 2 missing rows. Vertices of a fitted perfectly periodic underlying lattice are shows as white dots. (b)PDDs of hole center deviations from the vertices of a perfect lattice along 2 principal directions (solids) together with Gaussian fits (dotted lines): perpendicular to the waveguide σ1 (blue), and parallel to the waveguide σ₂ (red). (c) RMS deviations σ1,2 (along 2 principle directions) of hole centers from an underlying lattice against the number of features in a fit. Features in a fit are included one by one, row by row starting from the upper left corner of an image.

Download Full Size | PDF

where ā_1,2 are the basis vectors of an underlying perfect lattice, and $j_{\mathit{1}\mathit{,}\mathit{2}}^n$ are the integer lattice coordinates of an n’s hole center. It is relatively straightforward to find 2N _f integer coordinates in (21) analytically given a reasonable approximation to the basis vectors, thus leaving 4 continuous parameters in ā_1,2 to be fitted. As before, we perform a fit with multidimensional Newton method by minimizing the value of a fit variance function Q_lat . In Fig. 5(a) vertices of a fitted perfectly periodic underlying lattice are shows as white dots.

We now define a 2D random variable $\overline{\delta }$_c=${\overline{r}}_{\mathit{0}}^n$${\mathit{-}j}_{\mathit{1}}^n$ā₁${\mathit{-}j}_{\mathit{2}}^n$ā₂ which we assume to be 2D Gaussian distributed with PDD

 

Fig. 6. (a) Uniform square PC lattice [26]. (b) σ_1,2 as a function of the number of features in a fit. Distribution of feature centers around the vertices of an underlying perfect lattice is isotropic. (c) Triangular PC lattice with a waveguide and a bend [28]. (d) σ_1,2 as a function of the number of features in a fit. Distribution of feature centers around the vertices of an underlying perfect lattice is anisotropic.

Download Full Size | PDF

where $R=\left(\begin{array}{cc}\cos \left(\theta \right)& \sin \left(\theta \right)\\ -\sin \left(\theta \right)& \cos \left(\theta \right)\end{array}\right)$ is a 2D rotation matrix, and σ_1,2 are the variances along the two principle directions. One can deduce statistical parameters σ_1,2,θ by using the following averages of a 2D Gaussian random variable, <${\delta }_c^x$${\delta }_c^x$ >=${\sigma }_{\mathit{1}}^{\mathit{2}}$ cos² (θ)+${\sigma }_{\mathit{2}}^{\mathit{2}}$ sin² (θ), <${\delta }_c^y$ ${\delta }_c^y$ >=${\sigma }_{\mathit{1}}^{\mathit{2}}$ sin² (θ)+${\sigma }_{\mathit{2}}^{\mathit{2}}$ cos² (θ), <${\delta }_c^x$${\delta }_c^y$ >=2cos(θ)sin(θ)(${\sigma }_{\mathit{1}}^{\mathit{2}}$${\mathit{-}\sigma }_{\mathit{2}}^{\mathit{2}}$ ).

In Fig. 5(b) we plot PDD of $\overline{\delta }$_c along the two principle directions (θ=-1.6°) from the lattice fit (solid lines) and a corresponding Gaussian distribution (dotted lines). We find that a 2D distribution of feature center displacements from the vertices of an underlying perfect lattice indeed appears to be Gaussian and is highly anisotropic. The RMS of the hole center deviations from a perfect lattice is twice as large σ₁=6.4nm in the direction perpendicular to the waveguide than in the parallel direction σ₂=2.9nm. In Fig. 5(c) we investigate in more details the source of such an anisotropy. RMS deviations σ_1,2 (along 2 principle directions) of hole centers from an underlying lattice are plotted against the number of features in a fit. Features in consecutive fits are added one by one, row by row starting from the upper left corner of a structure Fig. 5(a). At leat two rows are needed to fit both basis vectors. When only few features are included in a fit, parameters σ_1,2 grow rapidly with each included feature, finally “saturating” when 50 features are included (2 rows). We observe, quite generally, that at least 10–20 features in each row are needed to determine parameters of a Gaussian fit reliably. When more then 5 layers are included in the fit, and while approaching a waveguide region, hole center deviations from an underlying lattice become highly anisotropic. As we have mentioned earlier, such an anisotropy can appear for non rectangular lattices. From analysis of various PC images we observe that anisotropy in the deviations of feature centers from an underlying perfect lattice is predominantly observed in the structures with symmetry breaking elements such as waveguides, bends, etc. For PC lattices with waveguides, for example, we find that frequentlyσ₁>σ₂, where σ₁ is RMS of feature center deviations in the direction perpendicular to a waveguide, while σ₂ is RMS of deviations in the direction parallel to a waveguide. Uniform rectangular PC lattices without functional elements are typically isotropic σ₁~σ₂. Most likely physical reason for such an anisotropy being varying along a non-uniform PC lattice e-beam proximity effects due to non-uniform local environments (see Appendix).

Image with resolution 1.63nm of a uniform square PC lattice containing 204 holes [26] is presented in Fig. 6(a). Dependence of RMS parameters σ ₁, σ₂ along two principal directions for increasing number of features in a fit is plotted in Fig. 6(b). Features in consecutive fits are added one by one, column by column starting from the upper left corner of a structure. One observes that σ₁~σ₂~1.6nm for any number of features in a fit, “saturating” to their stationary values after two rows (~20 features) are included. In Table 3 we present a complete statistical model of feature shapes and feature center distribution from ideal lattice for Fig. 6(a).

Table 3. Parameterization of features in Fig. 6(a), Air/Si membranes [26].

View Table | View all tables in this article

In Fig. 6(d) we present dependence of RMS parameters σ₁, σ₂ along the two principal directions for increasing number of features in a fit for a waveguide and a bend in a PC lattice of holes [28]. Features in consecutive fits are added one by one, row by row starting from the upper left corner of a structure Fig. 6(c). One observes thatσ₁~σ₂~3nm when first three rows of a structure (N_f <60) are included in the fit. Row number 4 is closest to a waveguide from the top and is made of the holes with sizes somewhat smaller than the bulk ones. When including waveguide edge row in a fit (60<N_f <80) first sign of anisotropy appears with σ₁<4nm>σ₂~3nm. Next twenty features introduce a bend into the structure (80<N_f <100) making anisotropy even stronger σ₁~5.5nm>σ₂~4nm. Finally, when the bulk of a PC lattice on the other side of a waveguide is added into the fit, values of RMS of hole center variations from an ideal lattice “saturate” to the values σ₁~5nm>σ₂~4nm, θ=23°. From Fig. 6(d) we also notice that RMS of distribution of feature centers around vertices of an underlying perfect lattice is location dependent. Thus, far from any structural element (such as waveguide or a bend) σ₁~σ₂, while for the rows bordering a waveguide σ₁>σ₂. Typically, PC regions directly bordering a waveguide are the ones determining scattering and absorption losses of radiation propagating through a waveguide, as light penetration into a PC lattice is limited to few periods. Thus, when modelling effects of waveguide non-uniformity on radiation loss one has to pay special care to derive a realistic statistical model of such non-uniformities in a region closest to a waveguide.

4. Discretized images and associated errors

In our edge detection algorithm, we first sort all the image pixels into two categories the ones that belong inside of a hole and the ones that belong to a substrate. We work with normalized grayscale images where the value of each pixel pix is between 0 (black-hole) and 1 (white-substrate). First, a 3–5 pixel convolution (“smoothing”) is applied to each image to reduce noise. To sort the pixels we compare their values to a threshold parameter tol, if pix<tol we consider a pixel to be in the hole and assign it a value 0, while in the opposite case we assign it a value 1. As different images have different contrast and noise levels to find a reasonable value of a tol parameter we use a histogram of pixel values. In Fig. 7(a) we show an image with resolution res=1.25nm of a hole [24] and a histogram of pixel values in the insert. This image features a moderate index contrast (pixel value ratio for the two maxima corresponding to white and black is 0.6:0.1) and a relatively high noise (region 0.2-0.4 of substantial number of pixels with values in between two maxima). For edge detection we try several values of a threshold parameter tol=0.37,0.40,0.43 around the local minima of a histogram in between two maxima (see insert). Once all the pixels are sorted into two groups (with values 0 or 1) hole edge is detected using a standard method of image convolution with a Sobel-like 3×3 matrix. In Figs. 7(b,c,d) we present edges detected by using three different threshold parameters. As seen from the images, detected edges are somewhat different from each other.

 

Fig. 7. (a) Image of a hole with a moderate contrast and a high noise level. Insert: histogram of pixel values. Hole edges are detected with: (b) tol=0.37 (c) tol=0.40 (d) tol=0.43

Download Full Size | PDF

An important question is how sensitive are the statistical parameters characterizing feature imperfections with respect to variations in tol. Detailed simulations show that for a Δtol/tol~15% change in a threshold parameter, statistical parameters only vary as ΔR_av/R_av ~2%, Δδ_el /δ_el ~12%, Δσ (2)/σ (2)<1%, ΔH_S,σ /H_S,σ ~15%, Δσ_1,2/σ_1,2~5%. Note that for this image, the parameter affected most by variations in tol is a Hurst exponent of a remaining wall roughness. This is generally expected in the case of a noisy image where detected roughness on the smallest case (several nm) is strongly affected by image noise, and hence a value of a tol parameter. Thus, when the same sensitivity analysis is repeated for a higher resolution, lower noise image Fig. 1(a) we find that for a~10% change in a tol parameter Hurst exponents changes only by a few percents. We conclude that while there is indeed an uncertainty associated with a choice of a threshold parameter, the resultant values of the statistical parameters are weakly sensitive functions of tol given a good quality image. Moreover, these uncertainties can be further reduced by lowering the image noise level and increasing contrast.

Next, we estimate the level of discretization noise due to finite resolution of an image (red curves in Fig. 4). Particularly, given an image we first fit its coarse parameters R ₀,X ₀,Y ₀,A_m,B_m for a certain small number of angular momenta m=[2,M], M=3-8. Then, we consider an analytical curve $r_{\mathit{\text{fit}}}^M$ (θ) with thus found fitted parameters (2) centered around X ₀=0,Y ₀=0. Next, we introduce a uniform mesh with the same resolution as that of an analyzed image and project an analytical edge onto a discrete mesh to get a discretized approximation of an edge $r_{\mathit{\text{mesh}}}^M$ (θ). Then, we take a difference between an analytical curve and its discretized version δ_noise (θ_i )=($r_{\mathit{\text{mesh}}}^M$ (θ_i )-$r_{\mathit{\text{fit}}}^M$ (θ_i )) to estimate the level and statistics of a nose due to image discretization (red curves in Fig. 4). For all the images analyzed we find that RMS of δ_noise (θ) is on the order of 0.4·res in the worst case, contributing most to the uncertainties in a wall roughness parameter σ (M). Moreover, when a discretized edge is fitted by minimizing (1), the error in the coarse parameters due to discretization is even smaller. Thus, when comparing the center coordinates $X_{\mathit{0}}^{\mathit{\text{mesh}}}$${\mathit{,}Y}_{\mathit{0}}^{\mathit{\text{mesh}}}$ by fitting a discretized edge with a circle and the exact center coordinates X₀=0,Y₀=0, we find that their maximum discrepancy over all the images and features did not exceed 0.1·res. Thus, uncertainties in the parameters σ_1,2 of RMS deviation of feature centers from a perfect lattice are also at most 0.1·res.

5. Discussion and conclusions

We find that at least three sets of parameters are necessary to create a minimal statistical model of 2D disorder in PC lattices. First set of parameters describes coarse properties of individual shapes persistent over all features such as radius, ellipticity and other low angular momenta, among which radius is the most important. Another set of parameters describes higher angular momenta plus random edge roughness by a set of correlation functions (14) with parameters ${\lambda }_c^M$ , σ (M) and H corresponding to correlation length, standard deviation and Hurst exponent. Typically, unless written deliberately, we find that even low angular momenta components (such as ellipticity) are not persistent from one feature to another and can be simply described as part of a random edge roughness. A final set of parameters describes deviations of feature centers from an ideal periodic lattice in terms of a 2D Gaussian distribution parameterized by two principal directions and two variances σ_1,2 along such directions. For PC lattices with symmetry breaking elements such as waveguides, bends, etc. we find that due to non-uniform e-beam proximity effects feature position disorder is frequently anisotropic σ₁≫σ₂.

Findings of this paper are based on the analysis of over 30 high resolution pictures of the “typical” e-beam written structures of various material combinations. In Table 4 we present the values (within 2 standard deviations from the averages) of various statistical parameters averaged over all analyzed images with resolutions 0.46nm-6nm. A somewhat surprising finding is that despite all the different material combinations from which these PC lattices are made a relatively narrow distribution of statistical parameters characterizing disorder is found.

Table 4. Ranges of statistical parameters over various PC lattices [21].

View Table | View all tables in this article