SOURCE WAVEFORM FOR FDTD SIMULATION
For FDTD simulation, sources are one of the necessary and important parameters which demand special attention. It is the source that excites electric and magnetic fields as a function of time. In general, waveform is selected depending on the problems under consideration. The wavelength spectrum of the source waveform must include all the wavelengths of interest for the simulation. In this chapter, either modulated Gaussian or Gaussian pulses are used for simulations to investigate the scattering characteristics of periodic structures. The field in the time domain is captured and frequencydomain response can be obtained by the Fourier transform of the captured field.
Gaussian Wave
FDTD simulations should provide the numerical results for all the wavelengths present in the spectrum. For simulation, the cell size of the problem space and the highest wavelength of the source waveform plays an important role. One may set the highest wavelength larger than 20 cell size to obtain a reasonably accurate result (to be discussed later in numerical dispersion section). A Gaussian wave as a function of time can be expressed as
Flere г determines the width of the Gaussian pulse and related with maximum frequency component as [65]
Once г corresponding to the maximum frequency is known, it is possible to construct the Gaussian wave for the FDTD simulation. One thing to be noted here is that for FDTD simulation, field values are initialized as zero. Therefore the source should also be zero at t = 0. For doing this, time shifted Gaussian wave is to be constructed so that the value of the waveform is zero at zero instant of time. A timeshifted Gaussian wave is represented as
Here,f_{0} is the time shift. Considering a negligible value of jc(r)(^^{20}) at zero instant of time, it can be shown that [65,66]
Considering the above equation for constructing a Gaussian waveform, the required waveform with a centre wavelength 1550 nm (frequency 193 THz) is constructed and shown in Figure 11.6. To make the value of the waveform zero at zero instant of time, the required amount time shift is 0.011658 ps also shown in the figure.
Modulated Gaussian Wave
A timeshifted cosine modulated Gaussian wave can be represented as
In FDTD simulation, this type of modulated Gaussian is required to find the spectral response for frequency band centred at w_{c} (w_{c} = 2irf_{t} ). To construct a modulated Gaussian wave, a Gaussian wave with required bandwidth is constructed first and then the modulating function with frequency/,, is multiplied with Gaussian waveform as shown in Equation. 11.45. г can be found from the relation [11.65]
FIGURE 11.7 Modulated Gaussian waveform with Gaussian envelop for FDTD simulation.
Figure 11.7 shows a modulated Gaussian wave along with the Gaussian envelope. Here, the central frequency of the modulating signal is 600 THz. Throughout this chapter, these type of Gaussian and modulated Gaussian waveforms are used as the source for FDTD simulations.
Total Field/Scatter Field Correction of the FDTD Source
FDTD method computes the numerical solution of the wave matter interaction located within the computational domain. For that, plane wave propagating in one direction is to be generated. This was quite challenging because FDTD domain is finite whereas plane wave is infinite in the direction perpendicular to the propagation. Total Field/Scattered Field (TF/SF) method (also known as the Huygens surface method), first reported in [67], and later in [68,69] is used to solve this problem. The concept of TF/SF is based on the linearity of Maxwell’s equations where total electric field and total magnetic field (#„„„,) is considered as
These values result from the wave matter interaction in space. Figure 11.8 shows the 1D FDTD lattice with TF/SF regions. Referring Figure 11.8, Equations 11.39 and 11.40 can then be written in terms of TF/SF zoning as
FIGURE 11.8 ID FDTD lattice with TF/SF regions.
These equations are inconsistent and incorrect because in the same equations, some fields are stored as total filed and some fields are stored as scatter fields in the memory of the computer. However as
Equation 11.48 can be made consistent and correct as
This correction fixes the inconsistency at grid point i_{L}, arises because of TF/SF boundary. Similarly, the inconsistency arising at grid point i_{L} 1/2 (Refer Figure 11.8), while calculating the magnetic field, can be avoided as
Similarly, TF/SF sources can be implemented in 2D FDTD formulation also. Details of these FDTD update equations including the TF/SF correction can be found in [70].
NUMERICAL DISPERSION AND STABILITY
FDTD method basically provides an approximate solution for the fields behaviours for the real physical behaviours of the fields. In this method, derivative of a continuous function is approximated using finite difference scheme as already discussed. It was found there the solution accuracy depends on the step size Ax. The error associated with the nonzero step size Ax is known as numerical dispersion [71]. Though, this numerical dispersion, depends on the wavelength, time step (Д/) and the direction of propagation of the wave also. The effect of numerical dispersion is equivalent to filling the medium with a material whose dielectric constant is different from the dielectric constant of the actual material. As a result, even in a homogeneous free space, the numerical velocity of the wave differs from the free space velocity of the wave. In a source free region. Maxwell’s curl equations for а ТЕМ wave propagating in x direction can be written as
These coupled firstorder differential equations can be decoupled and eliminating H., secondorder differential equation for £_{v} can be written as
The solution of monochromatic sinusoidal traveling wave can be written as
Equations 11.55 and 11.56 can be discretized and it can be shown that the dispersion equation can be written as
Equation 11.57 is the dispersion equation because it relates w and к. It is to be noted here that as Ax and At tend to zero, the wave propagation tends to become dispersion less because then the equation gives the exact value of the wave vector in the limiting case. This time step Дt = Ax/c is the magic time step and the phase velocity as (V_{p}) can be defined as
Considering the wave in traveling in free space, Equation 11.58 can be written as
where a (stability ratio) = c_{0}At/Ax and A_{0} is the free space wavelength. It can be easily seen that if a = 1, the time step becomes the magic time step and numerical dispersion can be avoided. Figure 11.9 shows the plot between the phase velocity and cell size Ax/X_{0}. From the figure, it is clear that as the value of a deviates from one, the introduced error increases, i.e. numerical phase velocity becomes more and more different from the actual phase velocity. From the figure, it is also clear that for a particular value of «(other than one), the error can be reduced, if the grid cell size is reduced. For FDTD simulations, the problem space is discretized into cells and then material properties are introduced. Smaller cell size ensures uniform distribution
FIGURE 11.9 Phase velocity plotted against cell sizes for different values of a showing numerical dispersion.
material properties and thereby uniform distribution of field values. But, as the cell sizes are reduced, the cost of computation increases. As a compromise between the accuracy and computational cost, the rule of thumb is that in no case the cell size should exceed Л/10 at the highest frequency. It was found that for an accuracy of about 1 %, the cell size should be //20 at the highest frequency. In FDTD simulation, cell size Av and time step At are closely related. If v_{max} is the maximum phase velocity then the maximum time step (Д/) permitted is Ax/v_{max}. If At > At_{max}, it can be said that the distance travelled by the wave in time At is more than Ac. It means that the wave will miss the next node in simulation which leads to instability. So the condition
is required to be satisfied for stability of the simulation. This is the Courant FriedrichsLewy (CFL) stability condition. Due to inhomogeneity of the medium, phase velocity may vary from cell to cell. Therefore At is chosen as
It means that « = 1/2 and the wave takes a time span of 2At to travel to the next node. From the following example the concept of numerical dispersion will be clear. Let us consider a wave (modulated Gaussian wave) of central frequency 600 THz, is traveling in free space. Then the free space wavelength (A_{0}) will be 0.5 //m. Table 11.1 shows the information regarding the numerical error introduced for different values of stability factor («) and cell size (Ac). Time step is calculated, maintaining the CFL condition. It can be seen from the table that for a = 1, phase velocity remains unchanged even for different values of cell sizes (Ax/A) and the leading edge of the wave travels a distance of 6//m in 20 fs. It is exactly equal to the distance what light should travel in 20 fs in free space. But as the cell sizes are reduced, time steps are also reduced and the results become more and more accurate. But the reduction in
TABLE 11.1
Numerical Dispersion of a Wave of Central Frequency 600 THz, Traveling in Free Space
Дх/Л 
a = 1 
a = 0.5 

1/5 
1/10 
1/20 
1/5 
1/10 
1/20 

Ax 
0.1 pm 
0.05 pm 
0.025 pm 
0.1 pm 
0.05 pm 
0.025 pm 
Problem space 
50 pm 
25 pm 
12.5 pm 
50 pm 
25 pm 
12.5 pm 
At 
0.33 fs 
0.167 fs 
0.083 fs 
0.167 fs 
0.083 fs 
0.042 fs 
v,Jc 
1 
1 
1 
0.9431 
0.9873 
0.9969 
Change in Phase velocity 
 
 
 
5.69% 
1.27% 
0.31% 
S (in 20 fs) 
6 pm 
6 pm 
6 pm 
5.659 pm 
5.9238 pm 
5.9814 pm 
FIGURE 11.10 Traveling wave showing no numerical dispersion for a = 1 even for different cell sizes.
cell size is associated with increased computational cost. Figure 11.10 shows the FDTD simulation results for the traveling wave discussed above for different values of Дх/Л considering a = 1. From the table, it is found that for Д.Г/Л = 1/5, the numerical phase velocity is 94.31 % of the actual phase velocity and the wave travels a distance of 5.659 /mi instead of 6 //m. This error can be reduced if the cell sizes are reduced. It can also be seen from the table that when the (Д.Ш) value is changed to l/20.the numerical phase velocity becomes 99.69 % of its actual value and the wave travels a distance of 5.9814 /mi which is very close to the actual value. Figure 11.11 shows the numerical error introduced in FDTD simulation for a = 0.5 for the same traveling wave. Results shown in Figures 11.10 and 11.11 are strictly in accordance with the results shown in Figure 11.9.