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

The pseudospectral time-domain (PSTD) method is a computationally efficient numerical method for calculating electromagnetic scattering [1,2]. The method appears to have been used first to solve the acoustic wave equation [3] and has been used widely in fields such as seismic imaging [4] to enable larger scattering volumes to be modelled. In electromagnetic scattering, the PSTD method is closely related to the finite-difference time-domain (FDTD) method with the main difference being how field quantities are spatially sampled. The FDTD method evaluates spatial derivatives using central differences and so field quantities must be sampled at least 10 times per wavelength. The PSTD method evaluates spatial derivatives using the discrete Fourier transform allowing the field sampling to approach two samples per wavelength. The result of this is that, for a problem of dimension D, the PSTD method is between 4D and 8D times more efficient than the FDTD method [1].

Being able to model significantly larger computational volumes is particularly important in biomedical optics applications. This is because sub-wavelength features, i.e., those on the sub-cellular scale, can be embedded in volumes with linear dimensions of thousands of wavelengths, i.e., on the tissue scale. The PSTD method has been used to model phenomena arising from scattering by tissue such as enhanced back scattering [5,6], scattering by non-spherical particles [7], and scattering by ensembles of discrete scatterers [8]. To the best of our knowledge, none of these works employed focused illumination. Adapting existing source conditions to this task would require a focused incident field to be specified on a sub-grid occupying between two [9] and 10 [10] planes of points arranged along the nominal direction of beam propagation. The goal of this paper is to develop a method of introducing focused light into a PSTD simulation via points residing on a single plane in the grid. This leads to a straight forward implementation which is very similar to the total-field scattered-field (TFSF) source condition which has been employed for over three decades by the FDTD community [11].

In the remainder of this paper, we discuss the field equivalence principle before giving a very brief introduction to the FDTD and PSTD methods, including the existing source conditions of both methods. We discuss details of evaluating derivatives using the discrete Fourier transform and the relevance of staggered and collocated grids. We then use these results to derive the new PSTD source condition. We conclude with examples of how the new source condition can be employed and an evaluation of the new source condition which demonstrates its efficacy.

2. Field equivalence

We briefly introduce the field equivalence principle as it is central to our new source condition. We employ a mathematical derivation of this principle by Stratton and Chu [12] which appears to us to be very clear. We repeat only as much of this work here as is necessary for our purposes. More detailed explanations of the field equivalence principle are given at length elsewhere by others [1316]. Stratton and Chu employ the following statement of Maxwell’s equations for time harmonic fields of angular frequency ω (assuming the exp(−iωt) convention):

×EiωμH=J*H=ρ*μ×H+iωεE=JE=ρε,
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:
E(x,y,z)=14πS[iωμ(n×H)G+(n×E)×G+(nE)G]da,
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, da 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:
Js=n×H,Js*=n×E,η=εnE
where Js is an equivalent electric surface current density, Js* is an equivalent magnetic surface current density and η an equivalent surface electric charge density. This result is usually referred to as the field equivalence principle.

 

Fig. 1 Diagrams illustrating the field equivalence principle. (a) A lens focuses light into a region bounded by fictitious surface S = S1S2S3S4. (b) The lens and light source are no longer present but the same field is produced within S due to equivalent sources on S.

Download Full Size | PPT Slide | PDF

3. The finite-difference time-domain method

It is beyond the scope of this paper to give a detailed account of the FDTD method, particularly as there are many good books on the subject [17, 18]. We present some salient features of the method to illustrate the innovation and significance of the source condition we are introducing. For brevity, we consider a two-dimensional vacuum where all field quantities are assumed to be uniform in the 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]:

Hyt=1μ[EzxExzJy*]Ext=1ε[HyzJx]Ezt=1ε[HyxJz],
where ε and μ are the electric permittivity and magnetic permeability of free space, respectively, Jy* is an equivalent source magnetic current density and Jx and Jz are electric source current densities. A set of recurrence relations is obtained from these equations by discretising Ex, Ez and Hy 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]:
Hy(i,k+1;nt+1)=Hy(i,k+1;nt)+ΔtμΔxz[Ez(i+1/2,k+1;nt+1/2)Ez(i1/2,k+1;nt+1/2)+Ex(i,k+1/2;nt+1/2)Ex(i,k+3/2;nt+1/2)ΔxzJy*(i,k+1;nt+1/2)]
Ex(i,k+1/2;nt+1/2)=Ex(i,k+1/2;nt1/2)+ΔtεΔxz[Hy(i,k;nt)Hy(i,k+1;nt)ΔxzJx(i,k+1/2;nt)]
Ez(i1/2,k+1;nt+1/2)=Ez(i1/2,k+1;nt1/2)+ΔtεΔxz[Hy(i,k+1;nt)Hy(i1,k+1;nt)ΔxzJz(i1/2,k+1;nt)],
where Δt and Δxz are the temporal and spatial grid spacings, respectively, which are subject to restrictions. The indexing used here is such that Hy(i, k; nt) refers to the value of Hy at position (iΔxz, kΔxz) at time ntΔt. Note also that the field components are offset by half a grid spacing in both space and time.

