## A compact source condition for modelling focused fields using the pseudospectral time-domain method |

Optics Express, Vol. 22, Issue 5, pp. 5599-5613 (2014)

http://dx.doi.org/10.1364/OE.22.005599

Acrobat PDF (2051 KB)

### Abstract

The pseudospectral time-domain (PSTD) method greatly extends the physical volume of biological tissue in which light scattering can be calculated, relative to the finite-difference time-domain (FDTD) method. We have developed an analogue of the total-field scattered-field source condition, as employed in FDTD, for introducing focussed illuminations into PSTD simulations. This new source condition requires knowledge of the incident field, and applies update equations, at a single plane in the PSTD grid. Numerical artifacts, usually associated with compact PSTD source conditions, are minimized by using a staggered grid. This source condition’s similarity with that used by the FDTD suggests a way in which existing FDTD codes can be easily adapted to PSTD codes.

© 2014 Optical Society of America

## 1. Background

5. M. Ding and K. Chen, “Staggered-grid PSTD on local Fourier basis and its applications to surface tissue modeling,” Opt. Express **18**, 9236–9250 (2010). [CrossRef] [PubMed]

6. S. H. Tseng, Y. L. Kim, A. Taflove, D. Maitland, V. Backman, and J. T. Walsh, “Simulation of enhanced backscattering of light by numerically solving Maxwell’s equations without heuristic approximations,” Opt. Express **13**, 3666–3672 (2005). [CrossRef] [PubMed]

7. G. Chen, P. Yang, and G. Kattawar, “Application of the pseudospectral time-domain method to the scattering of light by nonspherical particles,” J. Opt. Soc. Am. A **25**, 785–789 (2008). [CrossRef]

8. S. H. Tseng, J. H. Greene, A. Taflove, D. Maitland, V. Backman, and J. T. Walsh, “Exact solution of Maxwell’s
equations for optical interactions with a macroscopic random
medium,” Opt. Lett. **29**, 1393–1395
(2004). [CrossRef] [PubMed]

9. T. -W. Lee and S. C. Hagness, “A compact wave source condition for the pseudospectral time-domain method,” IEEE Antenn. Wirel. Pr. **3**, 253–256 (2004). [CrossRef]

10. X. Gao, M. Mirotznik, and D. Prather, “A method for introducing soft sources in the PSTD algorithm,” IEEE T. Antenn. Propag. **52**, 1665–1671 (2004). [CrossRef]

11. D. Merewether, R. Fisher, and F. Smith, “On implementing a numeric Huygen’s source scheme in a finite difference program to illuminate scattering bodies,” IEEE T. Nucl. Sci. **27**, 1829–1833 (1980). [CrossRef]

## 2. Field equivalence

12. J. A. Stratton and L. J. Chu, “Diffraction theory of electromagnetic waves,” Phys. Rev. **56**, 99–107 (1939). [CrossRef]

13. P. Petre and T. Sarkar, “Planar near-field to far-field transformation using an equivalent magnetic current approach,” IEEE T. Antenn. Propag. **40**, 1348–1356 (1992). [CrossRef]

*ω*(assuming the exp(−

*iωt*) convention): where a linear isotropic material of permittivity

*ε*and permeability

*μ*has been assumed,

**J**

^{*}is the magnetic current density and

*ρ*

^{*}the magnetic charge density. All other terms take their usual meaning. Neither

**J**

^{*}or

*ρ*

^{*}are required to be included but they make Maxwell’s equations symmetric and allow the equivalence principle to be developed. Perhaps the most widely known result from Stratton and Chu’s paper, developed using only vector identities and equations of continuity, is that the electric field within a closed surface,

*S*, containing no sources or sinks of radiation, is given by: where (

*x′*,

*y′*,

*z′*) is an interior point of

*S*,

*G*= exp(

*ikr*)/

*r*is the free space Green’s function,

*r*is the distance between (

*x′*,

*y′*,

*z′*) and a point on

*S*, d

*a*is a differential surface element,

**n**is an outward facing surface normal and

**E**and

**H**appearing inside the integral are the electric and magnetic fields on

*S*, respectively. A complementary result is also obtained for the magnetic field. A further important result of Stratton and Chu, which connects Eq. (2) to Eqs. (1), is that exactly the same field will result within

*S*, as illustrated in Fig. 1, if the field incident upon

*S*vanished but the following equivalent current

*surface*densities and charge existed on

*S*: where

**J**

*is an equivalent electric surface current density,*

_{s}*η*an equivalent surface electric charge density. This result is usually referred to as the field equivalence principle.

## 3. The finite-difference time-domain method

*y*direction. The concepts presented here are readily extended to three-dimensions. By choosing the so-called transverse magnetic (TM) independent subset of Maxwell’s equations, we obtain the following equations which govern the propagation of electromagnetic waves in the computational domain [17]: where

*ε*and

*μ*are the electric permittivity and magnetic permeability of free space, respectively,

*J*and

