Reflection of focused beams from opal photonic crystals

We present a robust method for computing the reflection of arbitrarily shaped and sized beams from finite thickness photonic crystals. The method is based on dividing the incident beam into plane waves, each of which can be solved individually using Bloch periodic boundary conditions. This procedure allows us to take a full advantage of the crystal symmetry and also leads to a linear scaling of the computation time with respect to the number of plane waves needed to expand the incident beam. The algorithm for computing the reflection of an individual plane wave is also reviewed. Finally, we find an excellent agreement between the computational results and measurement data obtained from opals that are synthesized using polystyrene and poly(methyl methacrylate) microspheres. © 2005 Optical Society of America OCIS codes: (000.3860) Mathematical methods in physics, (000.4430) Numerical approximation and analysis, (160.5470) Polymers

Photonic crystals have been simulated with various methods.Photonic dispersion curves are often solved using the plane wave method (PWM) [14][15][16] or the finite-difference frequencydomain method (FDFD) [17], whereas the density of photonic states and reflection coefficients can more readily be computed using some of the many implementations of the transfer matrix method (TMM) [18,19].The layer-multiple-scattering method [20] has proven to give very good prediction for three dimensional photonic crystals made of spherical particles, and other multiple-scattering methods have been successfully applied on 2D crystals [21,22].For problems where the frequency domain methods are not easily applicable, such as in cavity problems with leaky modes or transient phenomena, the finite-difference time-domain method (FDTD) is often the preferred choice [23,24].A review and comparison of existing computational methods can be found in Ref. [25].
Our approach for solving the electromagnetic fields in finite thickness photonic crystals is to divide the crystal slab into thin layers and expand the fields in them using plane waves.The plane wave expansions in adjacent layers are then related to each other through a finite difference scheme in computation.Fields above and below the slab are expanded in homogeneousmedium eigenvectors, which can be solved analytically and provides the boundary conditions of the computational domain.Finally, a system matrix is formed by requiring the continuity of the tangential fields on the interfaces between the three domains.Our approach has the following features: i) It leads to sparse matrices, which can be solved efficiently using iterative solvers, ii) the method is very stable, iii) the dielectric function can be dispersive and complex valued and is allowed to vary freely in all three coordinates inside the slab, and iv) the method can be used to solve for both excitation and eigenmode problems.
Our method borrows the plane wave basis from PWM and the finite difference technique from TMM and combines them in a way which is well suited to two-dimensionally periodic finite thickness structures.The biggest difference to TMM is that we do not use the finite difference scheme for eliminating the field variables inside the slab, but instead include them all explicitly into the system equation.This approach is completely free from the instabilities often encountered in TMM [18].The direct discretization technique of the dielectric function makes the implementation of our method simple.This is particularly pronounced in cases where the dielectric function cannot be divided into individual scatterers as required by the multiplescattering techniques, or where the scatterers are of arbitrary shape and the scattering matrices cannot be analytically solved.In this paper, we summarize the theoretical basis of the method and refer the reader to Ref. [26] for implementation details.
Opals are often characterized by measuring their reflection and/or transmission spectra.If the incident field is assumed to be a plane wave, the problem can be reduced to a single unit cell of the photonic crystal using Bloch periodic boundary conditions.However, to relate the calculated spectrum to the real-world experimental setups for these measurements, it is often desirable to analyze a focused beam in which case the Bloch periodic boundary conditions cannot be applied and the problem becomes significantly more complicated.Two approaches can be adopted: a large portion of the crystal can be considered as one entity, which is illuminated by the finite sized beam (see e.g.[27,28]) or the beam may be decomposed into its plane wave components, each of which are solved individually using Bloch periodic boundary conditions, thus considering only one unit cell of the lattice.In this paper we adopt the latter technique and compute the reflection from finite thickness opals using the presented computational method.The advantages of the decomposition technique include the following: i) Instead of solving a big problem once, one solves a small problem many times, resulting in a linear scaling of computational resources, ii) after the reflection coefficients for individual plane waves have been computed, the reflection coefficient of any arbitrary beam can be obtained with minor effort, iii) the decomposition technique parallelizes trivially since each of the incident plane waves creates an independent problem.

