Generation of an incident focused light pulse in FDTD

A straightforward procedure is described for accurately creating an incident focused light pulse in the 3-D finite-difference time-domain (FDTD) electromagnetic simulation of the image space of an aplanatic converging lens. In this procedure, the focused light pulse is approximated by a finite sum of plane waves, and each plane wave is introduced into the FDTD simulation grid using the total-field/scattered-field (TF/SF) approach. The accuracy of our results is demonstrated by comparison with exact theoretical formulas. © 2008 Optical Society of America OCIS codes: (050.1755) Computational electromagnetic methods; (260.1960) Diffraction theory; (080.3630) Lenses. References and links 1. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House, Boston, 2005). 2. H. K. Roy, Y. Liu, H. Subramanian, D. Kunte, P. Pradhan, R. K. Wali, J. Koetsier, M. J. Goldberg, Z. Bogojevic, and V. Backman, “Detection of the colorectal cancer (CRC) field effect through partial wave spectroscopic microscopy (PWS),” Gastroenterology 132, A169 (2007). 3. E. Wolf, “Electromagnetic diffraction in optical systems. I. An integral representation of the image field,” Proc. Roy. Soc. A 253, 349–357 (1959). 4. G. S. Smith, An Introduction to Classical Electromagnetic Radiation (Cambridge University Press, New York, NY, 1997). 5. R. W. Ziolkowski and J. B. Judkins, “Full-wave vector Maxwell equation modeling of the self-focusing of ultrashort optical pulses in a nonlinear Kerr medium exhibiting a finite response-time,” J. Opt. Soc. Am. B 10, 186–198 (1993). 6. D. B. Davidson and R. W. Ziolkowski, “Body-of-revolution finite-difference time-domain modeling of spacetime focusing by a 3-dimensional lens,” J. Opt. Soc. Am. A 11, 1471–1490 (1994). 7. B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems. II. Structure of the image field in an aplanatic system,” Proc. Roy. Soc. A 253, 358–379 (1959). 8. W. Sun, S. Pan, and Y. Jiang, “Computation of the optical trapping force on small particles illuminated with a focused light beam using a FDTD method,” J. Mod. Opt. pp. 2691–2700 (2006). 9. L. Cheng, H. Zhiwei, L. Fake, Z. Wei, D. W. Hutmacher, and C. Sheppard, “Near-field effects on coherent antiStokes Raman scattering microscopy imaging,” Opt. Express 15, 4118–4131 (2007). 10. W. A. Challener, I. K. Sendur, and C. Peng, “Scattered field formulation of finite difference time domain for a focused light beam in dense media with lossy materials,” Opt. Express 11, 3160–3170 (2003). 11. P. Török, P. R. T. Munro, and E. E.Kriezis, “High numerical aperture vectorial imaging in coherent optical microscopes,” Opt. Express 16, 507–523 (2008). 12. C. Liu and S.-H. Park, “Numerical analysis of an annular-aperture solid immersion lens,” Opt. Lett. 29, 1742– 1744 (2004). 13. J. Liu, B. Xu, and T. C. Chong, “Three-dimensional finite-difference time-domain analysis of optical disk storage system,” Jpn. J. Appl. Phys. Part 1 39, 687–692 (2000). (C) 2008 OSA 10 November 2008 / Vol. 16, No. 23 / OPTICS EXPRESS 19208 #101373 $15.00 USD Received 9 Sep 2008; revised 15 Oct 2008; accepted 15 Oct 2008; published 5 Nov 2008 14. K. Şendur, W. Challener, and C. Peng, “Ridge waveguide as a near field aperture for high density data storage,” J. Appl. Phys. 96, 2743–2752 (2004). 15. S.-Y. Sung and Y.-G. Lee, “Trapping of a micro-bubble by non-paraxial Gaussian beam: Computation using the FDTD method,” Opt. Express 16, 3463–3473 (2008). 16. K. Choi, J. W. M. Chon, M. Gu, and B. Lee, “Characterization of a subwavelength-scale 3D void structure using the FDTD-based confocal laser scanning microscopic image mapping technique,” Opt. Express 15, 10,767– 10,781 (2007). 17. M. Wang, J. Wu, J. Xu, D. Ge, H. Li, and J. Feng, “FDTD simulation on the interaction between Gaussian beam and biaxial anisotropic metamaterial slabs,” Int. J. Infrared Millim. Waves 29, 167–178 (2008). 18. Y.-F. Chau and D. P. Tsai, “Near-field optics imaging in silica waveguide using near-field scanning optical microscope,” Jpn. J. Appl. Phys. Part 1 46, 238–242 (2007). 19. J. B. Judkins and R. W. Ziolkowski, “Finite-difference time-domain modeling of nonperfectly conducting metallic thin-film gratings,” J. Opt. Soc. Am. A 12, 1974 (1995). 20. J. B. Judkins, C. W. Haggans, and R. W. Ziolkowski, “Two-dimensional finite-difference time-domain simulation for rewritable optical disk surface structure design,” Appl. Opt. 35, 2477 (1996). 21. M. Born and E. Wolf, Principles of Optics : Electromagnetic Theory of Propagation, Interference and Diffraction of Light, 7th ed. (Cambridge University Press, Cambridge, 1999). 22. W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes in C: The Art of Scientific Computing, 2nd ed. (Cambridge University Press, Cambridge, 1992). 23. S. Bochkanov and V. Bystritsky, “Computation of Gauss-Legendre quadrature rule nodes and weights,” Alglib.net – Web Resource. Date Accessed: 08/2008. URL http://www.alglib.net/integral/gq/glegendre.php.