_{x}*J*are electric source current densities. A set of recurrence relations is obtained from these equations by discretising

_{z}*E*,

_{x}*E*and

_{z}*H*in space and time and approximating temporal and spatial derivatives with central difference equations. The simulation proceeds from an initial condition and ceases when a final condition is reached. The difference equations resulting from Eqs. (4) are [17]:

_{y}*t*and Δ

*are the temporal and spatial grid spacings, respectively, which are subject to restrictions. The indexing used here is such that*

_{xz}*H*(

_{y}*i*,

*k*;

*n*) refers to the value of

_{t}*H*at position (

_{y}*i*Δ

*,*

_{xz}*k*Δ

*) at time*

_{xz}*n*Δ

_{t}*t*. Note also that the field components are offset by half a grid spacing in both space and time.

19. K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE T. Antenn. Propag. **14**, 302–307 (1966). [CrossRef]

*J*and

_{x}*J*, appearing in Eqs. (5)–(7) to introduce an incident wave at a fictitious surface at each iteration of the FDTD algorithm [11

_{z}11. D. Merewether, R. Fisher, and F. Smith, “On implementing a numeric Huygen’s source scheme in a finite difference program to illuminate scattering bodies,” IEEE T. Nucl. Sci. **27**, 1829–1833 (1980). [CrossRef]

*S*=

*S*

_{1}∪

*S*

_{2}∪

*S*

_{3}∪

*S*

_{4}. Alternatively, if

*S*

_{1}is sufficiently large such that

*S*

_{2}and

*S*

_{3}are pushed to infinity and

*S*

_{4}is also pushed to infinity, the source condition can be implemented on

*S*

_{1}only since

*S*

_{2},

*S*

_{3}and

*S*

_{4}contribute negligibly to the field in the total field region [20]. This enables a focused illumination to be accurately introduced along a single interface, as in the right hand diagram of Fig. 2 as long as the majority of the focused beam’s energy propagates through

*S*

_{1}. In both cases, the incident field is introduced to the update equations (Eqs. (5)–(7)) via the source magnetic and electric current densities along the surface

*S*, in the left hand case, and along the surface

*S*

_{1}in the right hand case. The region in which the incident field is present is denoted the total field region, whilst that in which only the scattered field is present is the scattered field region.

*H*(

_{y}*i*,

*k*;

_{s}*n*), Eq. (5), in the region of the interface between total and scattered field regions at

_{t}*k*Δ

_{s}*, as depicted for*

_{xz}*i*= 0 in Fig. 2. If Eq. (5) were evaluated directly without accounting for the fact that one

*E*term is a scattered quantity and the other a total quantity, the update equation would be inconsistent. However, since the incident field

_{x}*E*(

_{x,i}*i*,

*k*+ 1/2;

_{s}*n*+ 1/2) is known analytically, an additional update equation can be executed to make the original update equation consistent [17], taking the form: where we have used the notation of Taflove and Hagness [17] to indicate that

_{t}*H*(

_{y}*i*,

*k*;

_{s}*n*)

_{t}_{(5)}is the quantity arising from the evaluation of Eq. (5). This can also be interpreted as including a source term

11. D. Merewether, R. Fisher, and F. Smith, “On implementing a numeric Huygen’s source scheme in a finite difference program to illuminate scattering bodies,” IEEE T. Nucl. Sci. **27**, 1829–1833 (1980). [CrossRef]

*E*with respect to

_{x}*z*, at

*z*=

*k*Δ

_{s}*, is able to be corrected locally with knowledge of the analytic incident field. In this two-dimensional example, a similar correction must be made to Eq. (6), the update equation for*

_{xz}*E*.

_{x}## 4. The pseudospectral time-domain method

1. Q. Liu, “The PSTD algorithm: A time-domain method requiring only two cells per wavelength,” Microw. Opt. Techn. Lett. **15**, 158–165 (1997). [CrossRef]

5. M. Ding and K. Chen, “Staggered-grid PSTD on local Fourier basis and its applications to surface tissue modeling,” Opt. Express **18**, 9236–9250 (2010). [CrossRef] [PubMed]

*H*along a line of constant

_{y}*x*at an instant in time. This computational domain has a total field region for

*z*≥ 0 and a scattered region otherwise. In the FDTD case, the update equation for

*E*(Eq. (6)) implicitly calculates

_{x}*∂H*using central differences resulting in an error at

_{y}/∂z*z*= Δ

*/2, as shown in Fig. 3(b), since the difference equation uses one total field quantity and one scattered field quantity as illustrated in Fig. 2. However, this error can be precisely compensated for with the analytically known incident field using Eq. (8). A similar compensation step is required for the*

_{xz}*H*update equation. The crucial point is that derivatives are evaluated locally.

_{y}*∂H*is calculated on a collocated grid using discrete Fourier transforms in Fig. 3(c) which contains a non-zero field in the scattered region and artifacts throughout the computational domain. We explain in more detail why the staggered grid is superior to the collocated grid in Sec. (4.3). Focusing on the staggered grid case, Fig. 3(d) shows that

_{y}/∂z*∂H*is similar to that calculated using central differences but there are two main differences: (1) there are some artifacts around

_{y}/∂z*z*= 0 resulting from the sharp transition between total and scattered fields; and (2) the error at

*z*= Δ

*/2 is now a function of the field everywhere along this particular line of constant*

_{xz}*x*. The consequence of this is that to implement the TFSF source condition within the PSTD method, source update equations must be applied at as many cells adjacent to the interface between total and scattered field regions as is computationally feasible. These updates, however, depend on the field throughout the computational domain, i.e., non-local update equations are required, thus increasing computational and implementation complexity, as well as computer memory requirements.

### 4.1. Existing source conditions

2. Q. Liu, “Large-scale simulations of electromagnetic and acoustic measurements using the pseudospectral time-domain (PSTD) algorithm,” IEEE T. Geosci. Remote **37**, 917–926 (1999). [CrossRef]

9. T. -W. Lee and S. C. Hagness, “A compact wave source condition for the pseudospectral time-domain method,” IEEE Antenn. Wirel. Pr. **3**, 253–256 (2004). [CrossRef]

10. X. Gao, M. Mirotznik, and D. Prather, “A method for introducing soft sources in the PSTD algorithm,” IEEE T. Antenn. Propag. **52**, 1665–1671 (2004). [CrossRef]

21. Z. Lin and L. Thylen, “An analytical derivation of the optimum source patterns for the pseudospectral time-domain method,” J. Comput. Phys. **228**, 7375–7387 (2009). [CrossRef]

22. Q. Li, Y. Chen, and C. Li, “Hybrid PSTD-FDTD technique for scattering analysis,” Microw. Opt. Techn. Lett. **34**, 19–24 (2002). [CrossRef]

21. Z. Lin and L. Thylen, “An analytical derivation of the optimum source patterns for the pseudospectral time-domain method,” J. Comput. Phys. **228**, 7375–7387 (2009). [CrossRef]

10. X. Gao, M. Mirotznik, and D. Prather, “A method for introducing soft sources in the PSTD algorithm,” IEEE T. Antenn. Propag. **52**, 1665–1671 (2004). [CrossRef]

2. Q. Liu, “Large-scale simulations of electromagnetic and acoustic measurements using the pseudospectral time-domain (PSTD) algorithm,” IEEE T. Geosci. Remote **37**, 917–926 (1999). [CrossRef]

9. T. -W. Lee and S. C. Hagness, “A compact wave source condition for the pseudospectral time-domain method,” IEEE Antenn. Wirel. Pr. **3**, 253–256 (2004). [CrossRef]

### 4.2. Discrete Fourier transform operators

*x*is a sequence obtained by sampling a function of a continuous variable. We assume

_{k}*N*to be even since it results in a more computationally efficient evaluation of the discrete Fourier transform and define an index map,

*a*, as Using this index map, we define an operator which allows us to form the optimal interpolation scheme, in terms of minimum deviation [24], as: where we introduce a compact notation to be used in the remainder of this paper, noting that IDFT and DFT refer to inverse discrete Fourier transform and discrete Fourier transform, respectively. Note also that we assume a normalized discrete Fourier transform such that Δ refers to a spatial displacement ΔΔ

_{n}*, where Δ*

_{xz}*is the spatial sampling period.*

_{xz}### 4.3. Staggered grids

1. Q. Liu, “The PSTD algorithm: A time-domain method requiring only two cells per wavelength,” Microw. Opt. Techn. Lett. **15**, 158–165 (1997). [CrossRef]

3. D. Kosloff and E. Baysal, “Forward modeling by a Fourier method,” Geophysics **47**, 1402–1412 (1982). [CrossRef]

25. T. Özdenvar and G. A. McMechan, “Causes and reduction of numerical artefacts in pseudo-spectral wavefield extrapolation,” Geophy. J. Int. **126**, 819–828 (1996). [CrossRef]

26. H. -W. Chen, “Staggered-grid pseudospectral viscoacoustic wave field simulation in two-dimensional media,” J. Acoust. Soc. Am. **100**, 120–131 (1996). [CrossRef]

5. M. Ding and K. Chen, “Staggered-grid PSTD on local Fourier basis and its applications to surface tissue modeling,” Opt. Express **18**, 9236–9250 (2010). [CrossRef] [PubMed]

1. Q. Liu, “The PSTD algorithm: A time-domain method requiring only two cells per wavelength,” Microw. Opt. Techn. Lett. **15**, 158–165 (1997). [CrossRef]

2. Q. Liu, “Large-scale simulations of electromagnetic and acoustic measurements using the pseudospectral time-domain (PSTD) algorithm,” IEEE T. Geosci. Remote **37**, 917–926 (1999). [CrossRef]

7. G. Chen, P. Yang, and G. Kattawar, “Application of the pseudospectral time-domain method to the scattering of light by nonspherical particles,” J. Opt. Soc. Am. A **25**, 785–789 (2008). [CrossRef]

**52**, 1665–1671 (2004). [CrossRef]

22. Q. Li, Y. Chen, and C. Li, “Hybrid PSTD-FDTD technique for scattering analysis,” Microw. Opt. Techn. Lett. **34**, 19–24 (2002). [CrossRef]

27. C. Liu, R. L. Panetta, and P. Yang, “Application of the pseudo-spectral time domain method to compute particle single-scattering properties for size parameters up to 200,” J. Quant. Spectrosc. Ra. **113**, 1728–1740 (2012). [CrossRef]

3. D. Kosloff and E. Baysal, “Forward modeling by a Fourier method,” Geophysics **47**, 1402–1412 (1982). [CrossRef]

25. T. Özdenvar and G. A. McMechan, “Causes and reduction of numerical artefacts in pseudo-spectral wavefield extrapolation,” Geophy. J. Int. **126**, 819–828 (1996). [CrossRef]

26. H. -W. Chen, “Staggered-grid pseudospectral viscoacoustic wave field simulation in two-dimensional media,” J. Acoust. Soc. Am. **100**, 120–131 (1996). [CrossRef]

*x*may be obtained by evaluating: and for a staggered grid as: by substituting Eq. (13) into Eqs. (14) and (15) and evaluating the limits. We can evaluate spatial derivatives as: where and the superscripts

_{k}*c*and

*s*refer to collocated and staggered, respectively, and the ± denotes the derivative being evaluated and interpolated by plus or minus half a spatial sampling step. Equations (18) and (19) can be used to demonstrate a previously published result [25

25. T. Özdenvar and G. A. McMechan, “Causes and reduction of numerical artefacts in pseudo-spectral wavefield extrapolation,” Geophy. J. Int. **126**, 819–828 (1996). [CrossRef]

26. H. -W. Chen, “Staggered-grid pseudospectral viscoacoustic wave field simulation in two-dimensional media,” J. Acoust. Soc. Am. **100**, 120–131 (1996). [CrossRef]

28. G. J. P. Corrêa, M. Spiegelman, S. Carbotte, and J. C. C. Mutter, “Centered and staggered Fourier derivatives and Hilbert transforms,” Geophysics **67**, 1558–1563 (2012). [CrossRef]

*x*is differentiated by convolution with the discrete Fourier transforms of

_{k}*N*= 32. The plots demonstrate the general result that

*n*=

*N*/2, whereas

28. G. J. P. Corrêa, M. Spiegelman, S. Carbotte, and J. C. C. Mutter, “Centered and staggered Fourier derivatives and Hilbert transforms,” Geophysics **67**, 1558–1563 (2012). [CrossRef]

## 5. A new source condition for the PSTD method

28. G. J. P. Corrêa, M. Spiegelman, S. Carbotte, and J. C. C. Mutter, “Centered and staggered Fourier derivatives and Hilbert transforms,” Geophysics **67**, 1558–1563 (2012). [CrossRef]

**E**

*,*

_{i}**H**

*), is focused into a region of interest, bounded by the surface*

_{i}*S*. As has been explained, the equivalence principle allows us to generate the same field, (

**E**

*,*

_{i}**H**

*), within*

_{i}*S*using equivalent surface current densities and surface charge density on

*S*, defined by Eqs. (3), as shown in Fig. 1b). Consider, however, a scenario where the focused field contributes negligibly to

