1. INTRODUCTION

The finite-difference time-domain method (FDTD) [1] is a numerical electromagnetic simulation technique that is commonly used in a wide variety of applications [2–6]. One of the main advantages of the FDTD method is its ability to solve Maxwell’s equations in the time domain, thus allowing simulation of the propagation of broadband electromagnetic fields through scattering media. In the FDTD method, space and time are discretized, and the electric and magnetic fields are calculated on a 3D rectangular grid. The FDTD method is an initial value problem, where electric and magnetic fields known at an initial point in time are used to iteratively calculate the field at a later point in time. Focused beams are introduced into the simulation through equivalent electric and magnetic current densities that produce the desired beam profile [7]. This requires that the current densities, and therefore the field associated with the focused beam, are computed at each time step of the FDTD simulation [1,7] in a plane $z = {z_s}$ transverse to the optical axis called the source interface plane. This calculation can be considerable since it generally involves the calculation of a diffraction integral at multiple wave numbers, followed by the evaluation of a Fourier transform to calculate the field in the time domain.

The pseudospectral time-domain method (PSTD) [8], is a memory-efficient variant of FDTD that employs the discrete Fourier transform to calculate the spatial derivatives of the electromagnetic fields instead of central differences. The use of spectrally evaluated spatial derivatives allows fields to be sampled at near the Nyquist rate, which allows for sparser sampling than is possible using the FDTD method. This allows for coarser grids to be used, thus allowing physically larger computational volumes to be simulated using the PSTD method compared with the FDTD method.

In what follows, we refer to the field introduced at the source interface plane, at each time step, as the “incident wave source.” The simplest incident wave source used in the FDTD/PSTD method is the plane wave. However, many optical applications require the simulation of more complicated illumination beams. In this paper we focus on the simulation of focused Gaussian pulses in the focal region. This choice is motivated by the fact that pulsed beams allow the simulation of a broadband response using a single FDTD/PSTD simulation.

Several works propose approximate approaches to calculating the field associated with focused beams, such as those based on the Debye–Wolf integral [9,10], which describes the electromagnetic field distribution in the focal region of high-NA objective lenses [11–21]. As will be explained later, our goal is to develop a technique that does not require multiple monochromatic components to be stored in order to simulate pulsed focal fields. This is because such approaches lead to unfeasible storage and computational loads in FDTD/PSTD simulations.

Even though our interest is primarily in approximating time-domain simulations of focused broadband incident fields, it is instructive to mention some works where approximations and calculations of focused beams have been made in the frequency domain.

One such example introduces an eigenfunction expansion of the electric field in the focal region of the lens [11,12]. In another work, Török et al. [13] compared four approximations of the Debye–Wolf integral, which represent the directly evaluated integral as a sum of simpler integrals, and showed that the direct integration is the fastest method. Leutenegger et al. [14] calculated the focused beam by applying the two-dimensional Fourier transform, achieving good accuracy and maintaining a limited computational time. More recently, a generalization of the Debye–Wolf integral that includes any kind of aberration has been proposed by Wang et al. [15].

There have been a number of studies that employ the FDTD/PSTD method to simulate time-domain fields in the focal region of a converging lens. Davidson and Ziolkowski [16] created a model for introducing a focused beam for rotationally symmetric linear-optics problems, reducing the problem to two dimensions (if rotationally symmetric samples are employed). Çapoglu et al. [17] described a procedure to represent a focused pulse as a finite sum of plane waves by calculating an angular spectrum decomposition of each monochromatic component. They also generalized this method to general laser transverse-electric-magnetic (TEM) modes [18]. This method can be used to simulate electric and magnetic fields with good accuracy even if high-NA lenses are employed, and it is faster than direct integration. Despite those advantages, Çapoglu et al.’s approach remains computationally demanding, because several plane waves must be summed in order to accurately represent the Gaussian pulse. Çapoglu et al.’s method has inspired other subsequent works. Singh et al. made a model [19] that extended Çapoglu et al.’s work [17,18] by introducing dispersion and polarization compensation of the beam. Even Bessel beams have been analytically approximated with the same approach by Wu et al. [20,21].

