# A Unified Hamiltonian Solution to Maxwell-Schrödinger Equations for Modeling Electromagnetic Field-Particle Interaction

###### Abstract

A novel unified Hamiltonian approach is proposed to solve Maxwell-Schrödinger equation for modeling the interaction between classical electromagnetic (EM) fields and particles. Based on the Hamiltonian of electromagnetics and quantum mechanics, a unified Maxwell-Schrödinger system is derived by the variational principle. The coupled system is well-posed and symplectic, which ensures energy conserving property during the time evolution. However, due to the disparity of wavelengths of EM waves and that of electron waves, a numerical implementation of the finite-difference time-domain (FDTD) method to the multiscale coupled system is extremely challenging. To overcome this difficulty, a reduced eigenmode expansion technique is first applied to represent the wave function of the particle. Then, a set of ordinary differential equations (ODEs) governing the time evolution of the slowly-varying expansion coefficients are derived to replace the original Schrödinger equation. Finally, Maxwell’s equations represented by the vector potential with a Coulomb gauge, together with the ODEs, are solved self-consistently. For numerical examples, the interaction between EM fields and a particle is investigated for both the closed, open and inhomogeneous electromagnetic systems. The proposed approach not only captures the Rabi oscillation phenomenon in the closed cavity but also captures the effects of radiative decay and shift in the open free space. After comparing with the existing theoretical approximate models, it is found that the approximate models break down in certain cases where a rigorous self-consistent approach is needed. This work is helpful for the EM simulation of emerging nanodevices or next-generation quantum electrodynamic systems.

###### keywords:

Maxwell-Schrödinger equation, Hamiltonian, finite-difference time-domain, reduced eigenmode expansion, Rabi oscillation, and radiative decay.^{†}

^{†}journal: Computer Physics Communications\newpagestyle

main \setheadComputer Physics Communications. DOI: 10.1016/j.cpc.2017.02.006 \headrule

## 1 Introduction

Computational electromagnetics (CEM) Chew2001 () is nowadays indispensable in the modeling and simulation of electromagnetic (EM) effects (radiation, scattering, or propagation) in various electronic devices. The applications include, but not limited to, antenna design for communication, radar, or biomedical systems, electromagnetic compatibility or interference (EMC/EMI), signal integrity analysis, etc. As the working frequency of modern electronic systems approaches terahertz or even near-infrared, and the physical dimensions of objects scale down to micron or even nanometers, many challenges arise in the development of traditional CEM techonologies. One of the challenges is the multiphysics simulation due to the fact that the real problem is simultaneously governed by several different physical laws. For instance, the performance of a modern electronic device is not solely governed by the circuit or EM physics, but also depends on thermal, mechanical, or even quantum effects Shao2012 (); Peng2000 (); Chen2012 (); Wang2010 (). Among all these effects, quantum effect becomes increasingly important due to the rapid development of nanotechnologies. At nanoscales, the interaction of EM field and matter is far richer compared to normal radiowave or microwave applications. Moreover, the effect of EM field on a particle (atom, molecule, or quantum dot, etc) can no longer be modeled by the macroscopically-defined polarization and magnetization. Hence, it is imperative to model the quantum effect of particles directly by solving Schrödinger equation. Meanwhile, the particles, which generate quantum current sources, also affect the external EM fields governed by Maxwell’s equations. The particles absorb and emit electromagnetic waves due to electronic transitions between different energy levels. Based on the above description, the traditional CEM methodology should be modified to integrate the relevant quantum effect into the classical EM simulation. For this purpose, we develop a novel simulation approach for the coupled Maxwell-Schrödinger system in the context of EM field-particle interaction where quantum effect plays an important role.

The solution of the coupled Maxwell-Schrödinger system has attracted much interest in the last decade. In Pierantoni2008 () and Pierantoni2009 (), the joint simulation of electronic/electromagnetic characterization of emerging nanodevices has been conducted, where a 3-D transmission line matrix scheme combined with 1-D FDTD is proposed to model the interaction of external EM field with the carbon nanotubes. In Ahmed2010 (), a hybrid approach is proposed for the simulation of nano-devices, where the conventional 1-D FDTD and a locally one-dimensional (LOD) FDTD are implemented to solve the Schrödinger and Maxwell’s equations, respectively.