*S*

_{2},

*S*

_{3}and

*S*

_{4}in Fig. 1a). This may be achieved in practice by making

*S*

_{1}large enough such that the incident field is negligible in the vicinity of

*S*

_{2}and

*S*

_{3}whilst moving

*S*

_{4}to infinity where the field is also negligible. We are then left with the scenario depicted in the right hand image of Fig. 2 where

*S*is composed of a single plane. In this case, the surface normal is −

**k̂**, where

**k̂**is the unit vector associated with the

*z*Cartesian coordinate. This gives the equivalent sources on

*S*

_{1}, which, without loss of generality, we assume to be the plane

*z*=

*z*, as: It is, however, possible to eliminate all but the

_{i}*S*

_{1}from below (i.e.,

*z*<

*z*), the electric current and charge densities must vanish as a result of the infinite conductivity, leaving

_{i}**î**and

**ĵ**are the Cartesian unit vectors in the

*x*and

*y*directions, respectively. Summarising this result, the electric field in the region

*z*>

*z*will be equal to

_{i}**E**

*if an equivalent magnetic surface current density equal to 2(*

_{i}**E**

*·*

_{i}**ĵ**, −

**E**

*·*

_{i}**î**, 0) flows on

*S*

_{1}. Another consequence of this approach is that a field will now propagate back into the scattered field region

