《纺织复合材料》课程参考文献（Computational Materials Science，From Basic Principles to Material Properties）05 Simulation of Active and Nonlinear Photonic Nano-Materials in the Finite-Difference Time-Domain（FDTD）Framework

5 Simulation of Active and Nonlinear Photonic Nano-Materials in the Finite-Difference Time-Domain (FDTD)Framework A.Klaedtke,J.Hamm,and O.Hess Advanced Technology Institute,School of Electronics and Physical Sciences, University of Surrey,Guildford,Surrey,GU2 7XH,UK Abstract.A numerical method is presented that unites three-dimensional finite- differnce time-domain (FDTD)computer simulations of active,nonlinear photonic nano-materials with optical Bloch equations describing their microscopic spatio- temporal dynamics.The constituent equations are derived and the algorithm is dis- cussed.Computationally simulated Rabi oscillations closely correspond to analytic results.Fully three-dimensional simulations reveal the nonlinear spatio-temporal dynamics of high-finesse whispering gallery modes in microdisk lasers. 5.1 Introduction With Moore's famous law still holding at the time of writing of this chapter, processors are getting faster and the main memory capacity of computers is increasing steadily.This opens up the possibility to refine physical models to a detailed level unthinkable several decades ago.Closely linked to this remarkable progress in the development of nano-structuring technologies in micro-and nano-electronics,novel photonic nano-materials and structures can be fabricated today that even carry functionalities on sub-wavelength scales. Theoretical models for and computer simulation of photonic nano-mate- rials often rely on well-known simplifications that profit from the assumption that all optical effects may be described on spatial scales that correspond to several times the optical wavelength.However,in the novel nano-photonic materials the typical structural variation are significantly smaller than the wavelength.In these cases (and even close to the wavelength)these sim- plifications can longer be made any more.Indeed,Maxwell's equations de- scribing the spatio-temporal evolution of the electromagnetic field have to be considered in full detail.Compared with simpler models this enormously in- creases the computational complexity.Ironically,the steadily shrinking scales of novel electronic and photonic nano-structures make on the one hand com- puter simulations of their optical and materials properties more tedious but on the other hand provide the basis for the conception and realization of increasingly more powerful computing platforms on which the simulations run. A.Klaedtke,J.Hamm,and O.Hess,Simulation of Active and Nonlinear Photonic Nano-Mater- ials in the Finite-Difference Time-Domain Framework,Lect.Notes Phys.642,75-101(2004) http://www.springerlink.com/ C Springer-Verlag Berlin Heidelberg 2004

5 Simulation of Active and Nonlinear Photonic Nano-Materials in the Finite-Difference Time-Domain (FDTD) Framework A. Klaedtke, J. Hamm, and O. Hess Advanced Technology Institute, School of Electronics and Physical Sciences, University of Surrey, Guildford, Surrey, GU2 7XH, UK Abstract. A numerical method is presented that unites three-dimensional finite￾differnce time-domain (FDTD) computer simulations of active, nonlinear photonic nano-materials with optical Bloch equations describing their microscopic spatio￾temporal dynamics. The constituent equations are derived and the algorithm is dis￾cussed. Computationally simulated Rabi oscillations closely correspond to analytic results. Fully three-dimensional simulations reveal the nonlinear spatio-temporal dynamics of high-finesse whispering gallery modes in microdisk lasers. 5.1 Introduction With Moore’s famous law still holding at the time of writing of this chapter, processors are getting faster and the main memory capacity of computers is increasing steadily. This opens up the possibility to refine physical models to a detailed level unthinkable several decades ago. Closely linked to this remarkable progress in the development of nano-structuring technologies in micro- and nano-electronics, novel photonic nano-materials and structures can be fabricated today that even carry functionalities on sub-wavelength scales. Theoretical models for and computer simulation of photonic nano-mate￾rials often rely on well-known simplifications that profit from the assumption that all optical effects may be described on spatial scales that correspond to several times the optical wavelength. However, in the novel nano-photonic materials the typical structural variation are significantly smaller than the wavelength. In these cases (and even close to the wavelength) these sim￾plifications can longer be made any more. Indeed, Maxwell’s equations de￾scribing the spatio-temporal evolution of the electromagnetic field have to be considered in full detail. Compared with simpler models this enormously in￾creases the computational complexity. Ironically, the steadily shrinking scales of novel electronic and photonic nano-structures make on the one hand com￾puter simulations of their optical and materials properties more tedious but on the other hand provide the basis for the conception and realization of increasingly more powerful computing platforms on which the simulations run. A. Klaedtke, J. Hamm, and O. Hess, Simulation of Active and Nonlinear Photonic Nano-Mater￾ials in the Finite-Difference Time-Domain Framework, Lect. Notes Phys. 642, 75–101 (2004) http://www.springerlink.com/
c Springer-Verlag Berlin Heidelberg 2004