Introduction
The finite-difference time-domain method (or simply FDTD) is a numerical electromagnetic simulation technique that is commonly used in a wide variety of applications.The main advantage of FDTD is its ability to solve Maxwell's equations directly in the time domain, thereby providing accurate broadband solutions to electromagnetic wave propagation and scattering problems.In this method, the electric and magnetic fields are placed in discrete positions in space, and conventionally arranged in a 3-D rectangular mesh called the FDTD grid.The field components are also evaluated at discrete instants in time.The staggered positioning of the electric and magnetic fields in both space and time provides second-order accuracy, and allows for an explicit leap-frogging scheme for obtaining the fields directly in time domain [1].
There are several situations where the analysis of a focused light pulse in the time domain might convey more direct and useful information than that of a monochromatic focused light beam.A straightforward example is the modeling of an ultrafast optical system that generates a laser light pulse with duration comparable to the mean wavelength.This "ultrashort pulse" has a broadband spectrum; and when focused by a converging lens, constitutes a direct example of a focused light pulse.Another case in which a focused light pulse is of interest is the modeling of a white-light optical illumination system, such as a regular bright-field microscope.The broadband character of this system makes the focused light pulse an ideal tool for obtaining information in a wide spectral range.Some promising experimental results obtained using a spectroscopic microscopy technique for cancer research [2] has motivated our efforts toward the FDTD modeling of a broadband imaging system, which produced the results in this paper.
FDTD has been successfully applied to the electromagnetic analysis of many optical systems in which the effects of a nonzero wavelength are significant.The most important limiting factor in the FDTD simulation of optical systems is the need to discretize the spatial domain by a fraction of the smallest wavelength in the illumination or the smallest spatial feature in the geometry.This usually limits the size of the FDTD simulation space to a small portion of the optical system that is comparable in size to the wavelength of the light.In many cases, this portion covers the area over which the most important interactions between light and matter take place, such as the focal area of a lens.In this paper, we describe a technique that can be applied to such a problem.Specifically, we show how to create an incident focused light pulse in the FDTD simulation of the focal area of an aplanatic (aberration-free) converging lens.The geometry of this setting is illustrated in Fig. 1(a).
The method that we use to create the focused light pulse is known as the total-field/scatteredfield boundary (TF/SF) method.In this method, a special boundary is constructed in the FDTD grid, and field corrections (or consistency conditions) [1] are applied on the field components on the boundary to maintain a perpetual contrast between the fields on the interior (the totalfield region) and those on the exterior (the scattered-field region) of the boundary.The creation of an incident field using this method is often referred to as injection, because the incident field has the appearance of being injected into the total-field region as a result of the mentioned field corrections.If no scatterer is present within a theoretically perfect TF/SF boundary, the incident field is strictly confined within the boundary.If there are scatterers within the boundary, the field inside the boundary is the total field, and the field outside the boundary is the scattered field.
Although the TF/SF method is theoretically capable of injecting any incident field, it is mostly used to inject plane waves, owing to the simplicity of the required theory and implementation [1].The TF/SF plane-wave source is now a standard feature in any commercial FDTD software package.Therefore, it would be advantageous to express the electromagnetic field of a focused light pulse as a linear combination of plane waves, for the problem of injecting a focused light pulse then reduces to the well-studied problem of injecting a plane wave.Such a linear decomposition is not only theoretically possible, but unique as well; and is explicitly manifested in the plane-wave spectrum of the electromagnetic field, which is a 2-D integral representation of the field in terms of plane waves [3,4].We will present the plane-wave decomposition of the electric field in the image region of the lens in Section 3. The practical concern in the implementation of this decomposition in FDTD is to approximate this planewave spectrum by a finite collection of plane waves, each of which is to be introduced into the FDTD grid by a TF/SF plane-wave source.
The rest of the paper is organized as follows.In Section 2, we present a review of the literature on the problem.In Section 3, we describe the theoretical basis for our approach, and discuss its relationship with previous studies.In Section 4, we apply the theoretical results to the FDTD method, and elaborate on some implementation details.In Section 5, we present an example of the mentioned FDTD implementation, and provide some field snapshots for the visualization of the focusing mechanism.In Section 6, we demonstrate the validity and accuracy of our results by comparison to exact theoretical formulas in the literature, and present a rigorous error analysis.Finally, we summarize our work in Section 7.