In both methods, the procedures of coupling EM and quantum mechanical (QM) equations are essentially the same. The electric and magnetic fields ( and ) are first obtained by solving Maxwell’s equations based on initial conditions. The auxiliary vector (magnetic) and scalar (electric) potentials ( and ) can then be obtained. The two terms are injected into the Schrödinger equation, where the wave function () is calculated, and the quantum current is derived. This current is then re-injected into the Maxwell’s equations as a self-generated source. This self-consistent cycle is repeated until a steady solution is achieved. It is shown lately that by properly choosing the gauge, such cyclic process can be simplified. In the modeling of the interaction between a nanoplate and 2-D EM fields Ohnuki2013 (), it is reported that if a “length gauge” is adopted, the EM fields can be directly inserted into the quantum system to improve the computational efficiency. This practice is also implemented to model the interaction of 1-D EM field with electron confined in various 1-D potentials Takeuchi2014 ().

In this work, we will develop a novel unified Hamiltonian approach for numerically solving the Maxwell-Schrödinger system. The solution to this system is essential to EM-particle interactions in semi-classical framework. Here, the physical model considered is a particle coupled to external EM fields in a closed resonant cavity and in the open free space, which are important to electromagnetic sources (maser, laser, etc) Oxborrow2012 (); Ellis2011 () and cavity quantum electrodynamics Andreani1999 (); Hummer2013 (). Different from the existing works, a canonical and unified Hamiltonian is developed to derive the Hamilton’s equations of the Maxwell-Schrödinger system, which ensures the energy conservation of the entire system. For the EM part, the vector potential (not fields) with a Coulomb gauge is updated by the FDTD method that is suitable for arbitrary 3-D dynamic problems. For the QM part, the reduced eigenmode expansion is adopted instead of discretizing the time-dependent Schrödinger equation for alleviating the multiscale difficulty. The well-known Rabi oscillation for characterizing (stimulated) emission and absorption processes of the particle is calculated, where the schematic configuration and the process of Rabi oscillation are shown in Fig. 1. Moreover, the effects of radiative decay and radiative shift for a particle in free space are also studied. The numerical results are compared with the approximate analytic models Gerry2005 (). It is indicated that in some cases, the analytic model breaks down and a rigorous self-consistent simulation is needed.

The rest of the paper is organized as follows: in Section 2, a unified Hamiltonian is derived for the Maxwell-Schrödinger system. A self-consistent approach with the reduced eigenmode expansion is developed in Section 3. After that, in Section 4, several numerical results from the proposed approach are presented and compared with the approximate theoretical models. Particularly, different aspects (including field intensity, cavity loss, etc) affecting the Rabi oscillation are discussed in details. Finally, a brief conclusion is included and possible future works are suggested in Section 5.

## 2 Hamilton’s Equations of the Electromagnetic and Quantum System

For the sake of completeness, some basic equations are reviewed first. The curl-form Maxwell’s equations are given by Kong1975 ()

(1) |

(2) |

On the other hand, the time-dependent Schrödinger equation governing dynamics of a particle is of the form Miller2008 ()

(3) |

where is the mass of the particle, is the potential energy, is the reduced Plank’s constant, and is the momentum operator

(4) |

Under an EM wave illumination, the equation is modified to

(5) |

where is the electric charge of the particle, and are the magnetic and electric potentials associated with the EM field. The electrostatic potential energy (such as nuclear-electron coulombic potential or effective confining potential) is assumed to be independent of time.

To be consistent with quantum theory, the classical EM formulation in terms of and fields should be replaced with the - formulation Masiello2005 (), Chew2014 (). In the EM field-particle interaction problem, the Coulomb gauge (radiation or transverse gauge) within the nonrelativistic limit is adopted, i.e.

(6) |

The Coulomb gauge splits the EM fields into the transverse (electrodynamic) and longitudinal (electrostatic) parts; and the longitudinal part can be directly inserted into the Schrödinger equation as shown in Eq. (5). Therefore, Coulomb gauge greatly simplifies the quantum optics problem and is commonly adopted Kira2014 (). Next, we define an auxiliary variable

(7) |

After that, the total Hamiltonian of the Maxwell-Schrödinger system can be expressed as

(8) |

where

(9) |

(10) |

In the above and are real valued while is complex valued. In Eq. (8), consists of electric and magnetic energy stored in the EM field, and consists the kinetic and potential energy of the quantum system. By invoking the variational principle, the Hamilton’s equations of the EM and QM parts can be derived as