In this paper we introduce a new approximate method for calculating the time-domain fields associated with a focused pulse, which significantly improves the computational efficiency of the FDTD/PSTD method. This new method employs a sixth-order Taylor expansion of the inner function in the Debye–Wolf integral, and it is valid only for low numerical apertures.

The remainder of the paper is organized as follows. In Section 2, we introduce the rigorous definition of a focused beam as the incident wave source for FDTD/PSTD methods. In Section 3, we summarize an existing approximation of the incident wave source, which uses the field at the central frequency to simulate the focused beam of a broadband incident wave source [2]. In Section 4, we describe a new analytical approximation of the incident field in the time domain, which is based on calculating the Taylor expansion of the inner function of the Debye–Wolf integral. In Section 5, we compare the performances of the two approximations of the incident wave source for some example applications.

 

Fig. 1. Typical optical system to which the Debye–Wolf formula [Eq. (1)] can be applied. The function $\widetilde \phi ({{u_x},{u_y}})$ specifies the profile of the field incident upon the aperture, which is mapped by the lens to the function $\phi ({{s_x},{s_y}})$ on the Gaussian reference sphere. Each point $({{u_x},{u_y}})$ uniquely corresponds to the vector $\Big({{s_x},{s_y},\sqrt {1 - s_x^2 - s_y^2}}\Big)$, which defines the direction of a particular ray in the sample space (see [22] for more details).

Download Full Size | PDF

2. RIGOROUS DEFINITION OF THE INCIDENT FIELD FOR FDTD/PSTD METHODS

Although this paper is focused on two approximate methods, we must first start with an introduction to the rigorous formulation of the incident field. This will enable us to assess the relative accuracy and performance of the approximate methods. With the aid of Fig. 1, let us start by defining the monochromatic component of the electric field at a point $(x,y,z)$ in the vicinity of the focus of a lens, which is calculated with the Debye–Wolf integral [9,10,23]

Given that the FDTD/PSTD method simulates the propagation of the electric and magnetic fields in the time domain, the incident field must also be introduced in the FDTD/PSTD computational space in the time domain at every iteration of the method. The source fields are usually introduced on a plane transverse to the optical axis (usually the $z$ direction) in a region free of scatterers. Irrespective of whether the FDTD or PSTD method is employed, the monochromatic components of the focused pulse are calculated by modulating monochromatic fields found using Eq. (1). Namely, a pulse with a Gaussian shape in the time domain can be described by the following modulation terms:

In general, the calculation of Eq. (4b) has a large computational cost because for each point $(x,y,z)$ and each time step $t$, three integrals must be calculated [two integrals from Eq. (1) and one from the Fourier transform]. This calculation can be sped up by using the fast Fourier transform (FFT) [25] to avoid the integration related to the Fourier transform. However, in order to apply the FFT, all monochromatic components must be calculated and stored in advance for each grid point on the plane where the source is introduced, which requires a substantial amount of computer memory. This paper has been motivated by the need to reduce the computational load of Eq. (4b), which must be calculated for each time step within the FDTD/PSTD algorithm.

3. APPROXIMATION BASED ON A INCIDENT WAVE SOURCE THAT EMPLOYS THE CENTRAL FREQUENCY OF THE SPECTRUM

In this section we review an existing incident wave source, which is an approximation of the rigorous incident wave source that we will denote the central frequency approximation [2,26]. We consider this existing approximation for two reasons. The first reason is to understand the properties of this incident wave source. The second reason is in order to provide a comparison with the new incident wave source, which will be introduced in the next section. The central frequency approximation calculates only the field ${\boldsymbol E}({\boldsymbol r},z;{\nu _0})$ at the spectrum’s central frequency ${\nu _0}$. This considerably reduces the amount of computer memory required to perform the simulation. This single monochromatic component is used to generate an approximate broadband incident waveform, which is introduced into the FDTD/PSTD algorithm.