*z*<

*z*. This field is the mirror image of the incident field and may be removed by subtraction at the conclusion of the PSTD simulation [23].

_{i}*J*and

_{x}*J*) from Eqs. (6) and (7) since, as just explained, we have only a magnetic surface current density present. The second step is to recognise that

_{z}*S*

_{1}and has value −

**E**

*·*

_{i}**î**. The final step is to express the spatial derivatives implicit in Eqs. (5)–(7) using spectral derivatives. For example,

*∂H*is evaluated in Eq. (6) by the terms (

_{y}/∂z*H*(

_{y}*i*,

*k*;

*n*) −

_{t}*H*(

_{y}*i*,

*k*+ 1;

*n*))/Δ

_{t}*, which may be evaluated spectrally using Eq. (17) as*

_{xz}*k*, the PSTD update equations corresponding to Eqs. (5)–(7) may then be expressed as

_{s}*δ*is the Kronecker delta function. Division of

*is necessary to convert from a surface current density to a current density.*

_{xz}- The incident electric field,
**E**, ordinarily employed by the FDTD method should be multiplied by two._{i} - The incident magnetic field ordinarily employed by the FDTD method,
**H**, associated with_{i}**E**should be set to 0._{i} **E**should be evaluated at a distance Δ_{i}/2 in the negative_{xz}*z*direction relative to the location at which the field would be evaluated using the TFSF formulation. This is because the new source condition assumes a sheet of magnetic current density centered on the magnetic field component.