Preliminary considerations and definitions
We focus on problems with the following geometry: Region 1, z < 0, is homogeneous medium characterized by ε a .Region 2, 0 ≤ z < h, is a finite thickness photonic crystal defined by ε b (r, z) = ε b (r + ma 1 + na 2 , z), where a 1 and a 2 are lattice vectors on xy-plane and m and n are integers.Region 3, z ≥ h, is again homogeneous medium characterized by ε c .Dielectric functions in all three regions may be dispersive and complex valued.
In our terminology, a unit cell is defined by the two-dimensionally periodic volume defined by the planes z = 0 and z = h and the two transverse lattice vectors a 1 and a 2 .This should not be confused with the unit cell of a truly three-dimensionally periodic photonic crystal, defined by three lattice vectors.The semi-periodicity in the z-direction, as in finite thickness opals, can be obtained simply by filling the unit cell with a semi-periodic dielectric function.The time dependence is assumed to be harmonic f (t) ∼ exp(− jωt) throughout this paper and its explicit notation is omitted.

Diagonalized form
Maxwell's curl equations can be written in the following form [26]: which we refer to as being diagonalized with respect to z since all z-dependency is isolated to the right hand side of the equation [29].Here, Ψ = [ E x E y H x H y ] T is a vector consisting of the transversal field components and ∂ k denotes differentiation with respect to the given coordinate k = x, y, z.L is a matrix operator where the 2 × 2 sub-operator A is defined below: Here, j is the imaginary unit and parameters ε and µ are the (generally) complex valued and position dependent permittivity and permeability functions, respectively.The remaining suboperator B can be obtained from A by replacing ε with −µ and vice versa.The obvious consequence of Eq. ( 1) is that once the transversal field distribution in some z =constant plane is known, the derivative in the normal direction is easily evaluated and the transversal fields are then uniquely determined everywhere in the space.The two remaining field components are redundant and can be computed from the transversal fields a posteriori [26].

Field expansions
The periodicity on the xy-plane suggests expanding the fields in a plane wave basis where [G n ] is a truncated set of reciprocal lattice vectors and K is a Bloch wave vector restricted in the first Brillouin zone (BZ).The expansion coefficients Ψ n (z) are vectors of the z-dependent transversal field components, which will be expanded differently in the photonic crystal medium and the claddings.
In the photonic crystal slab, we discretize the fields and material parameters on z =constant layers and use finite differences for relating the fields in adjacent layers to each other.Substituting Eq. ( 4) into Eq.( 1) gives the z-derivative of the transversal fields and by using the standard first order finite difference scheme we get where ∆ is the finite difference step length and i is the index of the plane.A numerically useful formula is obtained by multiplying both sides of Eq. ( 5) with exp [− j (G m + K) • r] for each plane wave G m in the basis and integrating over the unit cell.The result is a matrix equation relating the plane wave expansion coefficients in adjacent layers to each other.In the actual implementation, we define electric fields in planes z = (i + 0.5)∆ and magnetic fields in z = i∆, for whole numbers i, since it follows from Eq. (1) that the derivative of electric field depends only on the magnetic field and vice versa.Equation ( 5) is then correspondingly separated into two parts.Fig. 1 illustrates the discretization scheme.
In the cladding regions 1 and 3, where ε is constant, we assume an exponential z-dependency: which allows an analytical solution for the eigenvectors and eigenvalues of Eq. ( 1).Principally, there are in these regions 4N eigenvalues for a set of N plane waves (due to the four-dimensional vector coefficients) but the system decouples into N independent eigensystems of dimension four due to the position independent material parameters.For a given K, the complete expression for the fields is where a ± K,n and b ± K,n are scalar coefficients.The two doubly degenerate eigenvalues of Eq. ( 1) A schematic representation of the discretization.The unit cell in the slab is limited by planes z = 0 and z = h and is periodic with a 1 and a 2 .Electric fields are expanded in plane waves on light blue planes and magnetic on dark brown.The plane wave expansions in adjacent planes are related to each other through finite differences.Outside the slab, for z < 0 and z > h, the fields are expanded in outgoing eigenvectors.The eigenvector expansion is used to terminate the finite difference grid by relating the electric field at z = −(1/2)∆ to the magnetic field at z = 0 and similarly for the other cladding.
where the positive (negative) sign corresponds to a wave which either propagates or decays in the direction of the positive (negative) z-axis.The corresponding four eigenvectors are where k x,y = u x,y • (K + G n ), with u x,y denoting the dimensionless unit vector in x-and ydirections.There is a great freedom in the selection of the eigenvectors since any linear combination of eigenvectors corresponding to the same eigenvalue is also an eigenvector.Our preferred choice is such that where E evaluates to the electric field part of the argument and H to the magnetic, u z is the unit vector in z-direction and the asterisk denotes complex conjugation.The normalization coefficients α K,n and β K,n are chosen such that where l = 1, 2 and ℜ denotes the real part.It should be noted that the conditions in Eqs.(10) and ( 11) are only valid for real-valued material parameters.