We use the superscript $C$ to denote fields related to the central frequency approximation. Namely, ${{\boldsymbol E}^C}({\boldsymbol r},z,t)$ and ${\tilde {\!\boldsymbol E}^C}({\boldsymbol r},z;\nu)$ are the time-domain and time-harmonic incident fields associated with the central frequency approximation, respectively. ${\tilde {\!\boldsymbol E}^C}({\boldsymbol r},z;\nu)$ is defined as the multiplication of the modulation in spectral domain $\tilde H(\nu)$ of the central monochromatic component [Eq. (1) with $\nu = {\nu _0}$]. Each Cartesian component $\tau = x,y,z$ of the approximate incident field in the frequency and time domains (at the source interface $z = {z_s}$) is given by

We emphasize that this approximation requires only a single frequency component of the field to be calculated and stored, thus avoiding the storage of several monochromatic components. This incident wave source [Eq. (5b)] is introduced into the FDTD/PSTD grid (at the source interface $z = {z_s}$) as a computationally efficient approximation of the rigorous case [see Eq. (4b)]. This approximation has the useful property that the temporal and spatial dependencies are separable, so we need to update only the time-dependent function and not the space-dependent function at each time step of the FDTD/PSTD algorithm, dramatically reducing the computational load relative to the rigorous case.

A. Angular Spectrum Analysis

In order to analyze how the central frequency approximation perturbs the focused beam relative to the rigorous case, we compare the angular spectrum of the rigorous and approximate fields at the source interface $z = {z_s}$ where the approximation is defined.

We refer readers to Supplement 1 for the angular spectrum definition [27] and for its calculation in the rigorous $({\cal{\tilde A}_\tau}^R({\boldsymbol \alpha},{z_s};\nu))$ and approximate $({ \cal{\tilde A}_\tau}^C({\boldsymbol \alpha},{z_s};\nu))$ cases. The following equation shows the relation between the two angular spectra (see Supplement 1):

Given that $\Big({{\alpha _x},{\alpha _y},\sqrt {1 - \alpha _x^2 - \alpha _y^2}}\Big)$ is the direction of propagation of a plane wave component of the angular spectrum, the meaning of Eq. (6) is that, for a given frequency $\nu$, each plane wave component of the angular spectrum of the approximate case is given by a rescaled version of the same plane wave related to the rigorous case, whose propagation direction has been modified from $\Big({{\alpha _x},{\alpha _y},\sqrt {1 - \alpha _x^2 - \alpha _y^2}}\Big)$ to $\bigg({\frac{\nu}{{{\nu _0}}}{\alpha _x},\frac{\nu}{{{\nu _0}}}{\alpha _y},\sqrt {1 - {{({\frac{\nu}{{{\nu _0}}}{\alpha _x}})}^2} - {{({\frac{\nu}{{{\nu _0}}}{\alpha _y}})}^2}}}\bigg)$.

Figure 2 shows a plot of the relative error at the source interface of the angular spectrum of the incident wave source of the central frequency approximation relative to the rigorous case. It is clear from the figure that the relative error increases as a function of the distance of the interface plane from the focus and also increases as $| {\lambda - {\lambda _0}} |$ (or $| {\nu - {\nu _0}} |$) increases, showing that the approximation is not reliable for values of $\nu$ far from ${\nu _0}$.

 

Fig. 2. Plot of the angular spectrum error [see Eq. (S5) of Supplement 1] as a function of $\lambda = \frac{c}{\nu}$, calculated in four cases of source interface plane $z = {z_s}$, where ${z_s} = 0\,\,\unicode{x00B5}{\rm m}, - 50\,\,\unicode{x00B5}{\rm m}, - 100\,\,\unicode{x00B5}{\rm m}, - 200\,\,\unicode{x00B5}{\rm m}$, for a numerical aperture ${\rm NA} = 0.097$. The smallest error is in focus ($z = 0\,\,\unicode{x00B5}{\rm m}$), and the error increases as ${z_s}$ increases when $\lambda$ is far from the wavelength related to the central frequency ${\lambda _0} = \frac{c}{{{\nu _0}}} = 1300\;{\rm nm} $. The angular spectrum for the rigorous and approximate cases have been calculated with Eqs. (S2) and (S4) of Supplement 1.