One way of introducing a field into an FDTD simulation is to begin with an initial source condition [19], where the incident field is used to initialise the FDTD grid before iterations commence. However, the most widely used source condition is the TFSF source condition which makes use of the source terms Jy*, Jx and Jz, appearing in Eqs. (5)(7) to introduce an incident wave at a fictitious surface at each iteration of the FDTD algorithm [11]. This is demonstrated in the left hand diagram in Fig. 2, which shows how a focused illumination could be introduced into the region bounded by the surface S = S1S2S3S4. Alternatively, if S1 is sufficiently large such that S2 and S3 are pushed to infinity and S4 is also pushed to infinity, the source condition can be implemented on S1 only since S2, S3 and S4 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 S1. 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 S1 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.

 

Fig. 2 Demonstration of how the TFSF source condition introduces a focused field into the “Total field” region bounded by S = S1S2S3S4 (left) and how the focused field can be introduced via a single planar interface S1 (center). Three field sample points from the FDTD grid, in the vicinity of the planar interface (right). PML refers to the perfectly matched layer which absorbs outgoing radiation, allowing unbounded simulations to be performed. The intensity map depicting the incident field represents the electric field at an instant in time.

Download Full Size | PPT Slide | PDF

The implementation of this source condition can be understood by considering the update equation for Hy(i, ks; nt), Eq. (5), in the region of the interface between total and scattered field regions at ksΔxz, as depicted for i = 0 in Fig. 2. If Eq. (5) were evaluated directly without accounting for the fact that one Ex term is a scattered quantity and the other a total quantity, the update equation would be inconsistent. However, since the incident field Ex,i(i, ks + 1/2; nt + 1/2) is known analytically, an additional update equation can be executed to make the original update equation consistent [17], taking the form:

Hy(i,ks;nt)={Hy(i,ks;nt)}(5)+ΔtμΔxzEx,i(i,ks+1/2;nt+1/2).
where we have used the notation of Taflove and Hagness [17] to indicate that Hy(i, ks; nt)(5) is the quantity arising from the evaluation of Eq. (5). This can also be interpreted as including a source term Jy*(i,ks;nt) in Eq. (5) [11]. The important point to note is that the spatial derivative of Ex with respect to z, at z = ksΔxz, 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 Ex.

4. The pseudospectral time-domain method

The PSTD method update equations may be derived in a similar way to the FDTD update equations except that spatial derivatives are evaluated using discrete Fourier transforms rather than central differences [1]. The majority of reported implementations of the PSTD method also differ from the FDTD method by employing a so-called collocated grid in which all field components are sampled on the same computational grid. Recently, however, implementations employing a staggered grid have been reported [5] which we will consider in more detail in Sec. (4.3). The problem of finding a suitable source condition for the PSTD method stems from the spatial derivatives being evaluated globally using the discrete Fourier transform. To understand this further, consider again a two-dimensional TM simulation, as depicted in Fig. 3, where plot (a) shows Hy along a line of constant 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 Ex (Eq. (6)) implicitly calculates ∂Hy/∂z using central differences resulting in an error at z = Δxz/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 Hy update equation. The crucial point is that derivatives are evaluated locally.

 

Fig. 3 (a) An example plot of Hy at an instant in time, along a line of constant x with an interface between total and scattered regions at z = 0. (b) The central differences approximation to ∂Hy/∂z. (c) ∂Hy/∂z evaluated using discrete Fourier transform on a collocated grid. (d) ∂Hy/∂z evaluated using the discrete Fourier transform on a staggered grid. In (b)–(d) the analytically evaluated ∂Hy/∂z is plotted.

Download Full Size | PPT Slide | PDF

The PSTD case is more complex. We show how ∂Hy/∂z 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 ∂Hy/∂z is similar to that calculated using central differences but there are two main differences: (1) there are some artifacts around z = 0 resulting from the sharp transition between total and scattered fields; and (2) the error at z = Δxz/2 is now a function of the field everywhere along this particular line of constant 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

To the best of our knowledge, there is currently no analogue of the TFSF source condition suitable for application in the PSTD method. As explained in the previous section, spatial derivatives are calculated non-locally and, thus, source conditions introduced at a single point or plane are not naturally compatible with the PSTD method [2,9,10,21]. Since most PSTD implementations employ a collocated grid, the jump discontinuity in the field at the interface between the total and scattered field regions introduces ringing in the spatial derivatives, as shown in Fig. 3(c). Remedies for this include the initial source condition [22] in which an incident pulse is initialised within the computational grid prior to beginning time stepping. This has the disadvantage that sample points must be added to the computational domain to accommodate the initial field. Another alternative is to introduce the incident field over a number of cells, thus, smoothing the field discontinuity which causes ringing in the calculated field. The optimal source patterns for sources spanning between two and six cells have been calculated [21]. Other implementations employ sources spanning between eight and ten [10], four and six [2], and two cells [9], respectively. These are very effective solutions yet entail more complex implementation and additional sample points within the computational domain.

4.2. Discrete Fourier transform operators

In order to develop the PSTD update equations, we define the discrete Fourier and inverse Fourier transforms according to:

Xn=k=0N1xkexp(ink2π/N),n=0,1,,N1
xk=1Nn=0N1Xnexp(ink2π/N),k=0,1,,N1,
where xk is a sequence obtained by sampling a function of a continuous variable. We assume N to be even since it results in a more computationally efficient evaluation of the discrete Fourier transform and define an index map, an, as
an={n,0nN2nN,N2<nN1.
Using this index map, we define an operator
Sn(Δ)={exp(ian2πΔ/N)),nN/2cos(an2πΔ/N),n=N/2,
which allows us to form the optimal interpolation scheme, in terms of minimum deviation [24], as:
xk+Δ=1Nn=0N1Sn(Δ)Xnexp(ink2π/N),k=0,1,,N1=IDFT{Sn(Δ)DFT{xk}},
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 ΔΔxz, where Δxz is the spatial sampling period.

4.3. Staggered grids

Although the PSTD method was introduced to the FDTD community by Liu [1], pseudospectral methods were developed before this for applications in acoustic modelling and, no doubt, in many other fields. Some early works on acoustic pseudospectral methods [3, 25, 26] showed that staggered grids offer advantages over collocated grids in terms of numerical accuracy, as has been noted recently [5]. Most early implementations of the PSTD method employed collocated grids [1, 2, 710, 22, 27], which are advantageous because electric and magnetic material boundaries coincide, and which is not the case for staggered grids. Since we do not intend to vary the magnetic properties of materials, this is not a disadvantage for our application. Furthermore, it has been shown that spectrally obtained derivatives possess significantly less ringing if they are evaluated on a staggered grid [3, 25, 26]. We now give a brief explanation of this property and simultaneously derive our formula for spectral differentiation. Equation (13) can be used to derive formulae for differentiation on collocated and staggered grids from first principles. In particular, the derivative sampled on a collocated grid, corresponding to samples xk may be obtained by evaluating:

xk=limΔ0xk+ΔxkΔ
and for a staggered grid as:
xk±1/2=limΔ0xk±1/2+Δxk±1/2Δ
by substituting Eq. (13) into Eqs. (14) and (15) and evaluating the limits. We can evaluate spatial derivatives as:
xk=IDFT{DncDFT{xk}}
xk±1/2=IDFT{Dns±DFT{xk}},
where
Dnc={ian2π/N,nN/20,n=N/2
Dns±={exp(±ianπ/N)ian2π/N,nN/2π,n=N/2,
and the superscripts 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, 26, 28], that derivatives of order one possess substantially less ringing when evaluated on a staggered grid rather than on a collocated grid. To understand this, consider evaluating Eqs. (16) and (17) using discrete convolution, in which case a signal with samples xk is differentiated by convolution with the discrete Fourier transforms of Dnc and Dns±. Figure 4 contains plots of Dnc, Dns+ and their discrete Fourier transforms for the case of N = 32. The plots demonstrate the general result that Dnc possesses a discontinuity at the Nyquist component, n = N/2, whereas Dns+ is everywhere continuous. This is the reason why DFT{Dnc} contains substantial global ringing and the reason why spectrally calculated derivatives of discontinuous signals also possess global ringing. This example demonstrates the general result that first derivatives of discontinuous signals may only be evaluated without global ringing by use of a staggered grid [28].

 

Fig. 4 The values of Dnc (a) and Dns+ (c) for the case N = 32 and their discrete Fourier transforms which are purely real ((b) and (d)). ℜ and ℑ refer to real and imaginary parts respectively.

Download Full Size | PPT Slide | PDF

Equations (15) and (19) allows the PSTD update equations to be derived by evaluating the spatial derivatives in Eqs. (5)(7) using Eq. (17) instead of central differences. The update equations for the two-dimensional system of Eqs. (5)(7) are derived in Eqs. (21)(23) in Sec. (5).

5. A new source condition for the PSTD method

So far we have shown why the TFSF source condition is not an effective source condition for the PSTD method. We have given two examples which demonstrate the general result that spectrally evaluated first derivatives possess significantly less ringing when evaluated on a staggered grid [28]. We have also shown that the field equivalence principle provides a physical interpretation for a jump in the value of an electromagnetic field in terms of equivalent current surface densities and charges. The source condition which we introduce here overcomes the problem of defining a discrete form of the field equivalence principle, by inserting an equivalent magnetic current density at a plane but not an equivalent electric current density. In this way, the globally defined derivative does not pose a problem since its behavior in the vicinity of the transition between total and scattered fields does not need to be explicitly accounted for in any other PSTD source update equations.

To proceed, we return to the field equivalence principle expressed in Eqs. (3) and illustrated in Fig. 1. We wish to simulate the scenario depicted in Fig. 1a), where incident light, (Ei, Hi), is focused into a region of interest, bounded by the surface S. As has been explained, the equivalence principle allows us to generate the same field, (Ei, Hi), within 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 S2, S3 and S4 in Fig. 1a). This may be achieved in practice by making S1 large enough such that the incident field is negligible in the vicinity of S2 and S3 whilst moving S4 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 −, where is the unit vector associated with the z Cartesian coordinate. This gives the equivalent sources on S1, which, without loss of generality, we assume to be the plane z = zi, as:

Js=k^×Hi,Js*=k^×Ei,η=εk^Ei.
It is, however, possible to eliminate all but the Js* term from Eqs. (20). In particular, if a perfect electric conductor (PEC) is brought into contact with S1 from below (i.e., z < zi), the electric current and charge densities must vanish as a result of the infinite conductivity, leaving Js*=k^×Ei as the only non-zero source term. Although we could include a PEC in the PSTD simulation, this is restrictive and may result in other implementation problems, as PEC objects are not well represented in the PSTD method. As an alternative, we employ image theory which enables us to remove the PEC by doubling the magnetic surface current density to give Js*=2k^×Ei=2(Eij^,Eii^,0) [14], where î and ĵ are the Cartesian unit vectors in the x and y directions, respectively. Summarising this result, the electric field in the region z > zi will be equal to Ei if an equivalent magnetic surface current density equal to 2(Ei · ĵ, − Ei · î, 0) flows on S1. Another consequence of this approach is that a field will now propagate back into the scattered field region z < zi. This field is the mirror image of the incident field and may be removed by subtraction at the conclusion of the PSTD simulation [23].

It now remains to derive the discrete PSTD source condition update equations. We begin by deriving the PSTD update equations for the two-dimensional system of Eqs. (5)(7). The first step in deriving the PSTD update equations is to remove the electric source current densities (Jx and Jz) 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 Jy* in Eq. (5) is non-zero only on the plane S1 and has value −Ei · î. The final step is to express the spatial derivatives implicit in Eqs. (5)(7) using spectral derivatives. For example, ∂Hy/∂z is evaluated in Eq. (6) by the terms (Hy(i, k; nt) − Hy(i, k + 1; nt))/Δxz, which may be evaluated spectrally using Eq. (17) as IDFTk{Dns+DFTk{Hy(i,k;nt)}}. Finally, assuming that our source will be introduced at grid index ks, the PSTD update equations corresponding to Eqs. (5)(7) may then be expressed as

Hy(i,k+1;nt+1)=Hy(i,k+1;nt)+ΔtμΔxz[IDFTi{DnsDFTi{Ez(i+1/2,k+1;nt+1/2)}IDFTk{Dns+DFTk{Ex(i,k+1/2;nt+1/2)}δks,k+1ΔxzJy*(i,k+1;nt+1/2)]
Ex(i,k+1/2;nt+1/2)=Ex(i,k+1/2;nt1/2)ΔtεΔxz[IDFTk{Dns+DFTk{Hy(i,k;nt)}}]
Ez(i1/2,k+1;nt+1/2)=Ez(i1/2,k+1;nt1/2)+ΔtεΔxz[IDFTi{DnsDFTi{Hy(i,k+1;nt)}}],
where the subscript of the DFT and IDFT operators denotes which index they operate on and δ is the Kronecker delta function. Division of Jy*=2Eii^ by Δxz is necessary to convert from a surface current density to a current density.

It is now possible to observe that the source update term in Eq. (21) is equivalent to an update term:

Hy(i,ks;nt)={Hy(i,ks;nt)}(5)+ΔtμΔxz2Ei(i,ks;nt+1/2)i^
which is very similar to the equivalent source update term for the TFSF source condition in the FDTD method expressed in Eq. (8). It is straight forward to extend the preceding treatment to the three-dimensional case which reveals that this new source condition employs the same update equations as the TFSF source condition in the FDTD method with only three modifications:
  • The incident electric field, Ei, ordinarily employed by the FDTD method should be multiplied by two.
  • The incident magnetic field ordinarily employed by the FDTD method, Hi, associated with Ei should be set to 0.
  • Ei should be evaluated at a distance Δxz/2 in the negative 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.

The similarity of Eq. (24) with the TFSF condition makes it very easy for existing FDTD codes to be converted to PSTD codes. Further, this source condition may be used in FDTD implementations where it is desired to avoid calculation of the incident magnetic field. This source condition cannot be used on a closed surface, such as 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

We begin with an example of a focused beam introduced into a volume large enough to model low-numerical aperture optical imaging modalities such as optical coherence tomography. Figure 5a) shows how a beam of numerical aperture 0.056 propagates into a computational volume 93μ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]. A perfectly matched layer (PML) of 10 cells surrounding the computational volume was used to simulate propagation in free space. The computational grid has 300×300×420 cells. A λ/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].

 

Fig. 5 a) Plots of an x-polarised beam with wavelength 1325nm focused by a lens with numerical aperture 0.056. The values have been normalized by the analytic in-focus value. The beam was introduced at its waist, just after the PML. (top) Plot of |Ex| in the plane y = 0, the position where the beam was introduced, coincides with the perturbation in the field near the left vertical axis of the plot. (lower left) Plot of |Ex| in the plane z = 1.66μm, just after the beam is introduced. (lower right) Plot of |Ex| in the plane z = 129μm just before the PML. b) Plot of |Ex|, normalized by the in-focus analytic value along a line which passes through the focus of the beam which is at z = 0. The plot demonstrates that ripples in the incident beam are restricted to the vicinity of where the beam is introduced.

Download Full Size | PPT Slide | PDF

The main result of Fig. 5a) is that the field is accurately propagated down the grid without ringing, except in the vicinity of the plane where the source is introduced. Figure 5b) shows a line plot of |Ex| 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.