Constructing and solving the system equation
The final step in the method is to relate the field expansions in the three domains by requiring the continuity of the transverse field components at the interfaces.Only 2N of the total 4N eigenvectors are needed to satisfy the interface conditions for each cladding and we select those which either decay or radiate energy away from the slab.The result is a homogeneous system of equations, whose solutions correspond to guided slab modes in the photonic crystal.Sources such as currents or incident fields can be easily included by replacing the right hand zero vector with the source vector, to obtain an equation of the form where M is a matrix containing the eigenvectors and all the finite-difference relations in the slab volume, f is a vector containing all the expansion coefficients for the transverse field components and b describes the excitation.
The finite difference scheme in Eq. ( 5) is similar to some transfer matrix methods [18] but there are two differences: i) We apply it directly to the plane wave expansions instead of the real space fields and ii) we do not use it recursively for the elimination of the the field variables.Instead we include all the field expansion coefficients throughout the slab volume explicitly as unknowns in f.This means that if each of the four field components is expanded in terms of N plane waves and the slab is divided into I finite difference planes, the number of elements in f becomes 2(2I + 1)N (the electric field is expanded in one more plane than the magnetic, as shown in Fig. 1).In 3D problems this leads to huge equations but using iterative solvers together with efficient preconditioners [26], M never has to be constructed explicitly and the system can be solved very efficiently.Our first approach was to create a transmission matrix by eliminating all the field variables inside the slab using Eq. ( 5) recursively but we found this procedure to be numerically unstable and inefficient due to the involved repeated matrix multiplications.The instability is related to the presence of the complex valued propagation constants in the zdirection [19], though not explicitly considered in the finite difference formalism.In the current method, this instability is completely eliminated as there are no matrix multiplications and the exponential growth and decay never accumulate.

Reflection of beams from periodic structures
Due to the linearity of Maxwell's equations, the fields excited by multiple sources can be computed individually and then added to give the complete field.We use this principle for superposing the reflection coefficients of individual plane waves to yield the reflection coefficient of the complete beam.This way the reflection of a finite sized beam can be divided into smaller problems, each of which can be solved by taking advantage of the crystal periodicity.A similar approach for electrostatic problems has been taken by Manzuri-Shalmani et al. [30].The achieved savings depend on the computational method used for solving the individual plane waves but most methods scale significantly worse than linearly with the number of unknowns and the savings can therefore be considerable.Dividing the big problem into several smaller ones also saves computer memory, which is often of concern in 3D computation.An additional advantage is that once the reflected fields of the individual plane waves have been computed, they can be easily combined in different ways to yield the reflection of different incident beams.
We limit ourselves to considering the reflection only as the transmission can be computed from a very similar formula, or in the case of transparent materials, from the conservation of energy.We will also assume that there is no absorption or gain for z > h, i.e. in the half space of the incident and reflected beams.This is not essential to the method, but for absorptive materials, the incident and reflected intensities are z-dependent which complicates the interpretation of the results.

