Statistical analysis of geometrical imperfections from the images of 2 D photonic crystals

High resolution images of planar photonic crystal (PC) optical components fabricated by e-beam lithography in various materials are analyzed to characterize statistical properties of common 2D geometrical imperfections. Our motivation is to attempt an intuitive, while rigorous statistical description of fabrication imperfections to provide a realistic input into theoretical modelling of PC device performance. © 2005 Optical Society of America OCIS codes: 130.3120 Integrated optics devices, 030.5770 Roughness, 100.2960 Image analysis, 100.5010 Pattern recognition and feature extraction References and links 1. K.Y. Bliokh, Y.P. Bliokh, and V. D. Freilikher, “Resonances in one-dimensional disordered systems: localization of energy and resonant transmission,” J. Opt. Soc. Am. B 21, 113-120 (2004). 2. V.M. Apalkov, M.E. Raikh, and B. Shapiro, “Almost localized photon modes in continuous and discrete models of disordered media,” J. Opt. Soc. Am. B 21, 132-140 (2004). 3. E. Lidorikis, M.M. Sigalas, et al. “Gap deformation and classical wave localization in disordered two-dimensional photonic-band-gap materials,” Phys. Rev. B 61, 13458-13464 (2000). 4. A.A. Asatryan, P.A. Robinson, et al. “Effects of geometric and refractive index disorder on wave propagation in two-dimensional photonic crystals,” Phys. Rev. E 62, 5711-5720 (2000). 5. M.A. Kaliteevski, J.M. Martinez, et al. “Disorder-induced modification of the transmission of light in a twodimensional photonic crystal, ” Phys. Rev. B 66, 113101 (2002). 6. K.-C. Kwan, X. Zhang, et al. “Effects due to disorder on photonic crystal-based waveguides,” Appl. Phys. Lett. 82, 4414-4416 (2003). 7. B.Z. Steinberg, A. Boag, and R. Lisitsin, “Sensitivity analysis of narrowband photonic crystal filters and waveguides to structure variations and inaccuracy,” J. Opt. Soc. Am. A 20, 138 (2003). 8. B.C. Guptaa, and Z. Yeb, “Disorder effects on the imaging of a negative refractive lens made by arrays of dielectric cylinders,” J. Appl. Phys. 94, 2173-2176 (2003). 9. S. Lan, K. Kanamoto, et al. “Similar role of waveguide bends in photonic crystal circuits and disordered defects in coupled cavity waveguides: An intrinsic problem in realizing photonic crystal circuits,” Phys. Rev. B 67, 115208 (2003). 10. T. N. Langtry, A.A. Asatryan, et al. “Effects of disorder in two-dimensional photonic crystal waveguides,” Phys. Rev. E 68, 026611 (2003). 11. A.G. Martyn, D. Hermann, et al. “Defect computations in photonic crystals: a solid state theoretical approach,” Nanotechnology 14, 177183 (2003). 12. N.A. Mortensen, M.D. Nielsen, et al. “Small-core photonic crystal fibres with weakly disordered air-hole claddings,” J. Opt. A: Pure Appl. Opt. 6, 221223 (2004). (C) 2005 OSA 4 April 2005 / Vol. 13, No. 7 / OPTICS EXPRESS 2487 #6498 $15.00 US Received 31 January 2005; revised 17 March 2005; accepted 20 March 2005 13. M. Skorobogatiy, “Modelling the impact of imperfections in high index-contrast photonic waveguides,” Phys. Rev. E 70, 46609 (2004). 14. W. Bogaerts, P. Bienstman, and R. Baets, “Scattering at sidewall roughness in photonic crystal slabs,” Opt. Lett. 28, 689-691 (2003). 15. M.L. Povinelli, S.G. Johnson, et al. “Effect of a photonic band gap on scattering from waveguide disorder,” Appl. Phys. Lett. 84, 3639-3641 (2004). 16. S. Fan, P.R. Villeneuve, and J.D. Joannopoulos, “Theoretical investigation of fabrication-related disorder on the properties of photonic crystals,” J. Appl. Phys. 78, 1415-1418 (1995). 17. V. Yannopapas, A. Modinos, and N. Stefanou, “Anderson localization of light in inverted opals,” Phys. Rev. B 68, 193205 (2003). 18. D.J. Whitehouse, “Surface Characterization and Roughness Measurement in Engineering,” Photomechanics, Topics Appl. Phys. 77, 413461 (2000). 19. D.J. Whitehouse, “Some theoretical aspects of structure functions, fractal parameters and related subjects,” Proc. Instn. Mech. Engrs. Part J 215, 207-210 (2001). 20. I. Arino, U. Kleist, et al. “Surface Texture Characterization of Injection-Molded Pigmented Plastics,” Polym. Eng. Sci. 44, 1615-1626 (2004). 21. Developed software and analysed PC images are available at http://www.photonics.phys.polymtl.ca/PolyFIT/ 22. A. Talneau, M. Mulot, et al. “Compound cavity measurement of transmission and reflection of a tapered singleline photonic-crystal waveguide,” Appl. Phys. Lett. 82, 2577-2579 (2003). 23. M.Mulot, S.Anand, et al. “Low-loss InP-based photonic crystal waveguides etched with Ar/Cl2 chemically assisted ion beam etching,” J. Vac. Sci. Technol. B21, 900-903 (2003). 24. A. Xing, M. Davanco, et al. “Fabrication of InP-based two-dimensional photonic crystal membrane,” J. Vac. Sci. Technol. B 22, 70-73 (2004). 25. C. Monat, C. Seassal, et al. “Two-dimensional hexagonal-shaped microcavities formed in a two-dimensional photonic crystal on an InP membrane,” J. App. Phys. 93, 23-31 (2003). 26. P.E. Barclay, K. Srinivasan, et al. “Efficient input and output fiber coupling to a photonic crystal waveguide,” Opt. Exp. 29, 697-699 (2004). 27. H. Altuga and J. Vuckovic, “Two-dimensional coupled photonic crystal resonator arrays,” Appl. Phys. Lett. 84, 161-163 (2004). 28. M. Augustin, H.-J. Fuchs, et al. “High transmission and single-mode operation in low-index-contrast photonic crystal waveguide devices,” Appl. Phys. Lett. 84, (2004).


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][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][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][4][5][6][7][8][9][10][11][12][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 2 /Nb 2 O 5 /SiO 2 [28].

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.

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 0 ,Y 0 ) to an edge point  2) is an m = 0 term, implicit m = 1 terms correspond to the feature center coordinates (X 0 ,Y 0 ), while m = 2 expansion coefficients in (2) define feature ellipticity as 2 with an angle of a major axis defined by cos(2θ el ) = B M 2 /δ 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 σ 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 M R 0 (nm) X 0 (nm) Y 0 (nm) δ el (nm) θ el (Deg) σ (M)(nm) Note that, in general, assumption that feature edge r edge (θ i ) can be fitted with a single valued analytical curve r M f it (θ ) 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 n edge 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 n 0 , σ 2 n (M), [C, Γ, S] M n (λ ) are the radius, variance and correlation functions of a feature n, then their averaged counterparts are defined as

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 Note that < δ M r >= 0 as it is proportional to the derivative of (1) with respect to a feature radius.Variance of δ M r is given by (4).In what follows, we find that for large enough values of M distribution of the remaining edge roughness δ M r can always be fitted by a Gaussian probability density distribution 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 δ M r 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 σ 2 (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 δ M r 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 δ M r 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 σ 2 (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 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.
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.

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 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 Values of δ M n (θ ) 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 δ M n (θ ) on a uniform mesh θ = [0, 2π −2π/N n edge ] with N n edge 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 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 δ M n (θ ) is a fractal curve with Hurst exponent H, from definition (10) it follows that when λ → 0, C M n (λ ) ∝ λ 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 λ > λ M nc , Γ M n (λ ) → 0, and consequently, C M n (λ ) → 2σ 2 n (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 M n (λ ) and Γ M n (λ ) can be used to describe statistics of noise-like edge deviation from a smooth fit  These parameterizations work especially well when λ ∼ < λ M nc , while for λ ∼ > λ M nc oscillating features persist due to aliasing effects.Note, from their definitions on a periodic domain , therefore we will only consider 0 ≤ λ ≤ πR n 0 .Moreover, for asymptotics (13) to hold correlation length of fit error has to be λ M nc πR n 0 .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 M n (λ ) plotted on a log-log scale when λ ∼ < λ M nc .For wavelengths larger than correlation wavelength λ ∼ > λ M nc , C M n (λ ) 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 M n (λ ) 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 λ 1 c = 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 larger M's.Correlation length can be determined from Fig. 3(a) using parameterization (14) from which it follows that C M n (λ M nc ) = 2σ 2 n (M)(e − 1)/e.Thus, for M = 1, 2, 4, 8 the values of correlation lengths are λ M nc = 35nm, 22nm, 11.2nm, 6.4nm.We notice that λ M nc 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 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 λ λ M nc are due to aliasing effectss.An alternative way 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 n 0 ,Y n and radius R n 0 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 δ 1 n (θ ) of an edge from a circle, interpolate it onto a uniform grid with the same number of N n edge points as in a discretized image of an edge θ = (0 : 2π/N n edge : 2π(1 − 1/N n edge ), 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 1 n (λ m ) of edge deviation from a circle as where λ m = 2πR n 0 /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 n edge )).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 , 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, δ M n (θ ) 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  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. 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), λ M c , H C,S,σ describe statistical properties of the remaining roughness (edge deviation from a coarse fit).

Statistical description of feature lattices
In this section, we investigate variations in the feature positions from an ideal periodic lattice.
In Fig. 5 R av (nm) δ el (nm) θ el ( o ) A 3 (nm) B 3 (nm) A 4 (nm) B 4 (nm) 178.5 ± 7.8 22.1 ± 1.9 87.3 ± 1.5 0.6 ± 1.1 −0.8 ± 1.1 −0.5 ± 0.9 3.6 ± 0.9 σ (4)(nm) λ 4 c (nm) H S H σ 1.67 19.1 0.3 0.28   and a waveguide made of two rows of missing holes [22].At first, coordinates of the hole centers rn 0 = (X n 0 ,Y n 0 ) 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 where ā1,2 are the basis vectors of an underlying perfect lattice, and j n 1,2 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 δc = rn 0 − j n 1 ā1 − j n 2 ā2 which we assume to be 2D Gaussian distributed with PDD where R = cos(θ ) sin(θ ) −sin(θ ) cos(θ ) 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, < δ ).In Fig. 5(b) we plot PDD of δc along the two principle directions (θ = −1.6 o ) 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 σ 1 = 6.4nm in the direction perpendicular to the waveguide than in the parallel direction σ 2 = 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 σ 1 > σ 2 , where σ 1 is RMS of feature center deviations in the direction perpendicular to a waveguide, while σ 2 is RMS of deviations in the direction parallel to a waveguide.Uniform rectangular PC lattices without functional elements are typically isotropic σ 1 ∼ σ 2 .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 σ 1 , σ 2 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 ∼ σ 2 ∼ 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).
In Fig. 6(d) we present dependence of RMS parameters σ 1 , σ 2 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 σ 1 ∼ σ 2 ∼ 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 σ 1 ∼ 4nm > σ 2 ∼ 3nm.Next twenty features introduce a bend into the structure (80 < N f < 100) making anisotropy even stronger σ 1 ∼ 5.5nm > σ 2 ∼ 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 σ 1 ∼ 5nm > σ 2 ∼ 4nm, θ = 23 o .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) σ 1 ∼ σ 2 , while for the rows bordering a waveguide σ 1 > σ 2 .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.

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 3x3 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.
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 0 , X 0 ,Y 0 , A m , B m for a certain small number of angular momenta m = analytical curve r M f it (θ ) with thus found fitted parameters (2) centered around X 0 = 0,Y 0 = 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 M mesh (θ ).Then, we take a difference between an analytical curve and its discretized version δ noise (θ i ) = (r M mesh (θ i ) − r M f it (θ 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 mesh 0 ,Y mesh 0 by fitting a discretized edge with a circle and the exact center coordinates X 0 = 0,Y 0 = 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.

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 λ M c , σ (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 σ 1 σ 2 .
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.

Appendix
We demonstrate physical processes responsible for disorder in PC lattices on the example of a particular fabrication process using direct e-beam writing in InP/InGaAsP/InP materials [22].
One starts with a semiconductor multilayer where top InP layer is 200nm thick, followed by a 450nm thick optically guiding GaInAsP layer, then a 2µm thick buffer InP layer on the bottom, and, finally, an InP substrate.Processing steps are as follows: first, on the top of an InP layer one deposits a 250nm thick SiO 2 layer, then on the top of a SiO 2 layer one deposits a 300nm thick PMMA polymer layer.After that, circular features are developed in PMMA with e-beam writing.The short penetration length of electrons precludes the use of a solid SiO 2 substrate as a mask directly.Next, SiO 2 layer is dry etched with CHF 3 using the PMMA layer as a mask.Finally, PMMA is removed and SiO 2 layer is used as a mask for chemically assisted ion beam etching [23] of holes.Resulting holes are typically etched 3 − 4µm deep.The resultant structure is the etched PC lattice where roughness presents an accumulated effect of several fabrication steps: PMMA development by e-beam writing, SiO 2 etching and semiconductor multilayer etching.Overall, it seems that the resultant roughness is less a function of a particular material combination but rather the details of a fabrication process.For example, an alternative process to create a PC lattice resulting in higher roughness would be a so-called lift-off process.In this process one deposits a PMMA layer directly on the top of an InP layer.Then, the complementary of holes are developed by e-beam writing and a continuous metallic layer is deposited.When washed by a solvent remaining PMMA dissolves by leaving a metallic mask for the holes while somewhat tearing the metallic mask layer near the hole edges.
A typical e-beam writing strategy for planar PC lattices is direct writing : a beam of electrons of a given diameter, moves pixel by pixel in x-y directions with a step size as small as 2.5nm.A feature boundary is typically coded as a polygon of many vertices (18 in this case).This polygon is then subdivided by software into elementary shapes such as rectangles and triangles for further exposure.If lattice is not rectangular than it becomes impossible to resolve exactly non-integer coordinates of the feature centers, which could introduce larger deviations of feature centers from an ideally periodic lattice along certain spatial directions.As the electrons penetrate into the resist material a considerable number of them experience large angle scattering leading to backscattering, thus causing additional exposure in the resist and what is called the electron beam proximity effect.Roughness introduced by e-beam proximity effects in the PMMA resist is theoretically estimated to be on the order of 3 − 10nm, however the measured roughness is typically smaller.If no software to compensate for proximity effects is used then exposure conditions for points in different local environments will be different.This can result in measurable distortions for the features located near the symmetry breaking features such as corners, waveguides, bends, resonators, etc. by comparison with features inside of a uniform periodic lattice.Finally, distortions in the shape of a feature (such as ellipticity of a hole) could also come from a non optimal setting of an imaging SEM in a form of a "residual astigmatism", which could also make roughness look erroneously larger along certain spatial directions.

Fig. 1 .
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 M f it (θ ).(d) Edge roughness is self-similar on very different scales suggesting fractal description.

Fig. 2 .
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.
(C) 2005 OSA 4 April 2005 / Vol. 13, No. 7 / OPTICS EXPRESS 2492 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.
Fig. 3. (a) "Height to height" correlation function and (b) auto-correlation function of an edge deviation from smooth fits with M angular components.

Fig. 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.

Fig. 5 .
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 σ 2 (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.

Fig. 6 .
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.