76 A.Klaedtke,J.Hamm,and O.Hess In this chapter we will first describe a well established combination of algorithms to solve the Maxwell curl equations with boundary conditions and discuss their implementation.The code snippets shown are in a pseudo programming language similar to C or C++.Subroutines are called without the global variables as parameters.In the second part of this report,a method is demonstrated which couples the optical Bloch equations in a special form to the electromagnetic field.The optical Bloch equations describe the quantum behaviour of a system interacting with the electric field on the basis of a dipole interaction Hamiltonian.With the resulting time-domain full vectorial Maxwell Bloch equations,the interplay of light with optically active nano- materials can be computationally investigated. 5.2 Finite-Difference in Time-Domain For the calculations involving the electromagnetic fields we will formulate the Maxwell curl equations in the so called Heaviside-Lorentz form.This unit system is conveniant to equalize the magnitude of the numerical values of the electric and magnetic field allowing higher precision when additions are performed on digital microprocessors using floating point numbers.Note that if we were using the MKSA system of units in contrast,we would introduce an imbalance with the dielectric constant of vacuum co and the permeability of vacuum uo being of significantly different magnitude.Furthermore,we conveniently set the speed of light in vacuum c to 1. Frequently,the Maxwell curl equations can be simplified taking the spe- cific properties of the nano-materials into account.Considering non-magnetic materials results in the equivalence of the magnetic induction and the mag- netic field H (in the Heaviside-Lorentz system of units).Furthermore we assume that the backround dielectric constant er should not vary in time (only be spatially dependent)and of scalar nature.These assumptions lead to the following equations describing the coupled spatio-temporal dynamics of the three fields E (the electric field),D(the electric displacement field) and H (the magnetic field): 7×H(r,)=c-10Dr, (5.1) Ot V×E(r,t)=-c-10H(m,) Ot D(r,t)=er(r)E(r,t)+p(r,t) (5.2) The polarisation density p is the link to the(spatio-temporal)material prop- erties.The implementation of those will be discussed in the second part. In 1966,Yee 5.2 presented a formulation for solving the Maxwell curl equations 5.1 in a discretised way.The algorithm named after him is highly efficient and still rather accurate.To translate the differentials in the Maxwell

76 A. Klaedtke, J. Hamm, and O. Hess In this chapter we will first describe a well established combination of algorithms to solve the Maxwell curl equations with boundary conditions and discuss their implementation. The code snippets shown are in a pseudo programming language similar to C or C++. Subroutines are called without the global variables as parameters. In the second part of this report, a method is demonstrated which couples the optical Bloch equations in a special form to the electromagnetic field. The optical Bloch equations describe the quantum behaviour of a system interacting with the electric field on the basis of a dipole interaction Hamiltonian. With the resulting time-domain full vectorial Maxwell Bloch equations, the interplay of light with optically active nano￾materials can be computationally investigated. 5.2 Finite-Difference in Time-Domain For the calculations involving the electromagnetic fields we will formulate the Maxwell curl equations in the so called Heaviside-Lorentz form. This unit system is conveniant to equalize the magnitude of the numerical values of the electric and magnetic field allowing higher precision when additions are performed on digital microprocessors using floating point numbers. Note that if we were using the MKSA system of units in contrast, we would introduce an imbalance with the dielectric constant of vacuum 0 and the permeability of vacuum µ0 being of significantly different magnitude. Furthermore, we conveniently set the speed of light in vacuum c to 1. Frequently, the Maxwell curl equations can be simplified taking the spe￾cific properties of the nano-materials into account. Considering non-magnetic materials results in the equivalence of the magnetic induction and the mag￾netic field H (in the Heaviside-Lorentz system of units). Furthermore we assume that the backround dielectric constant r should not vary in time (only be spatially dependent) and of scalar nature. These assumptions lead to the following equations describing the coupled spatio-temporal dynamics of the three fields E (the electric field), D (the electric displacement field) and H (the magnetic field): ∇ × H(r, t) = c−1 ∂ D(r, t) ∂t , ∇ × E(r, t) = −c−1 ∂ H(r, t) ∂t (5.1) D(r, t) = r(r) E(r, t) + p(r, t) (5.2) The polarisation density p is the link to the (spatio-temporal) material prop￾erties. The implementation of those will be discussed in the second part. In 1966, Yee [5.2] presented a formulation for solving the Maxwell curl equations 5.1 in a discretised way. The algorithm named after him is highly efficient and still rather accurate. To translate the differentials in the Maxwell