The incident and reflected fields
Consider an optical beam at a given frequency ω, propagating from z = +∞ towards the photonic crystal surface at z = h.A general expression of such a beam is given by a Fourier integral of the eigenvectors given in Eq. ( 9) but in the realm of numerical computation we will sample the wave vector and use a series representation where w m is the propagation constant in z-direction (the eigenvalue), corresponding to the transversal wave vector K m , and Ψ 1− m and Ψ 2− m are the corresponding eigenvectors.The shape and polarization of the beam are defined by the coefficients c 1 m and c 2 m .The difference to the homogeneous-medium expansion as given in Eq. ( 7) is that the transversal wave vectors K m are not required to be multiples of the photonic crystal reciprocal vectors and the incident beam is not necessarily related to the lattice symmetry.
Our aim is to treat each K m as a Bloch vector of the photonic crystal lattice and compute the reflection of each plane wave exp ( jK m • r − jw m z) Ψ l− m , l = 1, 2, individually from Eq. ( 12).The results are then added to give the reflection of the complete beam.Since each of the incident plane waves is a Bloch wave, it suffices to consider only one unit cell of the photonic crystal.Notice that K m need not be limited to the first BZ of the photonic crystal since we can always write it in the form K m = K m + G n , where K m is a vector in the first BZ and G n is a suitable reciprocal vector.
Once the reflection of each plane wave in Eq. ( 13) is computed, we can write the total reflected field as where Θ 1 m and Θ 2 m , defined in Eq. ( 15), are the lattice periodic reflection functions for individual eigenvectors Ψ 1− m and Ψ 2− m , respectively, Here, w m,n is the eigenvalue corresponding to the transversal wave vector (K m + G n ) and a l+ m,n and b l+ m,n are scalar coefficients obtained from the solution of Eq. ( 12).Qualitatively speaking, summing over m means summing over different incident angles and summing over n means summing over different Bragg orders.The wave vectors of the incident and reflected beams are illustrated in Fig. 2.