*S*in the left hand image of Fig. 2, since erroneous mirror image fields will, in general, propagate into the total field region. Instead, this method is designed to launch a focused incident field or an apertured incident field. This method is very accurate for focused fields if the plane where the source is introduced is sufficiently wide, such that the majority of the energy in the beam propagates through it [20]. The PSTD method is well suited to this source condition since its reduced sampling requirement enables computational volumes with large cross-sectional areas to be modelled.

## 6. Evaluation and analysis

*μ*m wide and 133

*μ*m deep, using a

*λ*/4 grid spacing and a wavelength of 1325nm. The full aperture of a lens of numerical aperture 0.056 was assumed to focus the beam. The incident beam was calculated using the Richards-Wolf integral [29

29. B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems ii. Structure of the image field in an aplanatic system,” P. R. Soc. A **253**, 358–379 (1959). [CrossRef]

*λ*/4 grid spacing was used, as the PML performance was found to degrade when using larger grid spacings. This simulation required approximately 12 Gb of memory, which is the principal limitation upon the physical volume which can be modelled. If the FDTD method were used to model the same physical volume using a necessarily finer grid spacing of

*λ*/15, over 750 Gb of memory would be required making it unfeasible. Note, however, that our current implementation is far from computationally optimal since we have implemented Berenger’s split fields [17] throughout the entire computational grid for convenience. Remedying this would immediately enable a physical volume 500

*μ*m deep to be modelled and increasing the memory use to 24 Gb enables a depth of 1mm to be reached, deep enough to model optical coherence tomography image formation in biological tissue. Note that the memory usage can be further reduced if the complex amplitudes of the field within the computational volume are not required, as is the case when modelling image formation, as only field values at the surface are necessary [20].

*E*| in which the ripples in the vicinity of the plane where the beam is introduced are apparent. These ripples, however, are restricted to a few cells either side of this plane.

_{x}*x*-polarised plane wave of wavelength 405nm was focused in air by a lens with numerical aperture 0.95. The focused field is calculated using the Richards-Wolf integral [29

29. B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems ii. Structure of the image field in an aplanatic system,” P. R. Soc. A **253**, 358–379 (1959). [CrossRef]

*λ*behind the focus of the beam and was temporally modulated by a Gaussian pulse with a spectral width (full-width at half-maximum) of 60nm so as not to spuriously introduce any high-spatial frequency components. The complex amplitude of the field at the center wavelength was evaluated as the simulation proceeded by discrete Fourier transform. Table 1 summarizes the tests which were performed and the base simulation parameters. The base PSTD simulation employed a grid spacing of

*λ*/2.5, a transverse simulation size of 30

*λ*× 30

*λ*and a time step of Δ

*t*= 1.8511×10

^{−17}s. This time step value is approximately 27 times smaller than the theoretical maximum value of Δ

*t*that guarantees stability but was selected on the basis of the results in Fig. 6, which show this value is necessary to minimize errors due to numerical dispersion, as will be discussed later in this section.

*λ*behind the focus, with the key simulation parameters summarised in Table 1. These results are plotted in Fig. 7 and demonstrate the accuracy of the source condition and that there is no ringing present in the field calculated using the PSTD method. The field is sampled on a grid of spacing

*λ*/2.5 which is too coarse to illustrate all features of the field patterns. Note, however, that this information is readily available using Fourier interpolation.

*E*| as “fringes” which run approximately normal to the contours of equal field magnitude. Detailed investigation revealed that these artifacts arise due to the low resolution at which the images are displayed - not the accuracy of the calculated fields. Furthermore, these visual artifacts are less prominent in the fields calculated using the PSTD method because the high-frequency diffraction rings are not as prominent as in the analytic case. The principal result from this plot is that the PSTD method is able to propagate the focused field accurately based upon a visual comparison with the analytically calculated field.

_{x}*λ*/2. Note that the relative error was calculated for each scenario without re-sampling even though each simulation, having a different grid spacing, samples the field at different locations. Testing, however, showed that any variations arising due to this phenomenon were insignificant compared with the overall trend observed in Fig. 8a), which shows that the estimation of the dominant field component,

*E*, begins to become inaccurate when Δ exceeds approximately

_{x}*λ*/2.5. The next most significant component,

*E*, soon follows whilst

_{z}*E*maintains its accuracy. This result may be understood by reviewing Maxwell’s equations and, in particular, ∇ ×

_{y}**H**+

*iωε*

**E**=

**J**, which may be expressed in the time domain as: We expect the greatest inaccuracy to occur in derivatives evaluated along the

*z*-direction, the nominal direction of propagation of the focused beam and the direction which crosses the plane where the field is introduced. In particular, where the sampling is near the Nyquist rate, any discontinuity in the field is sufficient to cause aliasing which will be more severe the closer the sampling is to the Nyquist rate. Furthermore, noting that the magnetic field components in decreasing order of significance are

*H*,

_{y}*H*and

_{z}*H*, we expect the derivative with the greatest absolute error to be

_{x}*∂H*, which is why

_{y}/∂z*E*exhibits the greatest sensitivity to grid spacing. The most significant term contributing to the

_{x}*E*update equation is

_{z}*∂H*, which introduces an inaccuracy because

_{y}/∂x*H*becomes inaccurate by the same mechanism as

_{y}*E*. The most significant term contributing to

_{x}*E*is

_{y}*∂H*, which becomes inaccurate due to the dependence that

_{z}/∂x*H*has on

_{z}*∂E*. It is, thus, reasonable that the degree of the dependence of each field quantity, upon the least accurately evaluated spatial derivative, determines the overall sensitivity to grid spacing. It is from these results that we conclude a grid spacing of

_{x}/∂y*λ*/2.5 is optimal for this particular simulation. Figure 8a) also reveals that, whilst a minimum error in