Download Full Size | PDF

4. APPROXIMATION BASED ON A INCIDENT WAVE SOURCE THAT EMPLOYS THE TAYLOR EXPANSION

In this section we introduce an alternative approximation of the rigorous wave source that employs a sixth-order Taylor expansion of the inner integrand function of Eq. (4b) to find an approximate analytical solution for focused fields in the time domain, which does not require any integrals to be evaluated. We will refer to this as the Taylor-based approximation. In order to construct this incident wave source we first need to explicitly rewrite the rigorous incident field in the time domain (details in Supplement 1):

In order to make Eq. (10) analytically integrable, we need to have a second-order polynomial in the exponent, in which case the integrand function becomes a multiplication of a polynomial and a Gaussian function having a complex argument. By manipulating the second order of the Taylor expansion to the exponent term $T_g^{\,6}({{s_x},{s_y}})$, we obtain

The exponent $T_g^2({{s_x},{s_y}})$ in the inner integral is a second-degree polynomial of the variable ${s_y}$ and can be rewritten as $T_g^2({{s_x},{s_y}}) = c({s_x}) - {({\alpha ({s_x}){s_y} + \beta ({s_x})})^2}$, where all parameters $c({s_x}),\alpha ({s_x}),\beta ({s_x})$ are polynomials in the variable ${s_x}$. By rearranging the variables in this way, we can rewrite the inner integral as

From the simulations that we have run we have seen that the Gaussian terms $[{{p_m}({s_x}){e^{- {{({\alpha ({s_x})z + \beta ({s_x})})}^2}}}}]_{z = - \sqrt {{{({\frac{{\rm NA}}{\eta}})}^2} - s_x^2}}^{z = + \sqrt {{{({\frac{{\rm NA}}{\eta}})}^2} - s_x^2}}$ are much smaller in amplitude than the terms related to the ${\rm erf}$ function, so we can neglect those terms. We approximate the ${\rm erf}$ function terms of Eq. (14) with their sixth-order Taylor expansion in the variable ${s_x}$, so that Eq. (14) becomes

By substituting Eq. (15) into Eq. (12), we can calculate an approximate analytical solution

Keeping only the terms having the exponent of $s_x^{n + m}$ smaller or equal to 6, and rearranging the terms of Eq. (16b) such that $k = n + h$, we find the analytical approximation of the electric incident field:

Given that the term $c({s_x})$ is a second-degree polynomial, the integrals appearing in Eq. (17) can be solved analytically by calculating the complex ${\rm erf}$ function as we have done in Eq. (14). Given that the evaluation of the complex ${\rm erf}$ function can be computationally expensive, we have approximated that function as follows:

A. Practical Implementation of the Taylor-Based Approximation

The accuracy of the approximation of Eq. (17) depends on both the spatial location and time instant at which the field is evaluated. The relative error between the rigorous and Taylor-based incident wave source is minimized when $x = 0$ or $y = 0$. In particular, the terms containing ${s_x}$ and ${s_y}$ to the first order vanish when $x = 0$ and $y = 0$, respectively, from ${T_s} = t - {t_0} - \frac{\eta}{c}\bigg({x{s_x} + y{s_y} + z\sqrt {1 - {{(s_x^2 + s_y^2)}^2}}}\bigg)$ in Eq. (7). Under this condition, ${T_s}$ is better approximated by the Taylor expansion. When $x = 0$ or $y = 0$, the variable of integration of the inner integral of Eq. (16a) can be chosen to maximize the efficiency of the approximation. In particular, the inner integral of Eq. (16a) should be the ${s_\tau}$ variable associated with the null coordinate ($x$ or $y$). Equation (16a) has been derived by choosing the inner integral related to ${s_y}$, which is more suitable for the case $y = 0$. On the other hand, in the case $x = 0$, it is better to choose the inner integral of Eq. (16a) related to ${s_x}$ and the outer integral related to ${s_y}$.