Reflection coefficient of the beam
The reflected power is computed by integrating the time averaged Poynting vector P = (1/2)ℜ (E × H * ) over a surface S. Taking the differential surface vector dS to be parallel to z, the projection dS • P can be written in terms of the transversal field components and the reflected power becomes (16) where the operators E and H have the same meaning as in Eq. (10).The reflectance is obtained as a quotient of P and the power in the incident beam.
The wave vectors of the incident and reflected fields.a) A single plane wave with a wave vector (K 0 − w 0 u z ) is incident from the homogeneous medium (HM) to the surface of the photonic crystal (PC) and b) its reflection is expressed in terms of the different Bragg orders, both propagating and evanescent (not shown).c) An incident beam is decomposed to plane waves with wave vectors [K m − w m u z ] and d) the reflection is expressed in terms of the corresponding Bragg orders.Equation ( 16) can be greatly simplified if we place a few restrictions on the transversal wave vectors K m .The aim is to select the vector set [K m ] such that the harmonic basis functions defined by vectors [K m + G n ] form an orthogonal and unique (that is, no two functions are the same) set over some surface S. The orthogonality can be achieved if the vectors K m are selected among the reciprocal vectors of a supercell defined by L 1 a 1 and L 2 a 2 , where L 1 and L 2 are integers and a 1 and a 2 are the lattice vectors of the photonic crystal.The requirement of uniqueness is satisfied if we include only vectors which are inside the first BZ of the photonic crystal.Then [K m ] together with the reciprocal vectors of the photonic crystal [G n ] form the set of reciprocal vectors for the supercell.These restrictions limit the beams which can be expressed using Eq. ( 13).The introduction of the supercell makes the incident beam (and the reflected field) periodic, which is usually not a problem since adjacent beams can be decoupled by giving L 1 and L 2 sufficiently large values.A more fundamental restriction is the upper bound on the length of K m , which defines a cone or a numerical aperture (NA) limiting the possible plane waves in the incident beam.This in turn sets a lower bound on the achievable localization on the crystal surface.Whether this is of practical concern or not, depends on the measurement setup.In our case, the NA of the optics used for the measurement is more restrictive than the NA originating from the limitations on K m .In any case, if these conditions are too restrictive, the reflected power can always be evaluated directly from Eq. ( 16), disregarding the simplified formula we are about to develop.
It should be noted that whenever the photonic crystal is made of completely transparent materials, the incident beam should not contain evanescent components, i.e. one should set Relaxing this condition opens up a possibility that the incident beam excites guided slab modes, which cannot be decoupled from each other simply by (C) 2005 OSA increasing the periodicity interval of the sources.These modes should in principle have an infinite amplitude because they are excited by an infinite number of sources and they do not decay, but in practice the amplitude is determined by numerical effects and is more or less random.Even if guided modes are excited, it is not impossible that the reflection coefficient of intensity is computed correctly, since guided modes, by definition, do not contribute to radiation.However, we do not have a proof on this and therefore we avoid computing the reflection in the presence of guided modes.
Having introduced the necessary restrictions, we now proceed with the simplification of Eq. ( 16).Performing the surface integral over the afore mentioned supercell, the integration of the harmonic basis functions evaluate to Kronecker delta functions with arguments (m, m ) and (n, n ) and the double summations cancel to single summations.The numerous vector products can be evaluated using the properties of the eigenvectors given in Eqs.(10) and (11).After a cumbersome but straightforward calculation we can show that the reflection coefficient of a beam defined by the coefficients c 1 m and c 2 m is given by where the parameters R 1 m , R 12 m and R 2 m do not depend on the coefficients c 1 m and c 2 m and can therefore be calculated a priori as Here δ (w m,n ) = 1, if the eigenvalue w m,n is real valued and 0 otherwise.The physical origin for the appearance of δ (w m,n ) is that the eigenvectors with an imaginary wave number are evanescent and do not carry energy in a direction parallel to the z-axis.

PMMA opals for visible light
The crystals under investigation were self-assembled [31] in the etched areas of double side polished 100 silicon substrates, from a 4.5 wt-% suspension of PMMA spheres in de-ionized water, using the vertical deposition technique [32,33].The PMMA spheres were fabricated with a median diameter of a = 268 nm, using the modified surfactant free emulsion polymerization technique described elsewhere [34].The main advantage of using PMMA as medium is that its material properties are well known and that it is commonly used also as an optical medium, making the optical design and subsequent theoretical evaluation process straightforward.Other advantages include its suitability for patterning with electron beam lithography [35].
Prior to the opal growth, the silicon substrates were cleaned for three hours in a 1:1 solution of sulfuric acid (95%) and hydrogen peroxide (30%).The substrates were then hydrophilized during three hours in a 1:1:5 agent of hydrogen peroxide (30%), ammonium hydroxide (25%), and de-ionized water, and finally blown dry with nitrogen.The opal samples were grown at room temperature and normal atmospheric pressure, at a vertical drawing speed of v 0 = 2.6 mm/h, resulting in films of approximately 18 monolayers, or a thickness of 4.2 µm.After the growth, the samples were sintered at 80 • C during 90 minutes.
The optical reflectance spectra were measured using a NanoSpec III spectrophotometer equipped with a microscope.The incident white light was supplied through a microscope objective such that the focus point was in front of the sample leading to an illumination of a relatively large area on the surface.The reflected spectrum was measured through the same objective but the light was collected from a small spot and focused on a small aperture of a spectral analyzer.Two magnifications were used in the measurements: 50x, which had a spot size of 3 µm on the sample, numerical aperture NA=0.55 and acceptance angle 33.4 • and 10x, with the corresponding values being 15 µm, 0.3 and 17.5 • .