History
There have been a number of studies on the FDTD simulation of the focal area of a converging lens; although a direct time-domain implementation of a focused pulse has not yet been introduced.Some indirect approaches that utilize Gaussian beams with time profiles have been developed [5,6].However, these approaches rely on the fact that the exponential beam profile of a Gaussian beam at the focal plane is known a priori, and is frequency-independent.In the focal area of a converging lens, however, the shape of the beam profile is usually frequencydependent, and not completely known in advance.It is determined by the properties of the incident and focused wavefronts, and some additional factors pertaining to the diffraction of the wavefront at the exit pupil.These have to be included as a whole in the solution for a focused light pulse created by a converging lens.
A comprehensive theoretical account of the behavior of the electromagnetic fields around the focus of an aplanatic converging lens has been given by Richards and Wolf [7].They assume that a linearly-polarized monochromatic plane wave is incident on a converging aplanatic lens, and the resulting spherical wavefront satisfies the Kirchhoff boundary conditions at the circular exit pupil.The electromagnetic field in the image space of the lens then follows from vector dif- fraction theory.Because their results are given for monochromatic excitation, they have formed the basis for the time-harmonic FDTD simulation of a number of optical focusing systems, where both the pure-scattered-field approach and a variation of the total-field/scattered-field (TF/SF) approach has been employed for introducing the incident monochromatic focused beam.In the pure-scattered-field approach [8][9][10], the vector diffraction integrals for the incident focused light beam are evaluated numerically at each point within a 3-D volume in the FDTD grid, and added to the fields at those points at each time step.This is the most straightforward way in which the results in [7] can be translated into FDTD modeling; however, this is also the most computationally intensive.Some authors have preferred a simplified variation of the TF/SF approach [11,12], which requires that the vector diffraction integrals be calculated and applied as field corrections only on a 2-D rectangular region over the focal area.Although less computationally intensive, this method lacks some key advantages offered by the full TF/SF method.First, the scattered fields can be observed only above the rectangular region, which limits the use of the conventional near-field-to-far-field transformer (NFFFT) [1] because the total-field region cannot be bounded by a closed surface.This might impose a severe limitation on the FDTD modeling of some light scattering problems in which the conventional NFFFT is a crucial element.Second, the planar rectangular region over which the corrections are applied cannot have an arbitrary lateral extent, because it needs to cover the waist of the incident focused beam at that position.The Gaussian focused beam has drawn considerable attention in the FDTD literature, owing to the simple analytical expressions for the time-harmonic fields around the focal region [4].Due to this simplicity, the electromagnetic fields belonging to the time-harmonic Gaussian beam can be easily computed and introduced into the FDTD grid using either the pure-scattered-field approach [13][14][15] or the TF/SF approach [16][17][18][19][20].However, the adoption of a time-harmonic Gaussian beam suffers from the same disadvantages as those summarized above for a pulsed Gaussian beam.
In this paper, we demonstrate a direct time-domain application of the results in [7] in the context of FDTD modeling.We neither employ the simplistic Gaussian beam assumption, nor incorporate any frequency-domain calculation followed by inverse Fourier transform.In this regard, our results constitute the first FDTD implementation of a full TF/SF source for a light pulse focused by a converging lens.