We evaluated the new source condition quantitatively by simulating the propagation of a tightly focused beam in free space. We consider this example as tightly focused beams are more sensitive to numerical dispersion than weakly focused beams because they possess a wider angular spectrum of plane waves. Using a more tightly focused beam, such as may be employed in confocal microscopy, allowed us to use a smaller computational grid and, thus, each simulation required less time and the errors due to the PML could be minimized by introducing it where erroneously reflected fields were insignificant, compared with other sources of error. An 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]. The focused field was introduced into the PSTD grid at a plane some 2λ 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−17s. 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.

 

Fig. 6 The relative error in Ex evaluated at several planes starting from the focus (z = 0) and for several PSTD simulations each having a different value of Δt. Δtmax is the maximum value of Δt for which the PSTD simulation is stable as given by the Courant stability criterion [17].

Download Full Size | PPT Slide | PDF

Tables Icon

Table 1. Summary of simulation parameters used in each test, where Δ is the grid spacing employed in all directions.

The first comparison we present is the magnitude and error of each field component at the plane located 4λ 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.

 

Fig. 7 Images of the magnitude and relative error of each field component for the analytic and PSTD cases on a plane 4λ beyond the focus. Each field magnitude is normalized by the analytic on-axis and in-focus value of Ex. The error is calculated according to |E{x,y,z}AnalyticE{x,y,z}Analytic|/max{|E{x,y,z}Analytic|}, where the numerator is the maximum in order to avoid spurious results when the field has a small magnitude. Best viewed on screen.

Download Full Size | PPT Slide | PDF

Close inspection of Fig. 7 reveals what appears to be artifacts in both the analytic and PSTD cases that are more severe in the analytic case. These artifacts appear in the images of |Ex| 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.

We proceed now to quantify the accuracy of the field propagated by the PSTD method. We use an error metric for this purpose, which we call relative error, defined as:

ε{x,y,z}(k)=ij|E{x,y,z}PSTD(iΔ,jΔ,kΔ)E{x,y,z}Analytic(iΔ,jΔ,kΔ)|2/ij|E{x,y,z}Analytic(iΔ,jΔ,kΔ)|2
that sums over a transverse plane in the PSTD simulation. The first result we show compares the relative error of each field component as the grid spacing is made larger, approaching λ/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, Ex, begins to become inaccurate when Δ exceeds approximately λ/2.5. The next most significant component, Ez, soon follows whilst Ey maintains its accuracy. This result may be understood by reviewing Maxwell’s equations and, in particular, ∇ × H + iωεE = J, which may be expressed in the time domain as:
Ext=1ε[HzyHyzJx];Eyt=1ε[HxzHzxJy];Ezt=1ε[HyxHxyJz].
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 Hy, Hz and Hx, we expect the derivative with the greatest absolute error to be ∂Hy/∂z, which is why Ex exhibits the greatest sensitivity to grid spacing. The most significant term contributing to the Ez update equation is ∂Hy/∂x, which introduces an inaccuracy because Hy becomes inaccurate by the same mechanism as Ex. The most significant term contributing to Ey is ∂Hz/∂x, which becomes inaccurate due to the dependence that Hz has on ∂Ex/∂y. 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 λ/2.5 is optimal for this particular simulation. Figure 8a) also reveals that, whilst a minimum error in Ex of less than 1% can be obtained, the minimum error which can be obtained in Ey and Ez is less than 2%. We believe that this is the result of there being less energy in these components which biases the error metric.

 

Fig. 8 a) The relative error in each field component, evaluated on a plane 4λ behind the focus as the grid spacing is varied; and b) The relative error in Ex evaluated at several planes starting from the focus (z = 0) and for several PSTD simulations each having a different lateral dimension.

Download Full Size | PPT Slide | PDF

The next test compares how the relative error in Ex (the principal field component) at particular planes along the 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.

The final simulation compared the relative error in Ex as Δ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]. These tests demonstrate that all of the three simulation parameters considered must be optimized to obtain a relative error in all field components of less than 2% after several wavelengths of propagation. However, the order of importance of these parameters is: grid spacing, transverse dimension and time step, since any reduction in the relative error by optimizing a less significant parameter will be limited by the more significant parameter(s).

7. Conclusions

We have introduced a new source condition for the PSTD method which employs a staggered grid. This source condition is implemented at a single plane in the PSTD simulation and is designed specifically for focused illuminations such as may be encountered in biomedical optics. This new source condition allows beams focused by high and low numerical aperture lenses to be efficiently introduced into computational volumes large enough to model image formation for systems employing such lenses, for example, in high-resolution microscopy and in optical coherence tomography.

We have shown, for a beam focused by a high numerical aperture lens, that the source condition is accurate to within 2% throughout the computational space when optimal parameters are employed. We have shown that a grid spacing not exceeding λ/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.

Finally, we note that this source condition will enable easy incorporation of focused beams into PSTD simulations. It may be employed in the FDTD method also and presents the opportunity to easily convert existing FDTD code into PSTD code. In particular, employing a staggered grid and a source condition compatible with the FDTD method means that the only modification required is that of evaluating the spatial derivatives spectrally rather than using central differences.

Acknowledgments

P.M. is supported by a Discovery Early Career Research Award from the Australian Research Council ( DE120101331).

References and links

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]  

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

4. J. Virieux and S. Operto, “An overview of full-waveform inversion in exploration geophysics,” Geophysics 74, WCC127 (2009). [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]  

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]  

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]  