(11) |

(12) |

(13) |

(14) |

where has the expression of

(15) |

which is essentially the QM current generated from the particle illuminated by the external EM fields. The expression of QM current is the same as that of probability current in the continuity equation Abers2003 (); but it can be derived by the variational principle directly. The generated quantum current will in return perturb and deform the EM field. In matrix notation, the Hamilton’s equations can be written as

(16) |

The above is not a canonical form of Hamilton’s equations. To this end, we apply change of variables and split the wave function into the real and imaginary parts in a normalized way

(17) |

Taking and as new (real-valued) variables, the total system can be re-derived as

(18) |

If we define “generalized coordinates and momenta” as

(19) |

the system is shown to be analogous to the standard Hamilton’s equations Miller2008 ()

(20) |

(21) |

Finally we have

(22) |

where

(23) |

Obviously, this system is well-posed, which means the behavior of solution changes continuously with the initial conditions. Also satisfies the following relation

(24) |

where

(25) |

Hence, is a symplectic matrix; and energy (volume) conservation property can be maintained during the time evolution of the Maxwell-Schrödinger system Feng2010 (); Ren2012 (). In other words, if the initial state in the phase space is

(26) |

and after time it becomes

(27) |

we always have

(28) |

The above exterior-derivative identity suggests that Hamiltonian flows preserve the phase-space volumes Feng2010 (). Physically, it means that total energy of the Maxwell-Schrödinger system is conserved.

In principle, the system in Eq. (22) can be solved by the conventional FDTD method Taflove2005 (), where all the four variables , , , are updated in each time step with initial values and proper boundary conditions. However, the wavelength of the EM wave is much larger than that of the particle wave (). To show this, the profiles and supports of EM wave in a cavity of nm and a particle wave at its ground state are plotted in Fig. 2. There is a distinct mismatch of characteristic length between them. In order to guarantee the discretization accuracy, a much refined grid is required around the position of the particle. This multiscale nature leads to a serious efficiency problem in FDTD simulation. In the following section, we will show that this problem can be mitigated by a reduced eigenmode expansion technique.

## 3 Self-Consistent Solution with Reduced Eigenmode Expansion

Without loss of generality, the particle is assumed to be a 3D isotropic quantum harmonic oscillator with two energy levels for ground and excited states. Quantum harmonic oscillator is one of the most important systems in quantum mechanics. It can be applied to model effective confining potential in atoms, molecules and quantum dots Merkt1991 (). The Schrödinger equation for the 3D isotropic quantum harmonic oscillator is Miller2008 ()

(29) |

where is the angular frequency of the oscillator. The two corresponding eigenstates (eigenmodes) are denoted as and with the eigenenergies and , which are respectively for the ground and excited states. The particle absorbs electromagnetic waves when the electron jumps from ground to excited states. Similarly, the particle emits electromagnetic waves when the electron drops from excited to ground states. According to the reduced eigenmode expansion technique, the time-dependent wave function can be expanded as:

(30) |

where and are the unknown slowly-varying expansion coefficients. The fast-varying terms and describe the time evolution of the eigenstates. As shall be seen, only the difference of and , namely the transition frequency matters. The squared magnitudes denote the probabilities of occupation of the corresponding quantum states, which satisfy the following probability conserving relation

(31) |

By applying the Galerkin test, we have

(32) |

(33) |

where the inner product is defined as

(34) |

According to orthonormality of the eigenmodes

(35) |

and selection rule due to the parity of eigenmodes (the integral of an odd function is equal to zero when it is integrated over the whole of space):

(36) |

we can arrive at the following two ordinary differential equations:

(37) |

(38) |

where the expectation value of quantum current can be expressed as

(39) |

Finally, Eqs. (11), (12), (37), (38) and (39) constitute a complete system and can be solved self-consistently by the FDTD method. The initial values of the expansion coefficients ( and ), and the vector potential or the auxiliary variable are needed to be predefined. To guarantee a unique solution, the boundary conditions for the tangential should also be specified. For a closed resonant cavity, the tangential vanishes at the boundary. For the open free space, however, the convolutional perfectly matched layer (CPML) should be employed to absorb the outgoing waves Gedney2000 ().