Theory
The geometry of the optical system and the simulation space is shown in Fig. 1(a).The origin of the rectangular coordinates (x, y, z) is centered at the focal point F of the converging lens, and the z axis is chosen so as to coincide with the optical axis of the system.The spherical coordinates r = (r, θ , φ ) are also defined with origin F as shown.(We use the overbar sign for vectors, and the hat sign for unit vectors throughout the paper.)The refractive indices of both the object and image spaces are assumed to be unity.An x-polarized pulsed plane wave whose electric field is denoted by E i (t) in time is incident in the −z direction on the aplanatic converging lens, which creates a perfectly spherical converging pulsed wavefront W with focus F. The spherical wavefront W is diffracted by the circular aperture A that represents the exit pupil of the lens system.The FDTD simulation space, denoted by S, comprises a region around the focus F, which is in the far field of the pupil plane.The TF/SF boundary T , placed inside S as shown in Fig. 1(a), injects the focused light pulse into the FDTD grid.This procedure will be explained in more detail in Section 4. The time-domain electric and magnetic fields around F are denoted by Ē(r,t) and H(r,t), respectively.
In The polarization of the electric field on an incident ray (AP) and the corresponding ray in the image space after refraction through the lens (PF).
vector êi = x.Assuming that the converging wavefront W is perfectly spherical and geometrical optics is valid [3,21], the electric field along the corresponding ray PF on the wavefront W can be written as in which (θ , φ ) are the angles of incidence of the ray PF, r = (r , θ , φ ) denotes the coordinates on PF, the unit vector êw (θ , φ ) determines the polarization of the electric field on PF, and a(θ , φ ,t) is the strength factor of the same ray [7].In the realm of geometrical optics, the polarization of the electric field êw (θ , φ ) and the time waveform of the strength factor a(θ , φ ,t) do not vary along a ray, except a time advance by r /c in the latter.The attenuation of the fields by the distance r on the ray is a result of the intensity law in geometrical optics.
The relation between the waveforms E i (t) and a(θ , φ ,t) is determined by the effect of the converging lens on the incident plane wave.It can be shown from the intensity law of geometrical optics that for a thin aplanatic lens [7], in which f is the focal length of the lens.Furthermore, if the incidence angles at each surface in the lens system are small, each ray stays on the same side of the meridional plane (the plane containing the ray and the optical axis) during its passage from the lens system, and the angle that the electric field vector on a ray makes with the meridional plane is unchanged by refraction [7].This is a very important result from the perspective of our derivation, since it furnishes a very convenient way to evaluate the electric field polarization vector ê w (θ , φ ).If we define the polarization angle ψ(θ , φ ) of a ray as the angle that the electric vector on this ray makes with the meridional plane as shown in Fig. 1(b), it follows from the mentioned invariance of polarization angles that ψ(θ , φ ) is given by the simple expression tant and convenient result that greatly facilitates the FDTD implementation, as will be seen in Section 4.
The diffraction of the spherical wavefront W defined by ( 1) with ( 2)-( 3) at the exit pupil A is simplified considerably by the approximate Kirchhoff boundary condition [3], which asserts that the total electric field at each point on the plane of A is equal to the unperturbed electric field of the spherical wavefront W .The plane-wave spectral decomposition of the electric field belonging to W on the exit pupil A, along with the Kirchhoff boundary condition on A, yields the following solution for the electric field in the image region of the aplanatic lens [3,7]: in which ( 2) is used for the second equation.Here, the dot sign above a quantity denotes differentiation in time, θ max is the half angular range of the aperture A [see Fig. 1(a)], t = t + r • r/c is the advanced (or retarded) time, and r is the unit vector in the direction of r .The expression in ( 4) is merely the time-domain version of the equations (2.2) and (2.13) in [7], with the frequency-domain factor (−ik) translated into the time-domain derivative operator (1/c) d/dt.This is the reason for the appearance of the time derivative Ėi (t) in (4).A similar time-differentiation effect has been documented in other studies [5,6] on Gaussian beam pulses.
The key observation here is that the expression ( 4) is a decomposition of the electric field in the image region of the lens in terms of plane waves incident from a range of directions (θ , φ ).The plane wave from (θ , φ ) has polarization êw (θ , φ ), incidence vector kw (θ , φ ) = −r [see Fig. 1(b)], and infinitesimal amplitude waveform (referred to the origin r = 0) The exciting feature of this plane-wave decomposition from our point of view is the fact that it is given directly in the time domain.This property permits the application of the decomposition in (4) directly to FDTD analysis, as explained in the following section.

