《纺织复合材料》课程参考文献（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

80 A.Klaedtke,J.Hamm,and O.Hess We should point out that the use of three different grids for the inverse er field is not mandatory.Only one grid is of real significance.Values of er at the half-step points in space can be calculated by interpolation of neighbouring points as needed. The Yee algorithm is highly efficient as only few and computationally fast operations (i.e.multiplications and sums)have to be processed.And it is still a rather accurate method to calculate the evolution of the electromagnetic fields in spite of using only approximations of second order accuracy in the differences.This fact can be attributed to the type of algorithm,called a half-time stepping or leap-frog type.This means that the algorithm leaps from the calculation of the new electric field,using the old magnetic field,to the calculation of the new magnetic field using the thus calculated magnetic field,and so on. It can be shown that in d spatial dimensions the algorithm requires to satisfy the Courant condition c.6t≤d支.6s (5.7) in order to be stable.Another restriction has to be made for the spatial stepping 6s.The wavelengths in the material with the highest diffraction index,which are to be investigated,should be represented by at least 12 grid points.This is a more strict application of the Nyquist criterion which has proven to be a good rule of thumb. It should be mentioned here that there is a drawback of the method. Energy should be conserved by the basic Maxwell equations in free space. However,it turns out that the Yee algorithm is only energy conserving with respect to averaged values over longer times 5.5.If strict energy conservation is required at all times higher orders of approximation for the differentials have been applied in simulations (see Taflove 5.6). 5.3 Uniaxial Perfectly Matching Layers (UPML) Boundary Conditions As the computational grid for the fields can not be infinite,they have to be terminated somewhere.The equations describing the termination are called boundary conditions.There are boundary conditions representing all kinds of real world systems.Metallic or closed boundary conditions are simple to implement,as the field value is just assumed to be vanishing at the position of the border.Electromagnetic waves are then totally reflected.A pseudo code implementation is given in Fig.5.5.Periodic boundaries represent a certain kind of Bloch condition. In our simulations we use open boundaries to model free space.Energy transmitted away from the object of interest by electromagnetic waves has to be lost.Moreover,there should be no feedback of energy from the borders

80 A. Klaedtke, J. Hamm, and O. Hess We should point out that the use of three different grids for the inverse r field is not mandatory. Only one grid is of real significance. Values of r at the half-step points in space can be calculated by interpolation of neighbouring points as needed. The Yee algorithm is highly efficient as only few and computationally fast operations (i.e. multiplications and sums) have to be processed. And it is still a rather accurate method to calculate the evolution of the electromagnetic fields in spite of using only approximations of second order accuracy in the differences. This fact can be attributed to the type of algorithm, called a half-time stepping or leap-frog type. This means that the algorithm leaps from the calculation of the new electric field, using the old magnetic field, to the calculation of the new magnetic field using the thus calculated magnetic field, and so on. It can be shown that in d spatial dimensions the algorithm requires to satisfy the Courant condition c · δt ≤ d− 1 2 · δs (5.7) in order to be stable. Another restriction has to be made for the spatial stepping δs. The wavelengths in the material with the highest diffraction index, which are to be investigated, should be represented by at least 12 grid points. This is a more strict application of the Nyquist criterion which has proven to be a good rule of thumb. It should be mentioned here that there is a drawback of the method. Energy should be conserved by the basic Maxwell equations in free space. However, it turns out that the Yee algorithm is only energy conserving with respect to averaged values over longer times [5.5]. If strict energy conservation is required at all times higher orders of approximation for the differentials have been applied in simulations (see Taflove [5.6]). 5.3 Uniaxial Perfectly Matching Layers (UPML) Boundary Conditions As the computational grid for the fields can not be infinite, they have to be terminated somewhere. The equations describing the termination are called boundary conditions. There are boundary conditions representing all kinds of real world systems. Metallic or closed boundary conditions are simple to implement, as the field value is just assumed to be vanishing at the position of the border. Electromagnetic waves are then totally reflected. A pseudo code implementation is given in Fig. 5.5. Periodic boundaries represent a certain kind of Bloch condition. In our simulations we use open boundaries to model free space. Energy transmitted away from the object of interest by electromagnetic waves has to be lost. Moreover, there should be no feedback of energy from the borders

