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 tangential 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 region 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 frequency space, the dielectric constant and the permeability have to be tensorial 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 discretisation 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