FDTD implementation
We had mentioned in Section 1 that, from the FDTD perspective, the most important advantage of expressing the electric field as a superposition of plane waves as in ( 4) is the simplicity and well-established nature of TF/SF plane-wave sources.However, we had also noted that some degree of approximation is inescapable in the FDTD implementation of such a superposition, because only a finite number of TF/SF plane-wave sources can be included in an FDTD simulation.Therefore, the practical problem at hand is to approximate the integral in (4) accurately by a finite collection of plane waves: in which 0 < θ n < θ max , 0 < φ m < 2π are the incidence angles of the plane waves, and t nm = t + r nm • r/c are the retarded times for these incidence angles.The amplitude factors α nm can be viewed as the numerical integration (or quadrature) factors for the double integral and other factors [22].Our double-quadrature problem is further complicated by the presence of the additional parameters r and t that determine the observation point and time, which modify the integrand of the double integral.If the waveform Ėi (t) is an oscillating function (as is the case for sine-modulated signals, for instance) the integrand also becomes oscillatory for large |r| (see the definition of t nm above), which introduces an additional difficulty.Because of these complexities, developing a specific strategy for approximating this double integral with controlled accuracy seems to be a challenging task.It might be a topic for future study to develop better quadrature schemes suited to this particular integral, but here we adopt a simpler, more heuristic approach to the problem.Specifically, we make use of the fact that the region of integration in ( 4) is rectangular, and perform one-dimensional quadratures along both θ and φ dimensions.In other words, we regard the θ integral in (4) for fixed φ as a one-dimensional integral (and vice versa), and utilize well-known one-dimensional quadrature rules for each of these integrals: where (θ n , a n ) and (φ m , b m ) are now regarded separately as one-dimensional quadrature positions and factors in their respective regions 0 < θ n < θ max and 0 < φ m < 2π.This method has the advantage of simplicity, but lacks global optimality in the two-dimensional (θ , φ ) region.
In Fig. 2, an example of this quadrature scheme is presented for N = 5, M = 7.The incidence angles θ n , φ m are determined by the specific quadrature rules used for each one-dimensional integral.In the following, we consider two prevalent one-dimensional quadrature schemes, namely, the extended midpoint and Gauss-Legendre quadrature rules [22].In the extended midpoint quadrature rule, the incidence angles θ n , φ m are equally-spaced in their respective intervals, and the amplitude factors a n , b m are simply the angular spacing widths: The Gauss-Legendre quadrature rule uses a more sophisticated algorithm for computing the incidence angles θ n , φ m and the amplitude factors a n , b m .The purpose of the algorithm is to (C) 2008 OSA render the integral exact when the integrand is a polynomial of order 2N − 1 or 2M − 1.The algorithm is usually carried out for the standard interval [−1, 1], and the results are then scaled to the desired interval.Computational routines for implementing the Gauss-Legendre quadrature rule are available both in print and online [22,23].Once the quadrature points and weights x j , w j for the standart interval [−1, 1] are found, they can be scaled to give the incidence angles θ n , φ m and the amplitude factors a n , b m for our original sum (7): In the example presented in the next section, we use the Gauss-Legendre quadrature rule for both sums in (7), and demonstrate its advantages over the extended midpoint rule by an error analysis in Section 6.
It is crucial to note that the injection of N × M plane waves into the FDTD grid is done in a single FDTD run.The linear character of the TF/SF correction procedure allows consecutive application of N × M TF/SF correction operations for N × M different plane waves in each main-grid time step.Therefore, the computational intensity involved in the TF/SF formulation is increased by a factor of N × M. The caveat is that this does not mean that the entire FDTD computational load goes up by the same amount.It is only the computational load resulting from the field-correction operations on the TF/SF boundary that increase by this factor.The percentage increase in the entire FDTD computational load due to this procedure gets smaller as the number of field components on the TF/SF boundary becomes less significant compared to the total number of the field components in the entire FDTD grid.Specifically, this happens when the average dimension d (in grid cells) of the FDTD grid is increased.Even if the average dimension of the TF/SF boundary is also kept proportional to d, the surface area of the boundary increases by d 2 , whereas the volume of the entire grid increases by d 3 .The resulting ratio 1/d between the number of field components on the boundary and in the entire grid still decreases as d is increased.In conclusion, it can be stated that the relative efficiency of the TF/SF focusedpulse source increases with increased grid size, for a fixed number of plane waves in (7).
It is worthwhile to mention in passing that the above linear combination operation is also valid in the presence of a nonlinear scatterer inside the TF/SF boundary, since it is not the scattered fields resulting from individual incident plane waves but the incident plane waves themselves that are being added in the above TF/SF formulation.The main grid only sees the focused pulse as the incident field, which is constructed as a linear combination of plane waves.The nonlinearity of the scattered field resulting from the incident focused pulse has no effect on the initial computation of the incident focused pulse, which is described here.
It should be emphasized that the time waveform Ėi (t) in ( 7) is referenced to the origin r = 0, namely, the focus F. However, TF/SF plane-wave sources are usually constructed such that the time waveform of the plane wave is referenced to one of the corners of the TF/SF boundary.Therefore, the plane waves should be advanced in time carefully to satisfy (5).In doing so, it is important to remember the velocity anisotropy property of the FDTD grid, especially for large grids where it might introduce a significant delay in the desired waveform.The effects of these technical details are discussed more quantitatively in Section 6.
The polarization of each plane wave in ( 7) can be easily determined using (3).The advantage of this vectorial view as opposed to the detailed consideration of all the field components (as in the formulas (2.26)-(2.27) in [7]) can be explained as follows.In many TF/SF plane-wave source implementations, the polarization of the electric field is specified by some kind of polarization angle that the electric vector makes with a reference plane that includes the incident ray.Since the z axis is the optical axis of our lens system, it would make sense to include this axis in the reference plane of the TF/SF plane-wave source.With this definition, the reference plane becomes identical to the meridional plane of the incident ray (see Section 3), and the definition of the polarization angle ψ(θ , φ ) in (3) can be used directly in the TF/SF plane-wave source implementation.