When the particle (or two-level system) is illuminated by electromagnetic waves in a cavity, it cyclically absorbs photons and re-emits them by (stimulated) emission, which is called Rabi oscillation. If the resonant cavity is ideal without any material loss and radiation (leaky) loss, the emitted electromagnetic waves will react on the particle and thus Rabi oscillation will be cyclic. Once the expansion coefficients in Eq. (30) are numerically obtained, the population inversion (factor) can be defined as:

(40) |

## 4 Numerical Results

### 4.1 Particle in a resonant cavity

Several numerical results are presented in this section to validate the proposed approach. A nanocavity of dimension nm with a spatial grid of nm is considered as the resonant cavity. For simplicity, the cavity works at the fundamental resonant mode (TE), with an initial condition

(41) |

where

(42) |

A particle resides at the center of this cavity and is initially prepared at a superposition state with and . Due to the influence of the external EM fields, the particle oscillates between its ground state and excited state periodically, according to the semi-classical description of the Rabi oscillation Gerry2005 (). This phenomenon can be theoretically predicted by the Rabi model, under certain approximations (see Appendix). We will show in the following that the theoretical approximations can be reproduced by our numerical results in some cases. While in other cases, the approximate model breaks down and a rigorous numerical solution is needed. Besides, by using the proposed numerical methods, the model can be extended to accommodate more complex boundaries and different loss channels.

The FDTD discretizations of Eqs. (11) and (12) are briefly summarized for clarity,

(43) |

(44) |

where and are defined on the same collocated grids. The time-stepping schemes of Eqs. (37) and (38) are similar and will not be repeated again.

#### 4.1.1 Effect from Field Intensity

We first consider the effect from EM field strength. In tuning case (), where , i.e. the resonant working frequency of the cavity is equal to the transition frequency of the particle (see Appendix). In this case, the Rabi frequency () is proportional to the intensity of field [Eqs. (56)–(57)]. We compare the weak field case with and strong field case with . Rabi oscillations of population inversion are calculated and shown in Fig. 3 and Fig. 4. Here, three methods including the theoretical Rabi model [Eqs. (50)–(51)], rotating wave approximation (RWA) [Eqs. (52)–(53)], and the proposed numerical approach, are adopted. For the weak field case, all the curves agree well with each other. For the strong field case, however, due to the fact that the Rabi frequency is comparable with the EM frequency, the RWA is no longer an appropriate approximation. The fast oscillatory term with a frequency of [Eqs. (50)–(51)] is observable. It is shown that our numerical result captures this high frequency modulation well and agrees with the Rabi model.

#### 4.1.2 Effect from Detuning

Then, we consider the effect from the detuning factor . We fix and modify and then . The results are shown in Fig. 5 and Fig. 6. For a small detuning, there is a slight deviation between theoretical models and our numerical results. As the detuning becomes large, the discrepancy becomes more obvious. The Schrödinger equation with the dipole approximation and “length” gauge [Eqs. (48)–(53)] breaks down in this situation.

#### 4.1.3 Effect from Loss

Next, we investigate the effect from ohmic loss via filling conductive lossy meterial into the cavity. Since loss is introduced into our system, the EM fields will be attenuated and therefore the theoretical Rabi models cannot be adopted. The ohmic loss can be introduced via adding another conduction current to Eq. (12). Under the configuration of and , the results with small loss of and large loss are shown in Fig. 7 and Fig. 8. When the loss is small, the Rabi oscillation can be maintained cyclic. However, as the loss becomes larger, the EM field decays quickly, and the periodic Rabi oscillation can no longer be maintained. To measure the effect of the back coupling of the QM current, the Maxwell’s equations are updated with and without the re-injection of the QM current (but the fore coupling is always considered, as is done in most QM calculations, namely, the Schrödinger equation is solved considering the EM illumination). The feedback of the QM current to the EM system is relatively small when the loss is small, and such feedback becomes more significant when the loss becomes larger. Since a real system is always lossy, the unified coupled solution proposed in this paper will be preferred in real applications.

#### 4.1.4 Effect from Inhomogeneous Electromagnetic Environment

Finally, we investigate the influence of inhomogeneous electromagnetic environment on the population inversion of a particle. The particle is placed at the center of a dielectric sphere embedded in the cavity center. The cavity size is and the sphere has the radius of and relative permittivity of . Except for the sphere, other regions of the cavity is filled with air. Here, the governing equation (11) should be modified as

(45) |