The most efficient way to use the Eq. (17) in a FDTD/PSTD simulation is to evaluate the Taylor-based approximation at the points $(x,0,z,t)$ and $(0,y,z,t)$, and we use the radially symmetric properties of the integral in Eq. (17) to evaluate the field at all points. The approach we follow is based on existing works [9,10] and is described in Supplement 1. This approach allows us to rewrite the rigorous field in Eq. (7) as

In order to maintain this symmetry in the Taylor-based approximation, we calculate the Taylor-based approximation at the points $(\rho ,0,z,t)$ and $(0,\rho ,z,t)$, and we approximate the functions ${I_0}$ and ${I_2}$ by using Eqs. (21a) and (21b), where $E_x^R$ is replaced with $E_x^T$. After that, we use Eq. (19) to calculate the Taylor-based approximation at each point in space.

We have seen that the efficiency of the Taylor-based approximation can be further enhanced by modifying the integration domain $\Omega = \{{({s_x},{s_y}) \in {\rm I}{{\rm R}^2}|\sqrt {{\rm s}_{\rm x}^2 + {\rm s}_{\rm y}^2} \lt \frac{{{\rm NA}}}{\eta}} \}$ during the calculation of the points $E_x^R(\rho ,0,z,t)$ and $E_x^R(0,\rho ,z,t)$. In the case $x = 0$, the integration domain is modified to

5. COMPARISON OF THE TWO APPROXIMATIONS OF THE RIGOROUS INCIDENT WAVE SOURCE

In this section we compare the approximations presented in the previous sections. There are some physical parameters related to the simulation of the broadband illumination that must be chosen in advance, including ${\nu _0}$, $\eta$, $W$, the frequency bandwidth of the beam $\Delta \nu$, and the numerical aperture NA. In all examples shown here, we have chosen the parameters related to a Thorlabs Telesto-II Spectral Domain OCT Imaging System [2]: ${\nu _0} = 2.3061 \times {10^{14}}\;{\rm Hz}$ (giving a wavelength in air related to the central frequency of ${\lambda _0} = \frac{c}{{{\nu _0}}} = 1300\;{\rm nm} $) and $\Delta \nu = 3.015 \times {10^{13}}\;{\rm s^{- 1}}$ (corresponding to a wavelength width of 170 nm in air), $\eta = 1.42$, ${\rm NA} = 0.097$, $W = 1.44 \times {10^{13}}\; {\rm Hz}$.

A. Relative Error at the Source Interface

The first comparison is made at the source interface (i.e., $z = {z_s}$), where all fields are known analytically. Thus, we do not use the FDTD/PSTD algorithm for this comparison. The central frequency approximation has the advantage that it agrees very closely with the rigorous field when it is calculated near ${\nu _0}$, and it is identical in the case $\nu = {\nu _0}$. However, the error can be significant at frequencies far from ${\nu _0}$ (as shown in Fig. 2). By comparison, the error of the Taylor-based approximation is small, and it is accurate even for frequencies far from ${\nu _0}$. For this reason, the incident wave source related to that approximation is superior to the wave source related to the central frequency approximation for broadband simulations. Figure 3 compares the approximations for several different values of ${z_s}$ using the following error metric:

 

Fig. 3. Plot of the error metric of the central frequency and Taylor-based approximations [see Eq. (23)] on a log scale, evaluated at the source interface plane $z = {z_s}$, for four different values of ${z_s}$.

Download Full Size | PDF