Example
In this section, we consider an example that demonstrates the functionality and accuracy of the proposed technique for injecting a focused light pulse into an FDTD grid.In our example geometry [see Fig. 1(a)], the focal length of the lens is f = 10 mm, and θ max = 23.58• , which corresponds to a numerical aperture (N.A.) of sin θ max = 0.4.The electric field of the plane wave incident on the lens is polarized in the x direction, and is a unit-amplitude sine-modulated Gaussian pulse in time: with τ = 3 fs, f 0 = 5.889 × 10 14 Hz, and t 0 = 13τ.The amplitude of the wavelength spectrum of this broadband pulse is within −40dB of its maximum between 400 nm and 700 nm, which corresponds to the visible range of the electromagnetic spectrum.As discussed in Section 1, this pulse can either represent the actual physical excitation as in the case of ultrashort pulses, or can be used as an auxiliary tool for analyzing the broadband behavior of the system as in the case of broadband illumination.The boundary of the FDTD grid used for our example simulation is represented by the solid rectangle S in Fig. 1(a).The TF/SF boundary T , shown in Fig. 1(a) as the dashed rectangle, is 1μm from the FDTD grid boundary S. The FDTD simulation parameters are grid size 10μm ×10μm ×12μm, grid spacing Δx = Δy = Δz = Δ = 20nm (resulting in a 500 × 500 × 600 grid), and time step Δt = (0.98/ √ 3)Δ/c.The computation domain is not truncated by any absorbing boundary, since the TF/SF boundary T contains the incident field to an accuracy greater than that of our example.The approximate electric field in (7) is constructed using 324 plane waves, with N = 9, M = 36.The plane wave component incident from direction (θ n , φ m ) has polarization êw (θ n , φ m ) and amplitude a n b m ( f /2πc) cos 1/2 (θ n ) sin(θ n ) Ėi (t) where Ėi (t) is the time derivative of the sine-modulated Gaussian pulse above.The incidence angles (θ n , φ m ) and the quadrature factors a n , b m are determined by the preferred quadrature rule for the numerical integral (7).An error analysis of different quadrature rules will be presented shortly; but we first assume that both the θ and φ sums in (7) are evaluated using the Gauss-Legendre quadrature rule (see Section 4.) Several field snapshots from the FDTD simulation taken between time instants t = 7τ and 17τ with 2τ intervals are shown in Fig. 3 (Movie 1), in which the x component of the electric field on the yz plane is plotted in dB on a gray scale, with black representing −40 dB and white representing 0 dB.The electric field is normalized by its maximum E max x = E x (0,t 0 ) at the focus, given by [see ( 12)-( 13) in the following section]