Comparison of simulations and measurements
Before presenting the results, we will make a few simplifying observations: 1) Since the incident light was unpolarized, we computed the reflectance spectrum for e-and h-polarizations (e-polarization has electric field parallel to the sample surface and correspondingly for h-polarization) separately and averaged.For a given polarization, the ratio of c 1 m and c 2 m is fixed and Eq. ( 17) can be written in terms of a single variable c m .
2) The direction of a plane wave is determined by its wave vector components K and w such that the angle between the z-axis and the plane wave is given by φ = tan −1 (|K|/w) and the angle between x-axis and the plane wave by θ = tan −1 (K y /K x ), where K x,y = u x,y • K.It turns out that for a sufficiently small φ , the reflectivity varies only little with θ .Particularly, at the maximum acceptance angle of the collecting optics, φ = 33.4• , the variation of reflectivity with θ is at most 1.5% for any wavelength considered in the measurement.Therefore, we can take advantage of the rotational symmetry of the incident beam and integrate Eq. ( 17) over θ with a little loss in accuracy.Then it is sufficient to sample K only over a line instead of a surface, thus reducing the computational burden significantly.The low dependence on θ can be addressed to the high, six fold rotational symmetry of the FCC-lattice about [111]-direction (z-axis) and low refractive index contrast between PMMA and air.However, it should be noted that if ω is high enough to allow more than one diffraction order, then the orientation of K may indeed have a significant effect on the reflectivity.In our samples, the second diffraction order arises at a vacuum wavelength of 464 nm, which is outside the spectral scope of the measurements.
Incorporating the afore mentioned simplifications into Eq.( 17) gives where ∆ K = |K m+1 − K m | is the uniformly spaced difference in the length between two consecutive K-vectors on the sampling line and R m is the appropriate linear combination of R 1 m , R 12 m and R 2 m , defined in Eqs.(18) and (19).Simulations were performed on a grid of 16 by 16 plane waves and 16 finite difference planes per sphere diameter.With 18 monolayers in the sample, this results in an equation with 243200 unknowns.Solving the system once, requires typically 25 to 35 matrix iterations and takes about 10 to 15 seconds on a standard 2 GHz PC (This claim applies to this particular geometry, the convergence of the iterative solver depends strongly on the spatial characteristics of the dielectric function).The spectrum was computed at 75 ω-points and K was sampled at 40 points for each value of ω.We assumed PMMA to be transparent and to possess a frequency independent refractive index n PMMA = 1.489.The dispersive and absorptive dielectric constant of the silicon substrate was taken according to Ref. [36].The shape of the incident beam on the sample was not precisely known but the optics used for collecting light had a smaller NA than the optics used for illumination.Therefore we set c m = 1 for all plane waves that are accepted by the collecting optics, i.e. sin(φ m ) <NA, and c m = 0 otherwise.
Figure 3 shows the measured reflection spectra from a typical sample and the corresponding simulations for the two different spot sizes.For the smaller spot size (NA=0.55)the reflection curve has gone through a considerable blue shift, the main maximum has broadened and the fringes have leveled out.These effects are well predicted by the simulation, even up to some very small details on the short wavelength side of the reflection peak.According to the simu- lation, the highest reflectivity should occur for the larger spot size (NA=0.3)but the opposite is observed in the measurements.This is most likely due to defects in the crystal falling under the larger beam.The good match between the simulation and the measurement verifies that the crystal is of good quality and the approximations made in the simulation are valid.