In order to further validate the accuracy of the approximations, we have added a subsection to Section 5 of Supplement 1 where two additional error metrics are considered. These alternative error metrics give similar results to Eq. (23) and suggest that the Taylor-based approximation is more accurate than the central frequency approximation.

B. Relative Error after Propagating the Incident Field with the PSTD Algorithm

Figure 4 shows the integrated relative error [see Eq. (23)] of the electric fields obtained using both approximations, after having propagated all of the three wave sources (${{\boldsymbol E}^R},{{\boldsymbol E}^T},{{\boldsymbol E}^C}$) with the PSTD algorithm, from the source plane $z = {z_s}$ to another transverse plane $z = {z_F}$ (see Fig. 5). The PSTD simulations employed a spatial step size ${\Delta _x} = \frac{{{\lambda _0}}}{4}$ and a time step of ${\Delta _t} = \frac{1}{{2\sqrt 3}}\frac{{\eta {\Delta _x}}}{c}$ (which is about 25% smaller than the maximum time step allowed by the stability criterion for PSTD algorithm [8]). This choice of ${\Delta _t}$ results in numerical dispersion, so that even the numerical propagation of the rigorous field ${{\boldsymbol E}^R}$ will be slightly different from the field at the plane $z = {z_f}$ calculated analytically using the Debye–Wolf integral [see Eq. (1)]. Numerical dispersion does not, however, complicate the comparison, since the numerical dispersion acts in a similar way for the rigorous and approximate cases.

 

Fig. 4. Plots of the integrated relative error [see Eq. (23)] of the two approximations after the incident field has been propagated by the PSTD algorithm (see Fig. 5) from the source interface $z = {z_s}$, where ${z_s} = - 50\,\,\unicode{x00B5}{\rm m}$ for a distance of $\Delta z\unicode{x00B5}{\rm m}$.

Download Full Size | PDF

 

Fig. 5. Representation of the rigorous incident field (red) and an its approximation (the image applies to both Taylor-based and central frequency approximation). The approximate field is introduced at the source plane ($z = {z_s}$), and both fields are propagated to and evaluated at the plane ($z = {z_F}$).

Download Full Size | PDF

For this test the field was introduced at the plane $z = {z_s} = - 50\,\,\unicode{x00B5}{\rm m}$. Then, the rigorous and the two approximate source fields have been propagated using the PSTD algorithm to 3 transverse planes located at $z = {z_F} = {z_s} + {\Delta _z}$, where ${\Delta _z} = 0,50\,\,\unicode{x00B5}{\rm m},100\,\,\unicode{x00B5}{\rm m}$. For each case, the relative error [see Eq. (23)] has been calculated. Figure 4 shows that the error does not increase with the propagation distance of the field. For this reason, the error of both approximations is principally related to the error of the fields at the source interface. As it can be seen in Fig. 4, the field related to the central frequency approximation remains inaccurate across the spectrum after the associated wave source is numerically propagated with the PSTD algorithm.

In Section 5 of Supplement 1, we have plotted similar figures to Fig. 3 and to Fig. 4 with two alternative error metrics, which give an estimation of the error as a function of the positive $x$ axis.

C. Computational Time and RAM Occupied to Calculate the Incident Field in the Rigorous and Approximate Simulations

In this section we compare the rigorous incident field and the two approximations in terms of computational time and memory occupied. For each case, the computational time is the total time needed to calculate the incident field for all time steps of the PSTD algorithm. The occupied memory includes the RAM related to all variables employed in the calculation of each incident field.

Table 1 shows data for the three different ways used to calculate the incident field, in particular, the rigorous, central frequency approximation, and Taylor-based approximation. For each simulation, the computational time and RAM required to perform the computation are shown. For both cases (${\Delta _x} = \frac{\lambda{_0}}{4}$ and ${\Delta _x} = \frac{\lambda{_0}}{{20}}$), computational time and occupied memory are much smaller for the approximations [(C) and (T)] than the rigorous case (R). In the case ${\Delta _x} = \frac{{{\lambda _0}}}{{20}}$, it is clear that the approximation can save many hours of simulation and many gigabytes (GBs) of RAM.