Error analysis
In this section, we present an error analysis of the example FDTD simulation in the previous section.We consider the x component of the electric field, for which exact theoretical formulas are given by Richards&Wolf [7] for time-harmonic dependence: where the tilde sign (∼) over the variables denotes time-harmonic quantities with time dependence exp(iωt), and the auxiliary variables I 0 and I 2 are given by I 0 (kr, θ ) = θ max 0 cos 1 2 θ sin θ (1 + cosθ ) J 0 (kr sin θ sin θ ) exp(ikr cos θ cos θ )dθ , (13) kr sin θ sin θ ) exp(ikr cos θ cos θ )dθ , (14) in which J 0 (•) and J 2 (•) are the Bessel functions of the first kind and order zero and two, respectively.
For the purposes of error analysis, we need to attain a high degree of accuracy in computing the θ integrals in ( 13)-( 14) numerically.This can be achieved by using a more sophisticated but time-consuming integration method, e.g., the adaptive Gauss-Kronrod quadrature rule [22].Since these integrals are for time-harmonic dependence, inverse temporal Fourier transformation into the time domain is also necessary for comparison with the FDTD results.It should be emphasized strongly that the aforementioned direct steps cannot be performed efficiently for practical use in FDTD, but they provide a highly-accurate reference for the error analysis of our approximate method.
The geometry of our error analysis is shown in Fig. 4(a).We consider the x component of Fig. 4. Normalized maximum dB-error for the numerical integral (7) describing the electric field belonging to the focused light pulse.The Gauss-Legendre quadrature rule is used for the numerical integral, and the error ( 15) is normalized globally by E max x .(a) The geometry of the error analysis.The elliptical focal spot [see Fig. 3(d)] is denoted by the gray ellipse.The error is calculated along line segments l x , l y , l z , and l xyz extending from the focus F in different directions.(b) Normalized maximum dB-error ε(r) (15) in the x component of the electric field on the line segments l x , l y , l z , and l xyz .the electric field in the FDTD simulation described in Section 5, recorded along line segments l x , l y , l z , and l xyz stemming from the focus F. If no x component exists at the exact position on the line segment, we take the average of the two nearest x components.
We define the normalized maximum dB-error ε(r) in the electric field waveform at a point r in space as follows: where E th x (r,t) and E FDTD x (r,t) are the time waveforms for the x component of the electric field obtained from theory ( 12)-( 14) and FDTD simulation, respectively, and E ref x is the field value by which the error is normalized.Depending on the value of E ref x , we define two kinds of error normalization.In global normalization, we define E ref x as the maximum field value at the focus F, namely, 11).In local normalization, we use the local maximum theoretical field value to normalize the error, namely, E ref x = max|E th x (r,t)|.In Fig. 4(b), we present the normalized maximum dB-error ε(r) on the line segments l x (marked •), l y (marked * ), l z (marked •), l xyz (marked ×) with respect to the distance d from the focus.For this data set, the Gauss-Legendre quadrature rule has been used in both the θ and φ sums, and the error has been normalized globally.It is seen that the normalized maximum error is below −35dB for the entire range.It is also clear from Fig. 4(b) that the normalized maximum error is larger on l z , which lies along the ẑ axis.We can explain this easily by inspecting Fig. 3: The energy of the focused pulse is concentrated mostly within an ellipsoidal region around the focus F with its major axis parallel to the ẑ direction.Because of the large amplitude of E th x (r,t) and E FDTD x (r,t) in (15), the error in this region is also comparatively larger.It is also of interest to investigate the impact of the quadrature rule and the error normalization on the error ε(r) in (15).One would expect a priori that the error in the numerical integration inherent velocity anisotropy of the FDTD grid [1], which is responsible for the dependence of the propagation velocity ṽ p (θ n , φ m ) on the direction of incidence and the excitation frequency.Velocity anisotropy becomes a more serious problem for larger FDTD grids, such as the one considered in Section 5.This can pose a serious problem for the FDTD implementation of the plane-wave decomposition in (7), since it may cause considerable differences between the expected waveform E d x (r,t) and the observed waveform E FDTD x (r,t) at a point.This problem was hinted at in Section 4, where velocity anisotropy was named as a potential source of error.The effects of velocity anisotropy can be lessened by using ṽ p (θ n , φ m ) instead of c while calculating propagation delays; see discussion regarding the TF/SF implementation in Section 4. For a broadband excitation, the center frequency in the excitation waveform can be used to calculate ṽp (θ n , φ m ).For the error analysis presented so far, this grid-velocity correction technique has been used.It is worthwhile to quantify the improvement that this correction offers in reducing the error.For this purpose, we revisit the case in Fig. 5(a) for Gauss-Legendre quadrature rule and global error normalization (4th dataset in legend), and investigate the effect of grid velocity correction on the error.In Fig. 5(b), the normalized error reproduced from Fig. 5(a) with grid-velocity correction is denoted by the * -line, and the normalized error recalculated without grid-velocity correction is denoted by the •-line.It is seen that grid-velocity correction offers an improvement of up to ∼15dB.This proves that grid-velocity correction is a very crucial component of the TF/SF-based implementation of the focused light pulse, especially for a large grid such as the one considered in Section 5 (500 × 500 × 600 grid with ∼20 cells per smallest wavelength).