14. C. Balanis, Advanced Engineering Electromagnetics (John Wiley and Sons, 1989).

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. 3, 355–391 (1950). [CrossRef]  

16. R. Harrington, Time-Harmonic Electromagnetic Fields (McGraw-Hill, 1961).

17. A. Taflove and S. Hagness, Computational Electrodynamics, Third Edition (Artech House, 2005).

18. A. Taflove, A. Oskooi, and S. Johnson, eds., Advances in FDTD Computational Electrodynamics. Photonics and Nanotechnology (Artech House, 2013).

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]  

20. 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).

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]  

23. A. Oskooi and S. G. Johnson, “Electromagnetic wave source conditions,” in Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology, A. Taflove, A. Oskooi, and S. G. Johnson, eds. (Artech House, 2013), pp. 65–100.

24. V. Čížek, Discrete Fourier Transforms and Their Applications (Adam Hilger, 1986).

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]  

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]  

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]  

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]  

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]  

References

  • View by:
  • |
  • |
  • |

  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]
  3. D. Kosloff and E. Baysal, “Forward modeling by a Fourier method,” Geophysics 47, 1402–1412 (1982).
    [Crossref]
  4. J. Virieux and S. Operto, “An overview of full-waveform inversion in exploration geophysics,” Geophysics 74, WCC127 (2009).
    [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]
  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]
  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]
  14. C. Balanis, Advanced Engineering Electromagnetics (John Wiley and Sons, 1989).
  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. 3, 355–391 (1950).
    [Crossref]
  16. R. Harrington, Time-Harmonic Electromagnetic Fields (McGraw-Hill, 1961).
  17. A. Taflove and S. Hagness, Computational Electrodynamics, Third Edition (Artech House, 2005).
  18. A. Taflove, A. Oskooi, and S. Johnson, eds., Advances in FDTD Computational Electrodynamics. Photonics and Nanotechnology (Artech House, 2013).
  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]
  20. 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).
  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]
  23. A. Oskooi and S. G. Johnson, “Electromagnetic wave source conditions,” in Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology, A. Taflove, A. Oskooi, and S. G. Johnson, eds. (Artech House, 2013), pp. 65–100.
  24. V. Čížek, Discrete Fourier Transforms and Their Applications (Adam Hilger, 1986).
  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]
  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]
  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]
  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]
  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]

2012 (2)

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]

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]

2010 (1)

2009 (2)

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]

J. Virieux and S. Operto, “An overview of full-waveform inversion in exploration geophysics,” Geophysics 74, WCC127 (2009).
[Crossref]

2008 (2)

2005 (2)

2004 (3)

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]

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]

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]

2002 (1)

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

1999 (1)

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]

1997 (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]

1996 (2)

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]

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

1992 (1)

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]

1982 (1)

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

1980 (1)

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]

1966 (1)

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]

1959 (1)

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]

1950 (1)

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. 3, 355–391 (1950).
[Crossref]

1939 (1)

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

Backman, V.

Balanis, C.

C. Balanis, Advanced Engineering Electromagnetics (John Wiley and Sons, 1989).

Baysal, E.

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

Carbotte, S.

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]

Chen, G.

Chen, H. -W.

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

Chen, K.

Chen, Y.

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

Chu, L. J.

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

Cížek, V.

V. Čížek, Discrete Fourier Transforms and Their Applications (Adam Hilger, 1986).

Corrêa, G. J. P.

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]

Ding, M.

Fisher, R.

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]

Gao, X.

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]

Greene, J. H.

Hagness, S.

A. Taflove and S. Hagness, Computational Electrodynamics, Third Edition (Artech House, 2005).

Hagness, S. C.

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]

Harrington, R.

R. Harrington, Time-Harmonic Electromagnetic Fields (McGraw-Hill, 1961).

Johnson, S. G.

A. Oskooi and S. G. Johnson, “Electromagnetic wave source conditions,” in Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology, A. Taflove, A. Oskooi, and S. G. Johnson, eds. (Artech House, 2013), pp. 65–100.

Kattawar, G.

Kim, Y. L.

Kosloff, D.

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

Kriezis, E. E.

Lee, T. -W.

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]

Levine, H.

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. 3, 355–391 (1950).
[Crossref]

Li, C.

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

Li, Q.

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

Lin, Z.

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]

Liu, C.

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]

Liu, Q.

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]

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

Maitland, D.

McMechan, G. A.

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]

Merewether, D.

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]

Mirotznik, M.

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]

Munro, P. R. T.

Mutter, J. C. C.

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]

Operto, S.

J. Virieux and S. Operto, “An overview of full-waveform inversion in exploration geophysics,” Geophysics 74, WCC127 (2009).
[Crossref]

Oskooi, A.

A. Oskooi and S. G. Johnson, “Electromagnetic wave source conditions,” in Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology, A. Taflove, A. Oskooi, and S. G. Johnson, eds. (Artech House, 2013), pp. 65–100.

Özdenvar, T.

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]

Panetta, R. L.

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]

Petre, P.

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]

Prather, D.

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]

Richards, B.

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]

Sarkar, T.

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]

Schwinger, J.

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. 3, 355–391 (1950).
[Crossref]

Smith, F.

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]

Spiegelman, M.

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]

