Abstract
Accurate simulation from digital, submicron, optical elements is obtained by finite difference time domain (FDTD) results that are phase analyzed as sources for Huygens wavelets on fine scales much shorter than the wavelength used. Results, from the MIT electromagnetic evaluation program, are renormalized by a method here called “refractive impulse.” This is valid for polarized responses from digital diffractive and focusing optics. The method is employed with plane wave incidence at any angle or with diverging or converging beams. It is more systematic, more versatile, and more accurate than commercial substitutes.
1. Introduction
Modern optical components train submicron wavelengths through yet smaller optical diffractive and imaging elements [1]. Meanwhile, digitized optics enables novel applications with asymmetric diffractive and focusing components and with holographic systems. Classical wave theory can indicate only coarse descriptions of the response, though several techniques have been described that address the problem. For example, attempts at coupled-wave approaches [2–4] are typically limited in configuration and approximate in application. Generally, the solution requires integrations or summations in three dimensions, including the profile of an optical component.
Commercial software, typically, attempts simulation of the optics of microelements by adapting dynamical diffraction [2–4]. This adaptation involves treating a two-dimensional grating as if it were a three-dimensional crystal. However the method is inappropriate for many reasons and the predictions for diffracted intensities cannot be better than 50% accurate. The law of diffraction of light, having wavelength , from a plane grating, of spacing , is given by [5] where is the order and is the scattering angle. Diffraction from three-dimensional crystals is quite different. Bragg’s law gives which approximates to when as in transmission electron microscopy. However, the interplanar spacing in Bragg diffraction becomes multivalued and is peculiar for each diffracted beam following the method of indexation, , , [4]. Moreover, the scattering angle is approximately twice the Bragg angle. Unlike diffraction from gratings, Bragg diffraction requires the angle of incidence to change when high orders are measured. There are yet further fundamental differences that are overlooked in the commercial software. Bragg diffraction is not refracted since there is only one zero order beam and it lies in the direction of incidence; all Bragg diffraction is specular, while each diffracted beam has its own peculiar interplanar spacing, so that analysis that uses a single spacing for all beams is inaccurate, whereas, in crystals, the Bragg condition is satisfied by the Ewald sphere construction and this condition does not apply in diffraction from two-dimensional (2D) gratings. The Ewald sphere applies weakly in thin specimens in transmission electron microscopy of crystals by application of the deviation parameter [6]. Conjoining the two methods described by (1) and (2) is misleading, because diffraction from a 2D grating is different from 3D crystal diffraction.
For these and further reasons, we revert to grating optics in order to simulate diffraction and focusing that are due to microoptics with profiled elements. The response to microoptical elements, where , is represented by finite difference time domain (FDTD) computation. Consistent calculations on a very fine grid are taken to average individual amplitudes from correspondingly small cells arranged on a plane behind the optical microsystem. Such a plane is represented in Figure 1 by the dashed line. The impulsive amplitudes at voxels are added together in the near or far field after accounting the phase changes due to ray paths. This method has been justified by studies on proximity effects about multiple elements. It provides both accuracy in computed response and calculated differentials in polarization dependence. When optical components are large, computational economy is gained by combining stitching [7, 8] of local finite difference time domain (FDTD) results to Huygens wavelet analysis on a fine scale. The method is versatile and able to describe parallel wavefronts as well as divergent ones, before and after digital focusing and diffractive beam splitting. Though FDTD is generally used for understanding relative changes, renormalized results are stabilized by a method that we call “refractive impulse.” Detailed obliquity effects are included in the calculation. The method also provides internal consistency, since scattering angles and diffractive interference are calculated internally and are not given as parameters for the calculations. Often, the calculations can be speeded up by perceptive use of symmetries. The method has been applied to systems containing transmissive optical elements.
The method can be used in either transmission or reflection. In either case it provides a versatile solution not only in classical optics with profiled grating structures [9–11] or in combined digital focusing and diffractive elements or in holographic patterning [12] but also in other diverse applications.
There are further practical reasons for adopting this combination of FDTD phase shifts on a plane with wavelet propagation. The combination is equally adapted to divergent, convergent, or parallel beams, whether for analysis on symmetric or asymmetric gratings, analysis on digitized Fresnel lenses or conventional optics, or polarized or nonpolarized beams, and so forth. “Refractive impulse” refers to the method used for interfacing the FDTD and digitized Fresnel microoptics, including normalization of consistently iterated beams.
2. Finite Difference Time Domain Simulations
On a scale of 1/40 μm cubic “voxel” elements, Maxwell’s equations were solved using the MIT electromagnetic evaluation program (MEEP) [13]. Radiation from a laser source, ~800 nm wavelength, was represented as a collimated plane wave or as a divergent wave from a point source [14]. Responses from submicron optical elements were simulated downstream for various polarization conditions. One example is illustrated in Figure 1. Notice the change in wavelength due to the refractive index of the glass to the right side of the optical elements. A plane is selected with mean amplitude zero and the amplitudes are normalized, on nearby planes, before conversion to phase shifts : This is valid where the shift . This formula depends on the fact that is modulo in , since this condition has the result that . The phase is used to advance or retard the Huygens wavelets to be described in Section 3.
Calculations are done in 2 or 3 dimensions (2D or 3D) as required. Greater resolution and speed of calculation are obtained at the lower dimension owing to reduced computing power required. Snapshots are recorded over series of moments. The example shows normal incidence, but a wide range of configurations have been simulated, including collimation and focusing by Fresnel type lenses.
Normal methods of finite element analysis (FEA) are applied. Instead of a strict boundary condition, the boundary layer is applied with exponentially increasing absorption towards the edge. Optimization is obtained after multiple trials. Transients are avoided by self-consistent iteration. Boundary absorption causes the amplitude of the wave to decrease in the direction of propagation. The typical consequence is that the simulations are relative and require renormalization, described in Section 4.
3. Huygens Wavelet Interference in Far Field
On a far field plane illustrated in Figure 1, each voxel near the microelement is the origin of a Huygens wavelet with phase change applied as just described. Illumination is plane polarized. Though the figure shows uniform illumination, variations in intensity over the field, due, for example, to a weakly divergent laser beam, are easily applied. Typically the wavelet phases are summed at a distant plane (as in Figure 2) after accounting for the wave oscillations over the path towards the plane. The 1D section shown is suitable for incident plane waves and regular gratings, but when the intensity varies or when focusing elements are used, 2D intensity plots are calculated.
Several details must be computed. Consider obliquity requirements, first at the optical elements and secondly at the far field plane. In conventional optics for an absorptive grating, the obliquity factor is “a function proportional to the amplitudes of secondary waves propagating in various directions according to Huygens’ principle. It is where is the angle between the normal to the original wave-front and the normal to the secondary wave-front” [3]. Since the FDTD wavefront has a varying normal, the local normal at each voxel is calculated and used in the obliquity calculation. We find that typical intensity outcomes due to the correct obliquity vary by only about ~1% from the approximate value due to applying the normal of the incident beam. The small difference is a consequence of cancellations of positive and negative components on either side of the normal. In some calculations, Huygens wavelet obliquity can be neglected by approximation. A second obliquity at the plane recorded at far field is represented by a cosine term when intensity is spread over longer space. This term can be completely neglected when refracted intensity is subsequently focused, for example, by conventional optics. Other features of FDTD calculations require more systematic treatment in accurate simulations including reflection and normalization.
4. Normalization for Reflected Rays
The simulated optical response of the system, whether represented in 1D as in Figure 2 or in alternative 2D, depends firstly on stitching phases derived from FDTD wave amplitudes to wavelet propagation at the near field and secondly on summing the amplitudes at far field. Prima facie, the FDTD amplitudes, provide valuable comparative information, as for polarization dependence, but limited accuracy in absolute calculation. We correct for the latter limitation by means of the method here called “refractive impulse” (RI) as a corrected average for the two types of simulation. The correction is applied to the wavelet response described in Section 3.
MEEP is iterative. A principal benefit is that, as a wavefront advances, reflected beams are incorporated into self-consistent conditions at elemental boundaries. Meanwhile, graduated absorbing layers at global boundaries serve to minimize unwanted reflections. The absorption results in a loss of current in the simulation. Without correction, resulting simulations will be relative. However, it is possible to derive accurate, quantitative results by applying corrections for reflected rays. Consider an arrangement illustrated in Figure 1, where a beam is incident on a glass block containing digital microelements. Suppose the simulation is required for a feedback loop to be used in stabilization of an optical light source and coupling optics. Accurate calculations require absolute comparison of the transmitted beam with the incident beam.
The current inside the block is reduced by two principal effects: firstly due to reflected wavefronts from the face of the glass block and secondly due to absorption at the global boundaries. Here, the theoretical graduation in refractive index at the boundaries is designed to minimize computed reflections. Typically, the matching of the incident beam at the global boundary is achieved by logarithmically graded refractive indices. The effect of absorption at the global boundaries can be computed by switching the refractive index to and summing intensities within boxes that are set up to measure the glass block reflection and transmission (Figure 3). This computes the absorption, . The widths of the blocks (dashed lines) are an integral number of wavelengths, in air or in glass of refractive index . From the four intensity summations, the corrected ratio of transmitted to incident beam intensities is calculated and used in the stabilization of instrumentation through refractive beam splitting.
5. Conclusion
The method described for simulating responses due to submicron optical components is versatile. Its accuracy is increased systematically by normalization. FDTD methods are then applied to simulate accurate simulations of response. For typical 1D solutions, the results can be obtained in a few minutes plus some editing in selection and formatting of optimum amplitude profiles. For 2D calculations, involving polarized diverging beams and focusing optics, the calculations are longer. However the time can often be reduced by perceptive use of symmetries and by approximations made to employ individual stitches multiple times (e.g., by joining horizontal stitches locally to cylindrical fields).
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.