where is the position-dependent relative permittivity in the inhomogeneous electromagnetic system. and is the electric flux. The fundamental mode of the system is calculated by the FDTD method. The local-material approach and the subcell based average-material scheme are adopted Sha2007 () to treat the air-dielectric interface. As seen in Fig. 9, the average-material scheme reduces the staircase error significantly; and the EM energy is confined at the sphere core.

To study the population inversion of the particle, the initial excitation is set to be the fundamental mode of the inhomogeneous cavity and the average-material scheme is employed. For comparisons, the results for the air-filled homogeneous cavity are also given and the initial excitation is set to be fundamental mode. For both cases, the particle is first located at the cavity center, where the EM excitation has the same maximum amplitude. The transition frequencies of the particle are matched to the fundamental eigenfrequencies of the two cavities calculated by the FDTD method, respectively. Additionally, the particle position with an offset by 5 grids from the center is also considered. Due to a faster decay of the EM field away from the sphere center, the population inversion shows more significant lower Rabi frequency at the offset point for the inhomogeneous system (See Fig. 10).

### 4.2 Particle in free space

A particle in the open free space is also investigated in this section. The computational domain for EM simulation occupies grids with a space increment of nm and a time increment of fs. The total time steps are with the CPU time around minutes. Ten CPMLs are employed to absorb radiated waves from the particle. The initial status of the particle is at its ground state with and . The particle is driven by a predefined external field and the induced field generated by the QM current. The external E-field is a cosine-modulated Gaussian pulse

(46) |

with an effective frequency range and . Here, we set V/m, radians per second, fs and fs. The transition frequency is detuned from the EM frequency, i.e. . The external A-field with the Coulomb gauge can be obtained by a direct integration of the E-field, which is given by

(47) |

Considering the radiative dephasing (including the radiative decay and shift) is a very slow process, the dipole moment of the particle is set to be times of the quantum harmonic oscillator for observing the significant change of population inversion. In real applications, one can select molecules with the large dipole moment.

Figure 11 shows the calculated population inversion by the self-consistent solution (with the back coupling of the QM current) and by the non-self-consistent solution (without the back coupling of the QM current). For the self-consistent solutions, the results by the CPML absorbing boundary condition and the perfect electric conductor (PEC) boundary condition are also given for comparisons. After the external pulse is decayed to a negligible value, the Rabi oscillation of population is clearly observed for the PEC case. The radiated wave from the particle is reflected back by the PEC walls and is reabsorbed by the particle. Differently, the self-consistent solution of the population in free space (CPML case) keeps a constant and has a lower value compared to the non-self-consistent solution due to the radiative decay. Figure 12 depicts the dipole moment dependent population of the particle in free space. The self-consistent solution is adopted. A remarkable radiative decay and shift are observed as the dipole moment increases. As a result, the particle moves to its ground state.

## 5 Conclusion

A novel unified Hamiltonian approach is proposed to solve the coupled Maxwell-Schrödinger equation self-consistently. This unified coupled system holds the energy conservation property during its time evolution. To overcome the multiscale issue caused by the distinct wavelength mismatch between the EM wave and electron wave, the reduced eigenmode expansion technique is adopted in the quantum system; and the relevant partial differential equations is cast into ODEs. Represented by the vector potential with a Coulomb gauge, Maxwell’s equations are updated with the incorporation of quantum current obtained from the ODEs. Several physical settings (including field intensity, detuning, and material loss) affecting the Rabi oscillation are investigated and discussed, with comparison to theoretical approximate models. Furthermore, radiative decay and shift are also studied for the particle in free space. In future work, we will consider the EM field-particle interaction in more complex environment where material and structure dependent radiation and ohmic losses will modify the transition dynamics significantly.

## Appendix A Rabi Model and Rotating Wave Approximation

Under dipole approximation and “length” gauge, the Schrödinger equation in an EM environment can be expressed as Gerry2005 ()

(48) |

where the EM field is assumed to be monochromatic and is assumed not to be perturbed by the quantum system (namely no back coupling from the QM system)

(49) |

Assuming a two-level system and applying reduced eigenmode expansion, we obtain a coupled set of ordinary differential equations:

(50) |

(51) |

If we expand in exponentials and drop the term with fast oscillation frequency , we have

(52) |

(53) |

This is called the rotating wave approximation (RWA) and the equations can then be solved analytically if the atom is initially prepared at the ground state Gerry2005 ():