Stratton, J. A.

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

Taflove, A.

Thylen, L.

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]

Török, P.

Tseng, S. H.

Virieux, J.

J. Virieux and S. Operto, “An overview of full-waveform inversion in exploration geophysics,” Geophysics 74, WCC127 (2009).
[Crossref]

Walsh, J. T.

Wolf, E.

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]

Yang, P.

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]

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]

Yee, K.

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]

Commun. Pur. Appl. Math. (1)

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. 3, 355–391 (1950).
[Crossref]

Geophy. J. Int. (1)

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]

Geophysics (3)

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]

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

J. Virieux and S. Operto, “An overview of full-waveform inversion in exploration geophysics,” Geophysics 74, WCC127 (2009).
[Crossref]

IEEE Antenn. Wirel. Pr. (1)

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]

IEEE T. Antenn. Propag. (3)

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]

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]

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]

IEEE T. Geosci. Remote (1)

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]

IEEE T. Nucl. Sci. (1)

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]

J. Acoust. Soc. Am. (1)

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

J. Comput. Phys. (1)

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]

J. Opt. Soc. Am. A (1)

J. Quant. Spectrosc. Ra. (1)

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]

Microw. Opt. Techn. Lett. (2)

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

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

Opt. Express (3)

Opt. Lett. (2)

P. R. Soc. A (1)

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]

Phys. Rev. (1)

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

Other (6)

R. Harrington, Time-Harmonic Electromagnetic Fields (McGraw-Hill, 1961).

A. Taflove and S. Hagness, Computational Electrodynamics, Third Edition (Artech House, 2005).

A. Taflove, A. Oskooi, and S. Johnson, eds., Advances in FDTD Computational Electrodynamics. Photonics and Nanotechnology (Artech House, 2013).

C. Balanis, Advanced Engineering Electromagnetics (John Wiley and Sons, 1989).

A. Oskooi and S. G. Johnson, “Electromagnetic wave source conditions,” in Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology, A. Taflove, A. Oskooi, and S. G. Johnson, eds. (Artech House, 2013), pp. 65–100.

V. Čížek, Discrete Fourier Transforms and Their Applications (Adam Hilger, 1986).

Cited By

OSA participates in Crossref's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (8)

Fig. 1
Fig. 1

Diagrams illustrating the field equivalence principle. (a) A lens focuses light into a region bounded by fictitious surface S = S1S2S3S4. (b) The lens and light source are no longer present but the same field is produced within S due to equivalent sources on S.

Fig. 2
Fig. 2

Demonstration of how the TFSF source condition introduces a focused field into the “Total field” region bounded by S = S1S2S3S4 (left) and how the focused field can be introduced via a single planar interface S1 (center). Three field sample points from the FDTD grid, in the vicinity of the planar interface (right). PML refers to the perfectly matched layer which absorbs outgoing radiation, allowing unbounded simulations to be performed. The intensity map depicting the incident field represents the electric field at an instant in time.

Fig. 3
Fig. 3

(a) An example plot of Hy at an instant in time, along a line of constant x with an interface between total and scattered regions at z = 0. (b) The central differences approximation to ∂Hy/∂z. (c) ∂Hy/∂z evaluated using discrete Fourier transform on a collocated grid. (d) ∂Hy/∂z evaluated using the discrete Fourier transform on a staggered grid. In (b)–(d) the analytically evaluated ∂Hy/∂z is plotted.

Fig. 4
Fig. 4

The values of D n c (a) and D n s + (c) for the case N = 32 and their discrete Fourier transforms which are purely real ((b) and (d)). ℜ and ℑ refer to real and imaginary parts respectively.

Fig. 5
Fig. 5

a) Plots of an x-polarised beam with wavelength 1325nm focused by a lens with numerical aperture 0.056. The values have been normalized by the analytic in-focus value. The beam was introduced at its waist, just after the PML. (top) Plot of |Ex| in the plane y = 0, the position where the beam was introduced, coincides with the perturbation in the field near the left vertical axis of the plot. (lower left) Plot of |Ex| in the plane z = 1.66μm, just after the beam is introduced. (lower right) Plot of |Ex| in the plane z = 129μm just before the PML. b) Plot of |Ex|, normalized by the in-focus analytic value along a line which passes through the focus of the beam which is at z = 0. The plot demonstrates that ripples in the incident beam are restricted to the vicinity of where the beam is introduced.

Fig. 6
Fig. 6

The relative error in Ex evaluated at several planes starting from the focus (z = 0) and for several PSTD simulations each having a different value of Δt. Δtmax is the maximum value of Δt for which the PSTD simulation is stable as given by the Courant stability criterion [17].

Fig. 7
Fig. 7

Images of the magnitude and relative error of each field component for the analytic and PSTD cases on a plane 4λ beyond the focus. Each field magnitude is normalized by the analytic on-axis and in-focus value of Ex. The error is calculated according to | E { x , y , z } Analytic E { x , y , z } Analytic | / max { | E { x , y , z } Analytic | } , where the numerator is the maximum in order to avoid spurious results when the field has a small magnitude. Best viewed on screen.

Fig. 8
Fig. 8

