Abstract
We extended the current density convolution finite-difference time-domain (JEC-FDTD) method to plasma photonic crystals using the Crank–Nicolson difference scheme and derived the one-dimensional JEC-Crank–Nicolson (CN)-FDTD iterative equation of plasma photonic crystals. The method eliminated the Courant–Friedrich–Levy (CFL) stability constraint and became completely unconditional stable form. The incomplete Cholesky conjugate gradient (ICCG) algorithm is proposed to solve the equation with a large sparse matrix in the CN-FDTD method as the ICCG method improves the speed of convergence, enhances stability, and reduces memory consumption. The JEC-CN-FDTD method is applied to study the characteristics of time domain and frequency domain in the plasma photonic crystal objects. The high accuracy and efficiency of the JEC-CN-FDTD method are confirmed by computing the characteristic parameters of plasma photonic crystals under different conditions such as the electric field distribution of electromagnetic wave, reflection coefficients, and transmission coefficients. Simulation study showed that the algorithm performed stably and could reduce memory consumption and facilitate computer programming.
1 Introduction
Among the methods that have been used to study the interaction between plasma and electromagnetic wave, the finite-difference time-domain (FDTD) method [1] is one of the most matured and widely used. The FDTD method applied time and space domain difference in Maxwell’s equations and simulated real propagation of electromagnetic field by calculating electric and magnetic fields alternately. This method is relatively simple, efficient and precise, and particularly suitable for parallel computation [2]. The FDTD method has been widely used in many fields of electromagnetism [1–3]. In recent years, many new FDTD methods have been proposed for the purposes of solving electromagnetics problems in the dispersive media, which include the recursive convolution (RC) method [3], the auxiliary differential equation (ADE) method [4], the piecewise linear recursive convolution (PLRC) method [5], the Z-transformation method [6], the current density convolution (JEC) method [7], the matrix exponential (ME) method [8], the piecewise linear JE recursive convolution (PLJERC) method [9], the discrete time shift operator (SO) method [10], and the Runge–Kutta (RKETD) method [11]. However, the FDTD method is not unconditionally stable, and if the selection of time step does not satisfy the Courant–Friedrich–Levy (CFL) stability conditions, the calculation results will be divergent, which greatly reduces its flexibility and the method can only be applied to few complex objects [12]. In recent years, many unconditionally stable FDTD algorithms including the alternating-direction-implicit finite-difference time-domain (ADI-FDTD) method [13], the Crank–Nicolson finite-difference time-domain (CN-FDTD) method [14], and the locally one-dimensional finite-difference time-domain (LOD-FDTD) method [15] have been proposed. The selection of time step in these methods is determined by the precision of calculation without CFL stability conditions, and as the time step can be increased, the computation time can be reduced correspondingly. Although these three methods can be reliably used for simulation study of the electromagnetic field, when time step in the three algorithms is identical and far greater than that is required to meet the CFL stability condition, the accuracy of CN-FDTD is much higher than that of the other two algorithms [12, 14]. As the CN difference scheme produces a large sparse matrix, development of efficient methods to solve equation with large sparse matrix is the key for wide utilisation of the CN-FDTD method [16]. There are numerous methods for solving large algebraic equations, such as the Jacobi method [17] for solving symmetric positive definite equations, the Gauss–Seidel method [18], the successive over-relaxation (SOR) method [19], the least-squares QR factorisation (LSQR) method [20], the symmetric successive over-relaxation (SSOR) method [21], the conjugate gradient (CG) method [22], and the preconditioned conjugate gradient (PCG) method [23], and among these, the CG method is very effective due to its small storage requirement and simplicity for programming; however, when the condition number of matrix becomes large, the iterative convergence speed of the CG method is usually slow [22]. Compared to the CG method, incomplete Cholesky conjugate gradient (ICCG) converges faster, is numerically more stable, and requires relatively smaller computation memory in solving large sparse matrix [24].
Plasma photonic crystals are an artificial periodic structure composed of plasma and medium or vacuum; since it was first proposed in 2004, numerous studies of plasma photonic crystals have been reported, for example, Hitoshi and Atsushi [25] established the dispersion relation of electromagnetic wave in one-dimensional plasma photonic crystals using the Kronig–Penny model; Shiveshwari and Mahto [26] calculated energy band structure of one-dimensional plasma photonic crystals using the transfer matrix method of layered medium; Liu et al. [27] calculated the transmission and reflection spectrum of one-dimensional plasma photonic crystals using the FDTD method, and they discussed the effect of plasma frequency and other factors on the transmission properties in one-dimensional plasma photonic crystals. Not only do plasma photonic crystals have the characteristics of the band gap of photonic crystals, but they also have the dispersive and dissipative properties of plasma, which the conventional dielectric photonic crystals do not have. Therefore, the properties of plasma photonic crystals can be regulated by changing the parameters of plasma artificially, which has significant impact on its application to engineering [28]. The application of the ICCG to the calculation of the FDTD method has been studied [29–31] throughout the electromagnetic field solving process; first of all, a sparse matrix is derived using the corresponding FDTD method, and then the ICCG method is used to solve the matrix equation for fast solution of the electromagnetic field and less memory consumption. The aims of this study are to incorporate the CN ideas into the JEC-FDTD method and to derive the JEC-CN-FDTD iterative formula of plasma photonic crystals in one-dimensional case so that the algorithm can overcome the limit on the time step size arising from the CFL stability condition and solve large sparse matrix equation in the CN-FDTD method using ICCG, which construct preconditioned matrix using incomplete Cholesky decomposition to reduce the condition number of the coefficient matrix [22], thereby improving the convergence speed and reducing processing time and memory usage.
2 JEC-CN-FDTD Algorithm
In order to simulate a medium with a dielectric constant other than one, which corresponds to free space, we have to add the relative dielectric constant εr to Maxwell’s equations; thus, Maxwell’s equations and related equations for non-magnetised collision plasma are given by
where E is the electric field intensity; H is the magnetic field intensity; J is the polarisation current density; ε0 and μ0 are the vacuum permittivity and permeability, respectively; ν is the plasma collision frequency; and ωp is the plasma angular frequency.
For the time harmonic electromagnetic field, (3) can be written as
where
Equations (4) and (5) are transformed with the inverse Fourier; thus, the convolution relationship between electric field intensity and polarisation current density is
where U(τ) is the unit step function, and if σ(τ) in (6) is replaced by (7), then
Use Yee’s notation to make t=nΔt in the one-dimensional case
where the iterative equation of current density of two-order precision can be obtained [8]:
In the CN-FDTD algorithm, Yee’s rectangular differential grid is used, and placement of E and H components of any grid in space follows the traditional FDTD method [32]. Spatial partial differential equations are based on the central difference scheme, and the time partial differential coefficient on the left side of discrete Maxwell’s equations is still based on the central difference scheme; the main difference between CN-FDTD method and traditional FDTD method is the time discretisation on the right side of Maxwell’s equations [16, 32]:
According to CN’s theory, the one-dimensional JEC-CN-FDTD in cartesian coordinates can be formulated as follows:
Here, by using the CN scheme, Δt can appropriately increase Courant-Friedrich-Levy Number (CFLN) times. The relative dielectric constant εr of the medium is introduced here, making the formula more versatile, and it is easy to build plasma photonic crystals composed of different media. Equation (12) cannot be computed directly, and (13) and (14) are used to eliminate
Ex can be obtained by solving (15) with the ICCG algorithm. Then, Ex will be substituted into (13) and (10), and Hy and Jx are obtained.
3 Solving the JEC-CN-FDTD Algorithm Based on ICCG
The method of solving one-dimensional plasma photonic crystals JEC-CN-FDTD algorithm based on ICCG is as follows.
After the derivation above, there will be a problem of calculation of a large sparse matrix. Its form can be written as follows:
This equation can be abbreviated as eixi−1+gixi+1=hi , where e1=0 and gn =0.
The values of ei , fi , and gi are fixed in the one-dimensional plasma photonic crystals. It meets the requirement of the ICCG method that the matrix is symmetric and positive.
The final question after derivation above can be attributed to solving the following matrix equation essentially:
In this article, the ICCG method is proposed to solve the large sparse matrix equation because it has a noticeable advantage in computation time and memory consumption over the other traditional solution methods. The ICCG method consists of three procedures, i.e., incomplete factorisation of the coefficient matrix, the CG solving, and the one-dimensional compressed storage of the coefficient matrix [24, 33].
3.1 Incomplete Factorisation of the Coefficient Matrix
Theoretically, a positive definite and symmetric coefficient matrix A can lead to stable convergence in the CG method; however, if A is ill-conditioned, the CG algorithm will be unsuitable for numerical calculation due to its slow convergence rate. Meijerink and Van Der Vost [34] proposed a simple method for calculating the incomplete Cholesky decomposition, which can improve the linear system of (17) effectively and accelerate the convergence rate of the CG method greatly.
The incomplete Cholesky decomposition is given by
where C is a lower triangular matrix, which can be obtained from the diagonal matrix D defined as
where ajk are the elements in matrix A, and the terms with ajk =0 are not included in the sum apparently. The matrix D can be computed easily, and the matrix C is given by
where U is a lower triangular matrix with ujj =djj for all j, and ujk =ajk for k<j. The incomplete Cholesky factorisation C is as sparse as the lower triangular of A. Thus, C can be computed rapidly without any additional storage requirement since its non-diagonal elements are the same as those of A [24, 33].
3.2 The Conjugate Gradient Solving
The CG procedure for solving (17) is shown as follows [33]:
Let r0=b−Ax0, p0=r0, then
where α and β are constants, and (ri , ri ) denotes a dot product. A is a symmetric positive definite matrix, x0 is an arbitrary given initial solution of (17), pi is an conjugate vector system with
From the incomplete Cholesky decomposition in (18), (17) can be rewritten as the preconditioned system:
Now apply the CG method to the modified (22), and after a little rearrangement of (21), we can derive the ICCG method, the recurrence sequence of which is given in the following:
Let r0=b−Ax0, p0=(CCT)–1r0, then
In the practical numerical procedure, the ICCG algorithm recurs, while convergence criterion |b–Axk |/|r0|<ε is satisfied, where |r0| is the L2 norm of the residual in the first iteration, and ε is a pre-specified tolerance. The success of the algorithm depends on how well CCT approximates A [24, 33]. If (CCT)–1 is a good approximation to A–1, C–1A(CT)–1 will be an approximate identity matrix, which has better condition than the original system, (17); therefore, the CG method will converge very rapidly when applied to the matrix C–1A(CT)–1.
3.3 One-Dimensional Compressed Storage of the Coefficient Matrix
The ICCG method for solving large sparse matrix adopts one-dimensional compressed storage. In order to describe matrix AN×N , three one-dimensional arrays, AAA, AND, and ANC, are introduced.
AAA(NNC) is to store non-zero elements of the original coefficient matrix, where NNC is the total number of non-zero elements in the matrix.
AND(N) is to store the serial number of the diagonal elements in coefficient matrix, where N is the total number of required unknowns, i.e., the number of the diagonal elements including the non-zero elements.
ANC(NNC) is to store the column number of each element in original coefficient matrix.
Since AN×N is symmetric, we only store the diagonal elements and lower triangular elements based on the rule of the line sequential storage technology; furthermore, AN×N is sparse, so we only need to store its non-zero elements, and this approach can greatly save the storage space of computation [33]. Thus, combined with the one-dimensional compressed storage of the coefficient matrix, the ICCG algorithm has a noticeable advantage of reducing storage and time consumption over other traditional solution methods in solving the electromagnetic problems of plasma photonic crystals.
4 Numerical Validation and Discussions
4.1 Numerical Verification of the Stability
It can be proved theoretically that the JEC-CN-FDTD method is inherently unconditionally stable; here, an simple experiment that simulate a pulse propagate in free space after 100 time steps was performed for a more intuitive verification, the computational domain had 200 cells, and the cell size with Δx=7.5 cm was used. As can be seen from Figure 1, when the conventional FDTD method and the JEC-CN-FDTD method were run with same time step Δt=260 ps exceeding the time step limit for the stable conventional FDTD method that is ΔtFDTDmax=Δx/c=250 ps in this case, the conventional FDTD method quickly becomes unstable (Fig. 1a), whereas the JEC-CN-FDTD method still gives stable results (Fig. 1b).
Simulation of a pulse propagate in free space using (a) the conventional FDTD method and (b) the JEC-CN-FDTD method with same time step Δt=260 ps exceeding the time step limit for the stable conventional FDTD method.
4.2 Application and Validation in Plasma Photonic Crystals
In this section, the JEC-CN-FDTD method based on ICCG will be validated and compared with the conventional FDTD method using three examples of plasma photonic crystals.
The basic structure of plasma photonic crystals in the three examples is the same one composed of alternating N layers of plasmas and (N+1) layers of dielectric or vacuum slabs. Calculation model is shown in Figure 2, where A represents medium or vacuum, B represents plasma, and N is the number of the repeating unit. The thickness of each plasma slab B is L, and the thickness ratio of plasma–medium slab is d=1. Ten cells of perfectly matched layer (PML) absorbing boundaries were used at the terminations of the space to eliminate unwanted reflections [35]. Space between the PML and plasma photonic crystals is free. As the CN scheme is free of CFL stability constraint, the time step can increase CFLN times. In this article, the time steps were set as 6Δt, i.e., CFLN=6, the incident wave used in the simulation is a Gaussian-derivative pulsed plane wave [27, 28].
Schematic diagram of one-dimensional plasma photonic crystals.
Example 1 is the plasma photonic crystals that composed of alternating three layers of plasmas and four layers of vacuum slabs. To clearly display the reflected and transmitted electric field, a large plasma slab with a thickness of L=24 cm, which occupies 160 cells, is chosen in this simulation. The thickness of each vacuum slab is the same as that of the plasma slab, which also occupies 160 cells. The plasma frequency is ωp=2π×3 Grad/s, the plasma collision frequency is νc =3 Grad/s, the mesh width is Δx=1.5 mm, and the time step is Δt=2.5 ps.
Table 1 shows the comparison of storage requirements and run time for the ICCG and direct methods in solving tridiagonal matrix of the example after 1500 time steps, and it can be seen that the ICCG method has more advantages in the computation speed and memory consumption than the direct method when solving the sparse matrix. All the calculated data were obtained on a personal computer with a 2.13 GHz Intel Core2 Duo processor with 2 GBs of RAM utilising the C/C++ languages.
Comparison of storage requirements and run time for the ICCG and direct methods in solving tridiagonal matrix of the example after 1500 time steps.
| Solving method | Calculation time (ms) | Storage requirements |
|---|---|---|
| ICCG method | 126 | 3600B |
| Direct method | 1015 | 7156B |
Figure 3a–c shows electric field amplitude versus position after 1500, 2000, and 2500 time steps, and the traditional FDTD solutions of Z-transform [6] and JEC-CN-FDTD solutions designed in this article will be compared together. JEC-CN-FDTD is still proved accurate when time step is as six times as Δt but has higher computational efficiency due to the adaptation of the CN-FDTD difference scheme and the ICCG method. As is shown in the graphs, the wave traveled from one slab to the other with a reflection at each interface. The time domain process of reflection and transmission of electromagnetic wave through the plasma photonic crystals can be seen easily from the figure.
Electric field versus cell position after (a) 1500, (b) 2000, and (c) 2500 time steps.
Example 2 is the plasma photonic crystals that composed of alternating five layers of plasmas and six layers of vacuum slabs. The plasma frequency is ωp=2π×3 Grad/s, and the plasma collision frequency is νc =10 Grad/s. The mesh width is Δx=7.5 μm, and the time step is Δt=0.1 ps. The thickness of each vacuum slab is the same as that of the plasma slab. Figure 4 is frequency domain reflection coefficient of plasma–vacuum photonic crystals when the thickness of the plasma slab is 0.675 and 1.2 mm, respectively. Comparison of the traditional FDTD solutions of Z-transform [6] versus JEC-CN-FDTD solutions shows that the results of the two algorithms are consistent and the JEC-CN-FDTD algorithm is more effective. As shown in the figure, the plasma photonic crystals mainly have the characteristics of plasma at that condition, even though its reflection spectra was periodic.
The frequency domain reflection coefficient of plasma–vacuum photonic crystal when the thickness of the plasma slab is (a) 0.675 mm and (b) 1.2 mm.
Example 3 is the plasma photonic crystal that composed of alternating five layers of plasmas and six layers of dielectric slabs. The relative permittivity of each dielectric slab used for simulation is εr=5, the plasma collision frequency is νc =10 Grad/s, the mesh width is Δx=0.06 mm, the time step is Δt=0.1 ps, and the value of thickness of each plasma slab is L=4.8 mm, occupying 80 cells. The thickness ratio of plasma and medium slab is d=2, so each dielectric slab occupies 40 cells. When the plasma frequency fp is 3, 12 and 24 GHz, respectively, the frequency domain transmission coefficient was computed by the JEC-CN-FDTD method and traditional Z-transform FDTD method [6] as shown in Figure 5. The curves showed that the JEC-CN-FDTD method is accurate and stable. In addition, significant pass band and stop band are observed in the plasma photonic crystals. When the properties of photonic crystal got enhanced and the plasma frequency was low, the effect of plasma frequency on the band gap structure of plasma photonic crystals can be observed.
The frequency domain transmission coefficient when plasma frequency is (a) 3 GHz, (b)12 GHz, and (c) 24 GHz.
5 Conclusion
In this article, we derived the one-dimensional JEC-CN-FDTD difference iterative formula in plasma photonic crystals based on the JEC-FDTD method of CN difference scheme. The method eliminates the CFL conditions and is of unconditional stability. In order to solve the large sparse matrix generated by the CN difference scheme, ICCG was used in the JEC-CN-FDTD algorithm, which can reduce the computer memory consumption, improve the convergence speed, and reduce computation time. Studying on time domain and frequency domain parameters of plasma photonic crystals showed that plasma photonic crystals can have the characteristics of either plasma or the photonic crystals depending on the setting of parameters. Thus, we can control the property of plasma photonic crystals by adjusting these parameters. Calculation of the electric field distribution, reflection coefficient and transmission coefficient of plasma photonic crystals under different conditions and comparison with the traditional Z-transform FDTD show that the method can increase the time step of FDTD, improve computational efficiency and ensure sufficient accuracy. It is anticipated that when this method can be applied to the analysis of 2D and 3D electromagnetic problems, its advantages will become more obvious.
References
[1] K. S. Yee, IEEE T Antenn. Propag. 14, 302 (1966).10.1109/TAP.1966.1138693Search in Google Scholar
[2] J. L. Young and R. O. Nelson, IEEE Antenn. Propag. M. 43, 61 (2001).10.1109/74.920019Search in Google Scholar
[3] A. Cala’ Lesina, A. Vaccari, and A. Bozzoli, Pr. Electromagn. Res. M. 24, 251 (2012).10.2528/PIERM12041904Search in Google Scholar
[4] M. A. Alsunaidi and A. A. Al-Jabr, IEEE Photonic. Tech. L. 21, 817 (2009).10.1109/LPT.2009.2018638Search in Google Scholar
[5] G. E. Atteia and K. F. A. Hussein, Pr. Electromagn. Res. B. 26, 335 (2010).10.2528/PIERB10083102Search in Google Scholar
[6] V. Nayyeri, M. Soleimani, J. Rashed-Mohassel, and M. Dehmollaian, IEEE T. Antenn. Propag. 59, 2268 (2011).10.1109/TAP.2011.2143677Search in Google Scholar
[7] L. Shao-Bin, M. Jin-Jun, and Y. Nai-Chang, Acta Phys. Sin-Ch Ed. 53, 783 (2004).10.7498/aps.53.783Search in Google Scholar
[8] S. Liu and S. Liu, J. Infrared Millim. Te. 30, 1020 (2009).10.1007/s10762-009-9528-0Search in Google Scholar
[9] S. Liu, N. Yuan, and J. Mo, IEEE Microw. Wirel. Co. 13, 187 (2003).10.1109/LMWC.2003.811663Search in Google Scholar
[10] X. Duan, H. W. Yang, X. Kong, and H. Liu, ETRI J. 31, 387 (2009).10.4218/etrij.09.0108.0698Search in Google Scholar
[11] S. Liu, S. Zhong, and S. Liu, Int. J. Infrared Milli. 29, 323 (2008).10.1007/s10762-008-9327-zSearch in Google Scholar
[12] T. Namiki and K. Ito, IEEE Antennas Prop. 1, 192 (1999).Search in Google Scholar
[13] I. Ahmed and Zhizhang (David) Chen, Int. J. Numer. Model El. 17, 237 (2004).10.1002/jnm.543Search in Google Scholar
[14] G. Sun and C. W. Trueman, IEEE T. Microw. Theory 54, 2275 (2006).10.1109/TMTT.2006.873639Search in Google Scholar
[15] V. E. do Nascimento, B.-H. V. Borges, and F. L. Teixeira, IEEE Microw. Wirel. Co. 16, 398 (2006).10.1109/LMWC.2006.877132Search in Google Scholar
[16] K. Xu, Z. H. Fan, D. Z. Ding, and R. S. Chen, Pr. Electromagn. Res. 102, 381 (2010).10.2528/PIER10020606Search in Google Scholar
[17] Y. Wei, Linear Algebra Appl. 439, 2774 (2013).10.1016/j.laa.2013.07.032Search in Google Scholar
[18] Y. Jiang and L. Zou, Procedia Eng. 15, 1647 (2011).10.1016/j.proeng.2011.08.307Search in Google Scholar
[19] G. Lei, J. Yang, Y. Chen, J. Yang, and J. Li, Numer. Meth. Part D E. 30, 1716 (2014).10.1002/num.21876Search in Google Scholar
[20] Y. Qiu and A. Wang, Appl. Math. Comput. 236, 273 (2014).10.1016/j.amc.2014.02.073Search in Google Scholar
[21] J. H. Yun, J. Comput. Appl. Math. 217, 252 (2008).Search in Google Scholar
[22] D. S. Kershaw, J. Comput. Phys. 26, 43 (1978).10.2307/3350835Search in Google Scholar
[23] D. J. Evans and G.-Y. Lei, BIT 35, 48 (1995).10.1007/BF01732978Search in Google Scholar
[24] X. Wu, Geophys. J. Int. 154, 947 (2003).10.1046/j.1365-246X.2003.02018.xSearch in Google Scholar
[25] H. Hitoshi and M. Atsushi, J. Plasma Fusion Res. 80, 89 (2004).10.1585/jspf.80.89Search in Google Scholar
[26] L. Shiveshwari and P. Mahto, Solid State Commun. 138, 160 (2006).10.1016/j.ssc.2005.11.024Search in Google Scholar
[27] S. Liu, W. Hong, and N. Yuan, Int. J. Infrared Milli. 27, 403 (2006).10.1007/s10762-006-9075-xSearch in Google Scholar
[28] L. Shao-Bin, Z. Chuan-Xi, and Y. Nai-Chang, Acta Phys. Sin-Ch Ed. 54, 2804 (2005).10.7498/aps.54.2804Search in Google Scholar
[29] Y.-J. Gao, H.-W. Yang, R. Weng, Q.-X. Niu, Y.-J. Liu, and G.-B. Wang, Mod. Phys. Lett. B. 29, 1550052 (2015).10.1142/S0217984915500529Search in Google Scholar
[30] Z. Peng, H. W. Yang, R. Weng, Y. Gao, and Z. K. Yang, Comput. Phys. Commun. 185, 2387 (2014).10.1016/j.cpc.2014.05.008Search in Google Scholar
[31] Z. Peng, H. W. Yang, R. Weng, and Z. K. Yang, Optik 125, 5669 (2014).10.1016/j.ijleo.2014.07.033Search in Google Scholar
[32] Y. Yang, R. S. Chen, and Z. Ye, Microw. Opt. Techn. Let. 51, 2378 (2009).10.1002/mop.24613Search in Google Scholar
[33] X. Wu, Y. Xiao, C. Qi, and T. Wang, Geophys. Prospect. 51, 567 (2003).10.1046/j.1365-2478.2003.00392.xSearch in Google Scholar
[34] J. A. Meijerink and H.A. Van Der Vorst, Math. Comput. 31, 148 (1977).10.2307/2005786Search in Google Scholar
[35] J.-P. Berenger, J. Comput. Phys. 114, 185 (1994).10.1006/jcph.1994.1159Search in Google Scholar
©2015 by De Gruyter