Polystyrene opals for infrared light
The second photonic crystal structure studied in this work was fabricated by self organization of polystyrene nano-particles using the vertical solvent evaporation technique [8], where a substrate is placed vertically in a container partially filled with a nano-sphere dispersion.The nano-particles form a lattice on the substrate as the solvent evaporates.A small piece (5 mm × 20 mm) of a polished n-type gallium arsenide (GaAs) wafer was used as a substrate and the dispersion consisted of 5 wt-% of polystyrene nano-spheres in water with an average sphere diameter of 460 nm.The solvent evaporation temperature was 55 • C.
The structure of the photonic crystal lattice was characterized by scanning electron microscopy (SEM).A typical image from the crystal surface is shown in figure 4. Angle resolved reflectance spectra were measured in a goniometer setup.Light from a halogen lamp was passed through a monochromator and slightly focused on the sample.The reflected intensity was measured using standard lock-in techniques through a small aperture resulting in an angular resolution of about 3

Results
We measured the reflectance spectra for several values of the incident angle φ using e-polarized light.Simulations were performed with the same resolution as in the previous section but the incident beam was so wide (about 1 mm in diameter) and slightly focused that we assumed a plane wave excitation.A frequency dependent and complex valued dielectric function was used for both GaAs [37] and polystyrene [38].Based on the measured fringe positions, as shown in figure 5, we concluded that there are 47 monolayers in the sample.In the simulation, we assumed the Poynting vector's projection to lay along the Γ − M direction of the triangular lattice on the crystal surface.Simulations predict the positions of the reflection peak and fringes well but the difference in amplitudes is quite pronounced.Especially, the linearly increasing trend in the reflection is not captured by the computations at all.On the contrary, the simulated reflection has a decreasing trend, which can be shown to be caused by the infrared absorption of polystyrene.We assume the increasing trend to originate from a wavelength dependent scattering but the exact nature of this process is not known to us.However, the trend is present in every sample we have prepared (using the same batch of polystyrene spheres), even if silicon was used as a substrate material.The details of the scattering process are a topic of our further research.
Even though the microspheres are expected to crystallize into an FCC-lattice, they may also crystallize into a hexagonal close packed (HCP) lattice or a random mixture of the two [11].In an attempt to determine the crystal structure, we simulated the reflectance from a mixed lattice, as shown in Fig. 5(d).Unfortunately, the difference between FCC and mixed lattices is so small that the properties of the crystal cannot be established in this way.

Conclusion
We presented an efficient algorithm for computing electromagnetic fields in doubly periodic, finite thickness structures and applied it for computing reflection coefficients from opal photonic crystals.It was found that together with efficient preconditioners, our method is capable of solving large scale three dimensional problems on a standard PC very quickly.We also used a scheme for decomposing an arbitrary beam into plane waves and showed how to compute the reflection coefficient of the beam by superimposing the reflection coefficients of individual plane waves.The computation time in this scheme scales as O(NT ), where N is the number of plane waves needed to expand the incident beam and T is the time required for solving the Bloch periodic problem of a single plane wave.We demonstrated the method by computing the reflection coefficient of a focused beam from PMMA opals on silicon substrates.The match between the measurement and simulation was found to be remarkably good, thus showing the power of our method and the high quality of the sample.
We also studied polystyrene opals on GaAs substrates, which were found to be of lower quality.However, the reflection was measured only using a wide beam, which may explain the low reflectance.Simulations were able to predict the positions of the reflection peak and fringes well but the amplitude was not modeled correctly.We also found a wavelength dependent scattering in polystyrene opals, which needs to be thoroughly investigated.

Figure 3 .
Figure 3. Measured and simulated reflection spectra from a typical PMMA on silicon opal.Numerical aperture of 0.55 corresponds to 50x magnifying optics with a spot size of about 3 µm and numerical aperture 0.3 corresponds to 10x magnification and a spot size of 15 µm.

Figure 4 .
Figure 4.A typical SEM image of the polystyrene on GaAs opal.

Figure 5 .
Figure 5. Reflection from a polystyrene on GaAs opal for different incident angles: a) φ = 20 • , b) φ = 30 • and c) φ = 40 • .d) The incidence angle and the measurement data is as in c) but the simulated reflection is computed using a random mixture of FCC and HCP lattices.Curves marked with stars are simulated and continuous lines are measured.