(54) |

(55) |

where is the detuning factor that measures the deviation of the atom transition frequency with the EM frequency, and is the Rabi frequency

(56) |

where

(57) |

## Acknowledgments

This work was supported in part by NSFC 61201002, 61201122, in part by HK GRF 711511, UGC AoE/P¨C04/08, and in part by US AOARD 124082, 134140.

## References

## References

- (1) W. C. Chew, J.-M. Jin, E. Michielssen, J. Song, Fast and Efficient Algorithms in Computational Electromagnetics, Artech House, Boston, 2001.
- (2) Y. Shao, Z. Peng, J.-F. Lee, Proceed. Royal Soc. A 468 (2012) 1652–1675.
- (3) P. Kuo-Peng, N. Sadowski, N. J. Batistela, J. P. A. Bastos, IEEE Trans. Magn. 36 (4) (2000) 1458–1461.
- (4) Y. P. Chen, W. E. I. Sha, W. C. H. Choy, L. Jiang, W. C. Chew, Opt. Exp. 20 (18) (2012) 20210–20221.
- (5) X.-P. Wang, W.-Y. Yin, S. He, IEEE Trans. Elec. Dev., 57 (6) (2010) 1382–1389.
- (6) L. Pierantoni, D. Mencarelli, T. Rozzi, IEEE Trans. Micro. Theory Tech. 56 (3) (2008) 654–662.
- (7) L. Pierantoni, D. Mencarelli, T. Rozzi, IEEE Trans. Micro. Theory Tech. 57 (5) (2009) 1147–1155.
- (8) I. Ahmed, E. H. Khoo, E. Li, R. Mittra, IEEE Antennas Wireless Propag. Lett. 9 (2010) 914–917.
- (9) S. Ohnuki, T. Takeuchi, T. Sako, Y. Ashizawa, K. Nakagawa, M. Tanaka, Int. J. Numer. Model. 26 (2013) 533–544.
- (10) T. Takeuchi, S. Ohnuki, T. Sako, IEEE J. Quan. Elec. 50 (5) (2014) 334–339.
- (11) M. Oxborrow, J. D. Breeze, N. M. Alford, Nat. 488 (2012) 353–357.
- (12) B. Ellis, M. A. Mayer, G. Shambat, T. Sarmiento, J. Harris, E. E. Haller, J. Vuc̆ković, Nat. Photo. 5 (2011) 297–300.
- (13) L. C. Andreani, G. Panzarini, J.-M. Gérard, Phys. Rev. B 60 (19) (1999) 13276–13279.
- (14) T. Hümmer, F. J. García-Vidal, L. Martín-Moreno, D. Zueco, Phys. Rev. B 87 (2013) 115419.
- (15) C. Gerry, P. Knight, Introductory Quantum Optics, Cambridge University Press, New York, 2005.
- (16) J. A. Kong, Theory of Electromagnetic Waves, Wiley Intersci., New York, 1975.
- (17) D. A. B. Miller, Quantum Mechanics for Scientists and Engineers, Cambridge University Press, New York, 2008.
- (18) D. Masiello, E. Deumens, Y. Öhrn, Phys. Rev. A 71 (2005) 032108.
- (19) W. C. Chew, Prog. Electromag. Res. 149 (2014) 69–84.
- (20) M. Kira and S. W. Koch, Semiconductor Quantum Optics, Cambridge University Press, 2012.
- (21) E. S. Abers, Quantum Mechanics, Addison-Wesley, 2003.
- (22) K. Feng, M. Qin, Symplectic Geometric Algorithms for Hamiltonian Systems, Zhejiang Publishing United Group, Hangzhou, 2010.
- (23) X. Ren, Z. Huang, X. Wu, S. Lu, H. Wang, L. Wu, S. Li, Comput. Phys. Comm. 183 (2012) 1192–1200.
- (24) A. Taflove, S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd Ed., Artech House, Boston, 2005.
- (25) U. Merkt, J. Huser, M. Wagner, Phys. Rev. B 43 (1991) 7320.
- (26) J. A. Roden and S. D. Gedney, Microw. Opt. Technol. Lett. 27 (2000) 334–339.
- (27) W. E. I. Sha, Z. X. Huang, X. L. Wu, M. S. Chen, J. Comput. Phys. 225(1) (2007) 33-50.