a) The relative error in each field component, evaluated on a plane 4λ behind the focus as the grid spacing is varied; and b) The relative error in Ex evaluated at several planes starting from the focus (z = 0) and for several PSTD simulations each having a different lateral dimension.

Tables (1)

Tables Icon

Table 1 Summary of simulation parameters used in each test, where Δ is the grid spacing employed in all directions.

Equations (26)

Equations on this page are rendered with MathJax. Learn more.

× E i ω μ H = J * H = ρ * μ × H + i ω ε E = J E = ρ ε ,
E ( x , y , z ) = 1 4 π S [ i ω μ ( n × H ) G + ( n × E ) × G + ( n E ) G ] d a ,
J s = n × H , J s * = n × E , η = ε n E
H y t = 1 μ [ E z x E x z J y * ] E x t = 1 ε [ H y z J x ] E z t = 1 ε [ H y x J z ] ,
H y ( i , k + 1 ; n t + 1 ) = H y ( i , k + 1 ; n t ) + Δ t μ Δ x z [ E z ( i + 1 / 2 , k + 1 ; n t + 1 / 2 ) E z ( i 1 / 2 , k + 1 ; n t + 1 / 2 ) + E x ( i , k + 1 / 2 ; n t + 1 / 2 ) E x ( i , k + 3 / 2 ; n t + 1 / 2 ) Δ x z J y * ( i , k + 1 ; n t + 1 / 2 ) ]
E x ( i , k + 1 / 2 ; n t + 1 / 2 ) = E x ( i , k + 1 / 2 ; n t 1 / 2 ) + Δ t ε Δ x z [ H y ( i , k ; n t ) H y ( i , k + 1 ; n t ) Δ x z J x ( i , k + 1 / 2 ; n t ) ]
E z ( i 1 / 2 , k + 1 ; n t + 1 / 2 ) = E z ( i 1 / 2 , k + 1 ; n t 1 / 2 ) + Δ t ε Δ x z [ H y ( i , k + 1 ; n t ) H y ( i 1 , k + 1 ; n t ) Δ x z J z ( i 1 / 2 , k + 1 ; n t ) ] ,
H y ( i , k s ; n t ) = { H y ( i , k s ; n t ) } ( 5 ) + Δ t μ Δ x z E x , i ( i , k s + 1 / 2 ; n t + 1 / 2 ) .
X n = k = 0 N 1 x k exp ( i n k 2 π / N ) , n = 0 , 1 , , N 1
x k = 1 N n = 0 N 1 X n exp ( i n k 2 π / N ) , k = 0 , 1 , , N 1 ,
a n = { n , 0 n N 2 n N , N 2 < n N 1 .
S n ( Δ ) = { exp ( i a n 2 π Δ / N ) ) , n N / 2 cos ( a n 2 π Δ / N ) , n = N / 2 ,
x k + Δ = 1 N n = 0 N 1 S n ( Δ ) X n exp ( i n k 2 π / N ) , k = 0 , 1 , , N 1 = IDFT { S n ( Δ ) DFT { x k } } ,
x k = lim Δ 0 x k + Δ x k Δ
x k ± 1 / 2 = lim Δ 0 x k ± 1 / 2 + Δ x k ± 1 / 2 Δ
x k = IDFT { D n c DFT { x k } }
x k ± 1 / 2 = IDFT { D n s ± DFT { x k } } ,
D n c = { i a n 2 π / N , n N / 2 0 , n = N / 2
D n s ± = { exp ( ± i a n π / N ) i a n 2 π / N , n N / 2 π , n = N / 2 ,
J s = k ^ × H i , J s * = k ^ × E i , η = ε k ^ E i .
H y ( i , k + 1 ; n t + 1 ) = H y ( i , k + 1 ; n t ) + Δ t μ Δ x z [ IDFT i { D n s DFT i { E z ( i + 1 / 2 , k + 1 ; n t + 1 / 2 ) } IDFT k { D n s + DFT k { E x ( i , k + 1 / 2 ; n t + 1 / 2 ) } δ k s , k + 1 Δ x z J y * ( i , k + 1 ; n t + 1 / 2 ) ]
E x ( i , k + 1 / 2 ; n t + 1 / 2 ) = E x ( i , k + 1 / 2 ; n t 1 / 2 ) Δ t ε Δ x z [ IDFT k { D n s + DFT k { H y ( i , k ; n t ) } } ]
E z ( i 1 / 2 , k + 1 ; n t + 1 / 2 ) = E z ( i 1 / 2 , k + 1 ; n t 1 / 2 ) + Δ t ε Δ x z [ IDFT i { D n s DFT i { H y ( i , k + 1 ; n t ) } } ] ,
H y ( i , k s ; n t ) = { H y ( i , k s ; n t ) } ( 5 ) + Δ t μ Δ x z 2 E i ( i , k s ; n t + 1 / 2 ) i ^
ε { x , y , z } ( k ) = i j | E { x , y , z } P S T D ( i Δ , j Δ , k Δ ) E { x , y , z } Analytic ( i Δ , j Δ , k Δ ) | 2 / i j | E { x , y , z } Analytic ( i Δ , j Δ , k Δ ) | 2
E x t = 1 ε [ H z y H y z J x ] ; E y t = 1 ε [ H x z H z x J y ] ; E z t = 1 ε [ H y x H x y J z ] .

Metrics