5 Simulation of Active and Nonlinear Photonic Nano-Materials 81 void calcPCW () for (int k 0;k kmax;k++){ for (int i=0;i<imax;i++){ Ez(1,0,k)=0.; Ex(1,0,k)=0.; 2 for (int j=0;j<jmax;j++){ E_z(0,j,k)=0.; E_-y(0,j,)=0.; d 2 for (int j=0;j<jmax;j++){ for (int i=0;i<imax;i++){ E_-y(i,j,0)=0.; E_x(1,j,0)=0.; } Fig.5.5.Pseudo code realisation of the metallic boundary conditions.The tangen- tial field values of the electric fields are set to 0.It is only necessary to do this on the lower bounds of the grid dimensions,as only these values are accessed.The global parameters imax,jmax and kmax are the upper bounds of the grid dimensions. of the simulation region.It turned out that the so called Uniaxial Perfectly Matching Layers(U-PML)boundary condition is best suited for us. The basic idea behind PML boundaries is to set up a region of conducting material which absorbs electromagnetic waves without reflection.As this re- gion with the special material properties has to be terminated on the edge of PML Fig.5.6.A wave coming from the left enters the PML region and is dampened.It is reflected on the right by a perfectly conducting wall and passes the PML region again being still further absorbed

5 Simulation of Active and Nonlinear Photonic Nano-Materials 81 void calcPCW () { for (int k = 0; k < kmax; k++) { for (int i = 0; i < imax; i++) { E_z(i,0,k) = 0.; E_x(i,0,k) = 0.; } for (int j = 0; j < jmax; j++) { E_z(0,j,k) = 0.; E_y(0,j,k) = 0.; } } for (int j = 0; j < jmax; j++) { for (int i = 0; i < imax; i++) { E_y(i,j,0) = 0.; E_x(i,j,0) = 0.; } } } Fig. 5.5. Pseudo code realisation of the metallic boundary conditions. The tangen￾tial field values of the electric fields are set to 0. It is only necessary to do this on the lower bounds of the grid dimensions, as only these values are accessed. The global parameters imax, jmax and kmax are the upper bounds of the grid dimensions. of the simulation region. It turned out that the so called Uniaxial Perfectly Matching Layers (U-PML) boundary condition is best suited for us. The basic idea behind PML boundaries is to set up a region of conducting material which absorbs electromagnetic waves without reflection. As this re￾gion with the special material properties has to be terminated on the edge of ￾￾ ￾￾ ￾￾ ￾￾ ￾￾ ￾￾ ￾￾ ￾￾ ￾￾ ￾￾ ￾￾ ￾￾ ￾￾￾￾ ￾￾





























PML Fig. 5.6. A wave coming from the left enters the PML region and is dampened. It is reflected on the right by a perfectly conducting wall and passes the PML region again being still further absorbed

82 A.Klaedtke,J.Hamm,and O.Hess the computational cube it is terminated with a totally reflecting boundary. The workings of such an enclosing PM layer is shown in Fig.5.6.A wave coming from the left enters the PML.On the right side,the wave is totally reflected and crosses the PML again being further reduced in magnitude. To describe such a behaviour in terms of the Maxwell equations in fre- quency space,the dielectric constant and the permeability have to be tenso- rial quantities.It turns out that both can be separated into a scalar and a tensorial part,with the tensorial part s being the same for both. VXH=-wers·E (5.8) 7XE=w山rs·H The ansatz requiring a tensor s that absorbs an incoming wave without re- flecting it leads to 5.9. Oi 11=h- (5.9) w The parameters si-with i being one of x,y or z-can be separated into a real()and an imaginary part (oi/w).The real part degrades evanescent waves.The imaginary part absorbs energy. The conversion of the differential equation 5.8 from frequency space to time-domain requires the introduction of two new fields D and B as shown in 5.10 to avoid the calculation of a convolution integral. 0 0 E B=Hr (5.10) With this substitution,it is possible to separate the real and imaginary parts occurring on the right side of 5.8.The transformation to time-domain is then simply the substitution of -ww by the time derivative. 0 aD 0 VxH= Kz ·D Ot 0 Kx】 0 Ky 0 0 (5.11) 8B 7XE=- Kz ·B 0 8t

82 A. Klaedtke, J. Hamm, and O. Hess the computational cube it is terminated with a totally reflecting boundary. The workings of such an enclosing PM layer is shown in Fig. 5.6. A wave coming from the left enters the PML. On the right side, the wave is totally reflected and crosses the PML again being further reduced in magnitude. To describe such a behaviour in terms of the Maxwell equations in fre￾quency space, the dielectric constant and the permeability have to be tenso￾rial quantities. It turns out that both can be separated into a scalar and a tensorial part, with the tensorial part s being the same for both. ∇ × H = − ıω rs · E ∇ × E = ıω µrs · H (5.8) The ansatz requiring a tensor s that absorbs an incoming wave without re- flecting it leads to 5.9. s =   sysz sx 0 sxsz sy 0 sxsy sz   , si = κi − σi ıω (5.9) The parameters si — with i being one of x, y or z — can be separated into a real (κi) and an imaginary part (σi/ω). The real part degrades evanescent waves. The imaginary part absorbs energy. The conversion of the differential equation 5.8 from frequency space to time-domain requires the introduction of two new fields D and B as shown in 5.10 to avoid the calculation of a convolution integral. D = r   sz sx 0 sx sy 0 sy sz   · E; B = µr   sz sx 0 sx sy 0 sy sz   · H (5.10) With this substitution, it is possible to separate the real and imaginary parts occurring on the right side of 5.8. The transformation to time-domain is then simply the substitution of −ıω by the time derivative. ∇ × H =   κy 0 κz 0 κx   · ∂ D ∂t +   σy 0 σz 0 σx   · D ∇ × E = −   κy 0 κz 0 κx   · ∂ B ∂t −   σy 0 σz 0 σx   · B (5.11)

5 Simulation of Active and Nonlinear Photonic Nano-Materials 83 8D 0 0 D= 0 Ot Kz z 0 8E 0 Ot 0 Oy. (5.12) 0 B 07 ·B= 0 8t 0 0」 0 Ky Fig.5.7.The faces (1),edges (2)and corners(3)a PML has to be separated into. An exemplary assembly of regions for the simulation is shown in Fig.5.7. The inner core is surrounded by three types of PML layers which are finally terminated on the surface of the cube.As only that part of a wave travelling in the direction away from the center of the simulation should be absorbed,the fields components which are perpendicular to the faces of the cube should be effected.The PML has got a finite thickness,therefore edges and corners have to have more than one absorption coefficient to be present.On the faces of the cube,the waves should only be absorbed in the perpendicular direction,along the edges in both transverse directions.And in the corners the waves are to be absorbed in all directions.The conductivities and absorption coefficients oi and Ki have to be set accordingly. Discretisation of 5.11 and 5.12 is straightforward with respect to the dis- cretisation rule(5.3).The results for the x component of the electric field and the displacement field update are shown in 5.13 and 5.14.Other components and field updates are of the same form and can easily be derived

5 Simulation of Active and Nonlinear Photonic Nano-Materials 83   κx 0 κy 0 κz   · ∂ D ∂t +   σx 0 σy 0 σz   · D = r      κz 0 κx 0 κy   · ∂ E ∂t +   σz 0 σx 0 σy   · E      κx 0 κy 0 κz   · ∂ B ∂t +   σx 0 σy 0 σz   · B = µr      κz 0 κx 0 κy   · ∂ H ∂t +   σz 0 σx 0 σy   · H    (5.12) 1 2 3 Fig. 5.7. The faces (1), edges (2) and corners (3) a PML has to be separated into. An exemplary assembly of regions for the simulation is shown in Fig. 5.7. The inner core is surrounded by three types of PML layers which are finally terminated on the surface of the cube. As only that part of a wave travelling in the direction away from the center of the simulation should be absorbed, the fields components which are perpendicular to the faces of the cube should be effected. The PML has got a finite thickness, therefore edges and corners have to have more than one absorption coefficient to be present. On the faces of the cube, the waves should only be absorbed in the perpendicular direction, along the edges in both transverse directions. And in the corners the waves are to be absorbed in all directions. The conductivities and absorption coefficients σi and κi have to be set accordingly. Discretisation of 5.11 and 5.12 is straightforward with respect to the dis￾cretisation rule (5.3). The results for the x component of the electric field and the displacement field update are shown in 5.13 and 5.14. Other components and field updates are of the same form and can easily be derived

84 A.Klaedtke,J.Hamm,and O.Hess D.4=- D.+ St Ky +Oy ot ku+3oy8t H+-H:4-达 6s (5.13) H,4.k+专-H,件生k-生 δs (5.14) The pseudo code realisation of the above described and discretized algorithm is shown in Fig.5.8.We restrict the presentation to the electric field part and leave the programming of the magnetic field calculation to the reader. What is missing in this code is the setup of the coefficient arrays cep(5.15 with p{,y,z})which is presented in Fig.5.9.A similar setup has to be worked out for the magnetic field coefficient arrays. 1 ce,(0,)=n(间+2mn(间6t e,1,)-（,间+间i） 1 (5.15) cep(2,)=Kn(-2n(间t ce,3,)=(园-支)t p()+是on()t The implementation of just one highly absorbing sheath with a high ab- sorption coefficient o results in a large numerical error at the interface to the inner region.This error expresses itself in an artificial reflection.It turns out that the use of several layers with a smoothly increasing absorption coeffi- cient resolves the problem.The same applies for the coefficient K.A potential rising has proven to be a good approach and is shown in 5.16 a()=dmax )=1+(x-(月) (5.16) The parameter d is the thickness in layers of the boundary.The exponent p and the values for omax and Kmax can be freely chosen.Taflove shows in his excellent book on FDTD methods [5.6 that a choice for p between 3 and 4 gives reasonable results

84 A. Klaedtke, J. Hamm, and O. Hess Dx| m+ 1 2 i+ 1 2 ,j,k =κy − 1 2 σy δt κy + 1 2 σy δt Dx| m− 1 2 i+ 1 2 ,j,k + δt κy + 1 2 σy δt · Hz| m i+ 1 2 ,j+ 1 2 ,k − Hz| m i+ 1 2 ,j− 1 2 ,k δs − Hy| m i+ 1 2 ,j,k+ 1 2 − Hy| m i+ 1 2 ,j,k− 1 2 δs , (5.13) Ex| m+ 1 2 i+ 1 2 ,j,k =κz − 1 2 σz δt κz + 1 2 σz δt Ex| m− 1 2 i+ 1 2 ,j,k + 1 r  κz + 1 2 σz δt−1   κx + 1 2 σx δt Dx| m+ 1 2 i+ 1 2 ,j,k −  κx − 1 2 σx δt Dx| m− 1 2 i+ 1 2 ,j,k, (5.14) The pseudo code realisation of the above described and discretized algorithm is shown in Fig. 5.8. We restrict the presentation to the electric field part and leave the programming of the magnetic field calculation to the reader. What is missing in this code is the setup of the coefficient arrays cep (5.15 with p ∈ {x, y, z} ) which is presented in Fig. 5.9. A similar setup has to be worked out for the magnetic field coefficient arrays. cep(0, i) =κp(i) + 1 2 σp(i) δt cep(1, i) =  κp(i) + 1 2 σp(i) δt−1 cep(2, i) =κp(i) − 1 2 σp(i) δt cep(3, i) =κp(i) − 1 2 σp(i) δt κp(i) + 1 2 σp(i) δt (5.15) The implementation of just one highly absorbing sheath with a high ab￾sorption coefficient σ results in a large numerical error at the interface to the inner region. This error expresses itself in an artificial reflection. It turns out that the use of several layers with a smoothly increasing absorption coeffi- cient resolves the problem. The same applies for the coefficient κ. A potential rising has proven to be a good approach and is shown in 5.16. σ(x) = σmax x d p κ(x)=1+(κmax − 1) x d p (5.16) The parameter d is the thickness in layers of the boundary. The exponent p and the values for σmax and κmax can be freely chosen. Taflove shows in his excellent book on FDTD methods [5.6] that a choice for p between 3 and 4 gives reasonable results