*E*of less than 1% can be obtained, the minimum error which can be obtained in

_{x}*E*and

_{y}*E*is less than 2%. We believe that this is the result of there being less energy in these components which biases the error metric.

_{z}*E*(the principal field component) at particular planes along the

_{x}*z*-axis changes with the lateral dimension of the PSTD grid. This is an important parameter since this incident source condition requires the transverse dimension of the grid to be large enough for the majority of the beam’s energy to fall within it. Otherwise, we are implicitly simulating how a focused beam, diffracted by a square aperture, propagates in the grid. The plot in Fig. 8b) shows that the improvement observed by increasing the lateral dimension of the grid diminishes for a grid size of 50

*λ*× 50

*λ*. One minor aspect of this plot is the crossing over of the error lines for the 30

*λ*× 30

*λ*, 40

*λ*× 40

*λ*and 50

*λ*× 50

*λ*cases values of

*z*exceeding 1.5

*μ*m, the cause of which remains a question for future work.

*E*as Δ

_{x}*t*was reduced to below that required for stability. The plot in Fig. 6 shows that a time step in the vicinity of 0.031 times that required for stability (black curve) is necessary to minimize errors due to numerical dispersion, as has been noted by Tseng

*et al.*[30

30. S. H. Tseng, J. H. Greene, A. Taflove, D. Maitland, V. Backman, and J. T. Walsh, “Exact solution of Maxwell’s equations for optical interactions with a macroscopic random medium: addendum,” Opt. Lett. **30**, 56–57 (2005). [CrossRef] [PubMed]

## 7. Conclusions

*λ*/2.5, a time step some 30 times smaller than that required for stability, and a transverse grid size of 50

*λ*× 50

*λ*should be employed to obtain optimal accuracy.

## Acknowledgments

## References and links

1. | Q. Liu, “The PSTD algorithm: A time-domain method requiring only two cells per wavelength,” Microw. Opt. Techn. Lett. |

2. | Q. Liu, “Large-scale simulations of electromagnetic and acoustic measurements using the pseudospectral time-domain (PSTD) algorithm,” IEEE T. Geosci. Remote |

3. | D. Kosloff and E. Baysal, “Forward modeling by a Fourier method,” Geophysics |

4. | J. Virieux and S. Operto, “An overview of full-waveform inversion
in exploration geophysics,”
Geophysics |

5. | M. Ding and K. Chen, “Staggered-grid PSTD on local Fourier basis and its applications to surface tissue modeling,” Opt. Express |

6. | S. H. Tseng, Y. L. Kim, A. Taflove, D. Maitland, V. Backman, and J. T. Walsh, “Simulation of enhanced backscattering of light by numerically solving Maxwell’s equations without heuristic approximations,” Opt. Express |

7. | G. Chen, P. Yang, and G. Kattawar, “Application of the pseudospectral time-domain method to the scattering of light by nonspherical particles,” J. Opt. Soc. Am. A |

8. | S. H. Tseng, J. H. Greene, A. Taflove, D. Maitland, V. Backman, and J. T. Walsh, “Exact solution of Maxwell’s
equations for optical interactions with a macroscopic random
medium,” Opt. Lett. |

9. | T. -W. Lee and S. C. Hagness, “A compact wave source condition for the pseudospectral time-domain method,” IEEE Antenn. Wirel. Pr. |

10. | X. Gao, M. Mirotznik, and D. Prather, “A method for introducing soft sources in the PSTD algorithm,” IEEE T. Antenn. Propag. |

11. | D. Merewether, R. Fisher, and F. Smith, “On implementing a numeric Huygen’s source scheme in a finite difference program to illuminate scattering bodies,” IEEE T. Nucl. Sci. |

12. | J. A. Stratton and L. J. Chu, “Diffraction theory of electromagnetic waves,” Phys. Rev. |

13. | P. Petre and T. Sarkar, “Planar near-field to far-field transformation using an equivalent magnetic current approach,” IEEE T. Antenn. Propag. |

14. | C. Balanis, |

15. | H. Levine and J. Schwinger, “On the theory of electromagnetic wave diffraction by an aperture in an infinite plane conducting screen,” Commun. Pur. Appl. Math. |

16. | R. Harrington, |

17. | A. Taflove and S. Hagness, |

18. | A. Taflove, A. Oskooi, and S. Johnson, eds., |

19. | K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE T. Antenn. Propag. |

20. | P. Török, P. R. T. Munro, and E. E. Kriezis, “High numerical aperture vectorial
imaging in coherent optical microscopes,”
Opt. Express |

21. | Z. Lin and L. Thylen, “An analytical derivation of the optimum source patterns for the pseudospectral time-domain method,” J. Comput. Phys. |

22. | Q. Li, Y. Chen, and C. Li, “Hybrid PSTD-FDTD technique for scattering analysis,” Microw. Opt. Techn. Lett. |

23. | A. Oskooi and S. G. Johnson, “Electromagnetic wave source conditions,” in |

24. | V. Čížek, |

25. | T. Özdenvar and G. A. McMechan, “Causes and reduction of numerical artefacts in pseudo-spectral wavefield extrapolation,” Geophy. J. Int. |

26. | H. -W. Chen, “Staggered-grid pseudospectral viscoacoustic wave field simulation in two-dimensional media,” J. Acoust. Soc. Am. |

27. | C. Liu, R. L. Panetta, and P. Yang, “Application of the pseudo-spectral time domain method to compute particle single-scattering properties for size parameters up to 200,” J. Quant. Spectrosc. Ra. |

28. | G. J. P. Corrêa, M. Spiegelman, S. Carbotte, and J. C. C. Mutter, “Centered and staggered Fourier derivatives and Hilbert transforms,” Geophysics |

29. | B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems ii. Structure of the image field in an aplanatic system,” P. R. Soc. A |

30. | S. H. Tseng, J. H. Greene, A. Taflove, D. Maitland, V. Backman, and J. T. Walsh, “Exact solution of Maxwell’s equations for optical interactions with a macroscopic random medium: addendum,” Opt. Lett. |

**OCIS Codes**

(110.4500) Imaging systems : Optical coherence tomography

(170.3660) Medical optics and biotechnology : Light propagation in tissues

(180.0180) Microscopy : Microscopy

(050.1755) Diffraction and gratings : Computational electromagnetic methods

**ToC Category:**

Scattering

**History**

Original Manuscript: January 2, 2014

Revised Manuscript: February 13, 2014

Manuscript Accepted: February 21, 2014

Published: March 4, 2014

**Virtual Issues**

Vol. 9, Iss. 5 *Virtual Journal for Biomedical Optics*

**Citation**

Peter R.T. Munro, Daniel Engelke, and David D. Sampson, "A compact source condition for modelling focused fields using the pseudospectral time-domain method," Opt. Express **22**, 5599-5613 (2014)

http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-22-5-5599

Sort: Year | Journal | Reset

### References

- Q. Liu, “The PSTD algorithm: A time-domain method requiring only two cells per wavelength,” Microw. Opt. Techn. Lett. 15, 158–165 (1997). [CrossRef]
- Q. Liu, “Large-scale simulations of electromagnetic and acoustic measurements using the pseudospectral time-domain (PSTD) algorithm,” IEEE T. Geosci. Remote 37, 917–926 (1999). [CrossRef]
- D. Kosloff, E. Baysal, “Forward modeling by a Fourier method,” Geophysics 47, 1402–1412 (1982). [CrossRef]
- J. Virieux, S. Operto, “An overview of full-waveform inversion in exploration geophysics,” Geophysics 74, WCC127 (2009). [CrossRef]
- M. Ding, K. Chen, “Staggered-grid PSTD on local Fourier basis and its applications to surface tissue modeling,” Opt. Express 18, 9236–9250 (2010). [CrossRef] [PubMed]
- S. H. Tseng, Y. L. Kim, A. Taflove, D. Maitland, V. Backman, J. T. Walsh, “Simulation of enhanced backscattering of light by numerically solving Maxwell’s equations without heuristic approximations,” Opt. Express 13, 3666–3672 (2005). [CrossRef] [PubMed]
- G. Chen, P. Yang, G. Kattawar, “Application of the pseudospectral time-domain method to the scattering of light by nonspherical particles,” J. Opt. Soc. Am. A 25, 785–789 (2008). [CrossRef]
- S. H. Tseng, J. H. Greene, A. Taflove, D. Maitland, V. Backman, J. T. Walsh, “Exact solution of Maxwell’s equations for optical interactions with a macroscopic random medium,” Opt. Lett. 29, 1393–1395 (2004). [CrossRef] [PubMed]
- T. -W. Lee, S. C. Hagness, “A compact wave source condition for the pseudospectral time-domain method,” IEEE Antenn. Wirel. Pr. 3, 253–256 (2004). [CrossRef]
- X. Gao, M. Mirotznik, D. Prather, “A method for introducing soft sources in the PSTD algorithm,” IEEE T. Antenn. Propag. 52, 1665–1671 (2004). [CrossRef]
- D. Merewether, R. Fisher, F. Smith, “On implementing a numeric Huygen’s source scheme in a finite difference program to illuminate scattering bodies,” IEEE T. Nucl. Sci. 27, 1829–1833 (1980). [CrossRef]
- J. A. Stratton, L. J. Chu, “Diffraction theory of electromagnetic waves,” Phys. Rev. 56, 99–107 (1939). [CrossRef]
- P. Petre, T. Sarkar, “Planar near-field to far-field transformation using an equivalent magnetic current approach,” IEEE T. Antenn. Propag. 40, 1348–1356 (1992). [CrossRef]
- C. Balanis, Advanced Engineering Electromagnetics (John Wiley and Sons, 1989).
- H. Levine, J. Schwinger, “On the theory of electromagnetic wave diffraction by an aperture in an infinite plane conducting screen,” Commun. Pur. Appl. Math. 3, 355–391 (1950). [CrossRef]
- R. Harrington, Time-Harmonic Electromagnetic Fields (McGraw-Hill, 1961).
- A. Taflove, S. Hagness, Computational Electrodynamics, Third Edition (Artech House, 2005).
- A. Taflove, A. Oskooi, S. Johnson, eds., Advances in FDTD Computational Electrodynamics. Photonics and Nanotechnology (Artech House, 2013).
- K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE T. Antenn. Propag. 14, 302–307 (1966). [CrossRef]
- P. Török, P. R. T. Munro, E. E. Kriezis, “High numerical aperture vectorial imaging in coherent optical microscopes,” Opt. Express 16, 507–523 (2008).
- Z. Lin, L. Thylen, “An analytical derivation of the optimum source patterns for the pseudospectral time-domain method,” J. Comput. Phys. 228, 7375–7387 (2009). [CrossRef]
- Q. Li, Y. Chen, C. Li, “Hybrid PSTD-FDTD technique for scattering analysis,” Microw. Opt. Techn. Lett. 34, 19–24 (2002). [CrossRef]
- A. Oskooi, S. G. Johnson, “Electromagnetic wave source conditions,” in Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology, A. Taflove, A. Oskooi, S. G. Johnson, eds. (Artech House, 2013), pp. 65–100.
- V. Čížek, Discrete Fourier Transforms and Their Applications (Adam Hilger, 1986).
- T. Özdenvar, G. A. McMechan, “Causes and reduction of numerical artefacts in pseudo-spectral wavefield extrapolation,” Geophy. J. Int. 126, 819–828 (1996). [CrossRef]
- H. -W. Chen, “Staggered-grid pseudospectral viscoacoustic wave field simulation in two-dimensional media,” J. Acoust. Soc. Am. 100, 120–131 (1996). [CrossRef]
- C. Liu, R. L. Panetta, P. Yang, “Application of the pseudo-spectral time domain method to compute particle single-scattering properties for size parameters up to 200,” J. Quant. Spectrosc. Ra. 113, 1728–1740 (2012). [CrossRef]
- G. J. P. Corrêa, M. Spiegelman, S. Carbotte, J. C. C. Mutter, “Centered and staggered Fourier derivatives and Hilbert transforms,” Geophysics 67, 1558–1563 (2012). [CrossRef]
- B. Richards, E. Wolf, “Electromagnetic diffraction in optical systems ii. Structure of the image field in an aplanatic system,” P. R. Soc. A 253, 358–379 (1959). [CrossRef]
- S. H. Tseng, J. H. Greene, A. Taflove, D. Maitland, V. Backman, J. T. Walsh, “Exact solution of Maxwell’s equations for optical interactions with a macroscopic random medium: addendum,” Opt. Lett. 30, 56–57 (2005). [CrossRef] [PubMed]

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.