Conclusion
In this paper, we introduced a new technique for creating an incident focused light pulse for the finite-difference time-domain (FDTD) numerical analysis of the image space of an aberrationfree lens.The theory of the technique was based on the time-domain version of earlier results by Richards&Wolf [7], and applied to the FDTD method using the total-field/scattered-field (TF/SF) formalism.This was achieved by interpreting the expression (2.26) in [7] for the image field as a plane-wave spectrum, and approximating this spectrum by a finite collection of plane waves.The error analysis in Section 6 confirmed our assertion that the time-domain method introduced in Sections 3-4 for creating a focused light pulse in FDTD offers a fairly accurate and efficient alternative to the brute-force calculation of the integrals in ( 12)-( 14) followed by inverse Fourier transformation into the time domain.
The mentioned error analysis revealed several important aspects of the numerical integral (7).In Fig. 4(b), it was shown that the Gauss-Legendre quadrature rule provides reasonable accuracy for a large grid [see Section 5], with 10 degrees per quadrature point in both θ and φ .A comparison of the Gauss-Legendre and extended midpoint quadrature rules was presented in Fig. 5(a), and the general superiority of the former was established.Grid-velocity correction was introduced for alleviating the effects of velocity anisotropy in the FDTD implementation of our method, and was shown to play a vital role in reducing the errors due to spatial and temporal discretization [see Fig. 5(b)].

Fig. 1 (Fig. 1 .
Fig. 1.The geometry of the problem.(a) The optical system and the simulation space.(b)The polarization of the electric field on an incident ray (AP) and the corresponding ray in the image space after refraction through the lens (PF).

(Fig. 3 .
Fig. 3. Grayscale plots of the x component of the electric field (in dB) belonging to the focused light pulse evaluated in the x = 0 plane at different time instants.(a) -(f) t = 7τ to t = 17τ in 2τ intervals.(Media 1)