Table 1. Table Showing the Computational Time Required to Calculate the Introduction of the Entire Incident Field in the PSTD Algorithm at the Source Interface $z = - 50\,\,\unicode{x00B5}{\rm m}$, for the Rigorous Case (R), Central Frequency Approximation (C), and Taylor-Based Approximation (T)^a

View Table

D. Limitations of the Taylor-Based Approximation

The principal limitation of the Taylor-based approximation is related to the filling parameter $F$ of the Gaussian beam at the aperture plane, which specifies the fraction of the beam that is inside the physical aperture of the objective lens. From the Gaussian component of the Debye–Wolf integral [see Eq. (1)], we calculate the filling parameter as $F(\nu)$ such that ${e^{- {{({\frac{\nu}{W}|{\boldsymbol s}|})}^2}}} = {e^{- {{({\frac{F}{{{\rm NA}/\eta}}|{\boldsymbol s}|})}^2}}}\; \Rightarrow \; F(\nu) = \frac{{{\rm NA}}}{\eta}\frac{\nu}{W}$. If $F(\nu) \gt 1$, the aperture is underfilled, which means that a small diameter Gaussian fits well within the aperture. If $F(\nu) \lt 1$, the aperture is overfilled, which means that much of the beam is not transmitted through the aperture. In Fig. 6 we have plotted the averaged error of Eq. (23) over all the frequencies of the spectrum $({\overline {{\rm Err}} (z): = \frac{1}{{\Delta \nu}}\int_{{\nu _0} - \frac{{\Delta \nu}}{2}}^{{\nu _0} + \frac{{\Delta \nu}}{2}} {\rm Err}({\nu ,z}){\rm d}\nu})$ for several values of NA and $F({\nu _0})$. Each pair $(F({\nu _0}),{\rm NA})$ is related to a single simulation, and the error is showed as a function of $\frac{{{\rm NA}}}{\eta}$, for four cases of $F({\nu _0})$. Both fields (rigorous and Taylor-based) have been calculated at the focus (${z_s} = 0$). Figure 6 shows that the rigorous incident field is well approximated by the Taylor-based approximation in the underfilled case, and the error increases as a function of the numerical aperture. The error of the central frequency approximation (right side of Fig. 6) is roughly independent of the choice of NA, but even in this case the filling parameter $F({\nu _0})$ significantly affects the accuracy of the approximation.

 

Fig. 6. Plots of the average over the frequencies of the error [Eq. (23)] of the Taylor-based [on the left (T)] and central frequency [on the right (c)] approximations as a function of $\frac{{\rm NA}}{\eta}$, for four cases of filling parameter. Each field has been calculated in focus (${z_s} = 0$).

Download Full Size | PDF

Another important limitation of our approximation is that it can be employed only for time-domain focused pulses whose monochromatic components are fundamental Gaussian beams $({{{\rm TEM}_{(0,0)}}})$ on the Gaussian reference sphere of the lens, unlike the central frequency approximation or Çapoglu et al.’s approach [18], which can be generalized to any transverse-electric-magnetic mode $({{{\rm TEM}_{(n,m)}}})$.

6. CONCLUSION

In this paper we have studied two incident wave sources that are approximations of the rigorous incident field for FDTD/PSTD methods for simulations of focused Gaussian pulses. We have analyzed the central frequency approximation that employs the complex amplitude calculated at only the central frequency of a spectrum to approximate a focused broadband beam in the time domain. We have also introduced a new approximation which employs the Taylor expansion to approximate the rigorous incident field. This new analytical approximation is accurate for underfilled apertures but remains accurate in the overfilled case for low numerical apertures. We showed that the error related to the central frequency approximation increases for frequencies far from the central value, while the Taylor-based approximation well approximates the rigorous field for all frequencies in the spectrum, and therefore the Taylor-based approximation is a reliable approximation of the rigorous field.