5 Simulation of Active and Nonlinear Photonic Nano-Materials 77 equations,Yee made use of a straightforward approach for the difference ap- proximation of second order accuracy in 6s,exemplified in 5.3.The differen- tial at the position p is approximated by the difference quotient of the two neighbouring values at half-step distances 0f(s) lim f(p+s/2)-fp-6s/2) (5.3) os 68+0 0s 8=p In these(from the point of view of the Maxwell equations "auxiliary")equa- tions for the description of material properties,the use of an averaging rule is necessary,as a field value is not calculated for a certain time step.The following equation will be used in this case: f(e)→1im f(t+t/2)+f(t-t/2) (5.4) 6t→0 2 When carrying out the differencing it is essential to avoid as many averaging procedures as possible to achieve a high numerical stability. The field values in the equations in discrete form will be labeled in a special way to make the appearance short and concise.Therefore we introduce a notation for field values F at certain positions in Cartesian space r,y,z and time t.The discrete field values should be pinned on an evenly spaced grid with a spatial stepping of os and a temporal stepping of ot.Offsets can be changed as necessary. F昭k=F(i6s,j6s,ki,m60）=F(z,2,) (5.5) + ijy X Fig.5.1.The "Yee cube"showing the distribution of the field vector components in one grid cell at space point i,j,k.Shown are the Cartesian components of the electric (E)and the magnetic (H)fields

5 Simulation of Active and Nonlinear Photonic Nano-Materials 77 equations, Yee made use of a straightforward approach for the difference ap￾proximation of second order accuracy in δs, exemplified in 5.3. The differen￾tial at the position p is approximated by the difference quotient of the two neighbouring values at half-step distances ∂ f(s) ∂s      s=p → lim δs→0 f(p + δs/2) − f(p − δs/2) δs . (5.3) In these (from the point of view of the Maxwell equations “auxiliary”) equa￾tions for the description of material properties, the use of an averaging rule is necessary, as a field value is not calculated for a certain time step. The following equation will be used in this case: f(t) → lim δt→0 f(t + δt/2) + f(t − δt/2) 2 (5.4) When carrying out the differencing it is essential to avoid as many averaging procedures as possible to achieve a high numerical stability. The field values in the equations in discrete form will be labeled in a special way to make the appearance short and concise. Therefore we introduce a notation for field values F at certain positions in Cartesian space x, y, z and time t. The discrete field values should be pinned on an evenly spaced grid with a spatial stepping of δs and a temporal stepping of δt. Offsets can be changed as necessary. F| m i,j,k := F(i δs, j δs, k δs, m δt) = F(x, y, z, t) (5.5) j−1 j i+1,j+1 i k+1 k k−1 x y z E H z y H E H Ex z y x i−1 Fig. 5.1. The “Yee cube” showing the distribution of the field vector components in one grid cell at space point i, j, k. Shown are the Cartesian components of the electric (E) and the magnetic (H) fields

78 A.Klaedtke,J.Hamm,and O.Hess void calcH () for (int i =isc;i iec;i++){ for (int j=jsc;j<jec;j++){ for (int k ksc;k kec;k++){ H_x(i,j,k)-=cdtds *E_z(i,j+1,k)-E_z(i,j,k)- E_y(1,j,k+1)+E-y(i,j,k)); H_y(i,j,k)-=cdtds *E_x(i,j,k+1)-E_x(i,j,k)- E_z(1+1,j,k)+E_z(i,j,k)); H_z(i,j,k)-=cdtds *E_y(i+1,j,k)-E_y(i,j,k)- E_x(i,j+1,k)+E_x(1,j,k)); Fig.5.2.Pseudo code realisation of the algorithm calculating the evolution of the magnetic field.The global constants isc,jsc and ksc define the lower bound of the core region of the simulation.iec,jec and kec define the upper bounds of the core region.The global constant cdtds is the factor c.6t/6s.The global grid variables H and E represent the electric and magnetic fields. The result of the translation to difference form of the Maxwell curl equa- tions 5.1 is given in 5.6.For brevity,we only show the equation for the first Cartesian component r of the updating equation for the electric displacement field D.The other five scalar equations that complete the curl equation set are of similar form and can be easily deduced. cot D==D+[ H册专计5k-H:肝-k (5.6) -(,册5+传-，mk-） Translating this into a pseudo programming language is straight forward. The algorithm is composed of two parts.One part is the calculation of the new electric displacement field (Fig.5.3)and the other results in the new magnetic field (Fig.5.2). A very important aspect to consider in practice is the order of storage of the three dimensional arrays in memory.The innermost loop(k in the exam- ples)should belong to the fastest index in the array.This means consecutive iterations for the index of this dimension should access neighbouring posi- tions in memory.The cache architecture of the computer can then preload several array elements from the slower main memory needed in successive calculations all at once.The following iterations can then be executed with fast cache access. An interesting detail of the Yee algorithm is the spatial distribution of the grid points as well as their temporal distribution.These spatial positions are

78 A. Klaedtke, J. Hamm, and O. Hess void calcH () { for (int i = isc; i < iec; i++) { for (int j = jsc; j < jec; j++) { for (int k = ksc; k < kec; k++) { H_x(i,j,k) -= cdtds * ( E_z(i,j+1,k) - E_z(i,j,k) - E_y(i,j,k+1) + E_y(i,j,k) ); H_y(i,j,k) -= cdtds * ( E_x(i,j,k+1) - E_x(i,j,k) - E_z(i+1,j,k) + E_z(i,j,k) ); H_z(i,j,k) -= cdtds * ( E_y(i+1,j,k) - E_y(i,j,k) - E_x(i,j+1,k) + E_x(i,j,k) ); } } } } Fig. 5.2. Pseudo code realisation of the algorithm calculating the evolution of the magnetic field. The global constants isc, jsc and ksc define the lower bound of the core region of the simulation. iec, jec and kec define the upper bounds of the core region. The global constant cdtds is the factor c · δt/δs. The global grid variables H and E represent the electric and magnetic fields. The result of the translation to difference form of the Maxwell curl equa￾tions 5.1 is given in 5.6. For brevity, we only show the equation for the first Cartesian component x of the updating equation for the electric displacement field D. The other five scalar equations that complete the curl equation set are of similar form and can be easily deduced. Dx| m+ 1 2 i+ 1 2 ,j,k = Dx| m− 1 2 i+ 1 2 ,j,k + c δt δs   Hz| m i+ 1 2 ,j+ 1 2 ,k − Hz| m i+ 1 2 ,j− 1 2 ,k −  Hy| m i+ 1 2 ,j,k+ 1 2 − Hy| m i+ 1 2 ,j,k− 1 2   (5.6) Translating this into a pseudo programming language is straight forward. The algorithm is composed of two parts. One part is the calculation of the new electric displacement field (Fig. 5.3) and the other results in the new magnetic field (Fig. 5.2). A very important aspect to consider in practice is the order of storage of the three dimensional arrays in memory. The innermost loop (k in the exam￾ples) should belong to the fastest index in the array. This means consecutive iterations for the index of this dimension should access neighbouring posi￾tions in memory. The cache architecture of the computer can then preload several array elements from the slower main memory needed in successive calculations all at once. The following iterations can then be executed with fast cache access. An interesting detail of the Yee algorithm is the spatial distribution of the grid points as well as their temporal distribution. These spatial positions are

5 Simulation of Active and Nonlinear Photonic Nano-Materials 79 void calcD () for (int i =isc;i iec;i++){ for (int j=jsc;j<jec;j++) for (int k ksc;k kec;k++){ D_x(i,j,k)+=cdtds H_z(i,j,k)-H_z(i,j-1,k)- H_y(1,j,k)+Hy(i,j,k-1)); D_y(i,j,k)+=cdtds (H_x(i,j,k)-H_x(i,j,k-1)- Hz(i,j,k)+Hz(1-1,j,k)); D_z(i,j,k)+=cdtds H_y(i,j,k)-H_y(i-1,j,k)- H_x(1,j,k)+Hx(i,j-1,k)） Fig.5.3.Pseudo code realisation of the algorithm calculating the evolution of the electric displacement field.In addition to the global parameters introduced in calcH (Fig.5.2),the global grid D represents the electric displacement field. usually visualized on the so called Yee cube (Fig.5.1),a three dimensional rectangular grid displaying the positions of the field values. The electric displacement D and the magnetic induction B are at the positions of their appertaining fields E and H,which are shown in the Yee cube.Similarly,the polarisation density p is situated at the position of the electric field and the electric displacement field. The discretisation of the material equation 5.2,relating the electric field E to the electric displacement field D and the polarisation density p represents no further difficulty,as all field values are given at the same positions in space and time.A pseudo code implementation is shown in Fig.5.4. void calcE () for (int i=isc;i<iec;it+){ for (int j=jsc;j<jec;j++){ for (int k ksc;k kec;k++){ E_x(i,j,k)=InvEps_x(i,j,k)D_x(i,j,k) E_y(i,j,k)=InvEps_y(i,j,k)*D_y(i,j,k); E_z(i,j,k)=InvEps_z(i,j,k)*D_z(i,j,k); Fig.5.4.Pseudo code realisation of the algorithm calculating the evolution of the electric field.The global grid InvEps represents the inverse of the er material index. The different Cartesian indices x,y and z account for the spatial offset according to the Yee-cube

5 Simulation of Active and Nonlinear Photonic Nano-Materials 79 void calcD () { for (int i = isc; i < iec; i++) { for (int j = jsc; j < jec; j++) { for (int k = ksc; k < kec; k++) { D_x(i,j,k) += cdtds * ( H_z(i,j,k) - H_z(i,j-1,k) - H_y(i,j,k) + H_y(i,j,k-1) ); D_y(i,j,k) += cdtds * ( H_x(i,j,k) - H_x(i,j,k-1) - H_z(i,j,k) + H_z(i-1,j,k) ); D_z(i,j,k) += cdtds * ( H_y(i,j,k) - H_y(i-1,j,k) - H_x(i,j,k) + H_x(i,j-1,k) ); } } } } Fig. 5.3. Pseudo code realisation of the algorithm calculating the evolution of the electric displacement field. In addition to the global parameters introduced in calcH (Fig. 5.2), the global grid D represents the electric displacement field. usually visualized on the so called Yee cube (Fig. 5.1), a three dimensional rectangular grid displaying the positions of the field values. The electric displacement D and the magnetic induction B are at the positions of their appertaining fields E and H, which are shown in the Yee cube. Similarly, the polarisation density p is situated at the position of the electric field and the electric displacement field. The discretisation of the material equation 5.2, relating the electric field E to the electric displacement field D and the polarisation density p represents no further difficulty, as all field values are given at the same positions in space and time. A pseudo code implementation is shown in Fig. 5.4. void calcE () { for (int i = isc; i < iec; i++) { for (int j = jsc; j < jec; j++) { for (int k = ksc; k < kec; k++) { E_x(i,j,k) = InvEps_x(i,j,k) * D_x(i,j,k); E_y(i,j,k) = InvEps_y(i,j,k) * D_y(i,j,k); E_z(i,j,k) = InvEps_z(i,j,k) * D_z(i,j,k); } } } } Fig. 5.4. Pseudo code realisation of the algorithm calculating the evolution of the electric field. The global grid InvEps represents the inverse of the r material index. The different Cartesian indices x, y and z account for the spatial offset according to the Yee-cube