Previously the NL-FDFD method was applied for one-dimensional structures both for depleted and undepleted pump cases. In this section simulation results are presented for a two-dimensional structure: lattices of nonlinear dielectric cylinders. The geometry is adopted from [31,122]. The cylinders form a one-dimensional lattice with a lattice constant a along the x axis, and are surrounded by vacuum. The structure is finite in the y direction, and infinite and periodic in thexdirection with period a (it was denoted by Λ throughout this Thesis, but we apply herea for best comparison to former findings in [31,122]). The nonlinear coefficient is chosen as χ(2) = 2 pm/V, the refractive indices nω =n2ω =√
2, and the radius of the cylinders isδ= 0.1a. The structure is irradiated by an incident plane wave at frequency ω with an angle of incidence θ0 in the x−y plane, the incoming field is polarized in the z direction and its amplitude is E0 = 106 V/m. Two special cases will be examined; first one cylinder oriented along the z axis and infinitely
repeated in the x direction (due to the periodic boundary condition), and second, five cylinders with the same orientation, infinitely repeated in the x direction (see Fig. 4.6).
In the second case the distance between the cylinders along the y axis is equal to the periodicity a, so the cylinders form a square lattice, which is finite in the y direction.
The structure is mostly transparent but strong reflection peaks can be observed for some wavelengths; the positions of the peaks can be influenced by the angle of incidence.
Fig. 4.6. Simulation layout for five cylinders to be infinitely repeated in the xdirection.
One array of cylinders is examined first, where the angle of incidence is chosen to be θ0 = 8.25◦. For a better comparison with the reflectance spectrum in [31], the normalized frequency is introduced: γ = ωa/(2πc) = k0a/(2π) = a/λ0, with a = 1000 nm selected for the repetition period in the simulations. ∆x = ∆y = 5 nm resolution is applied in both directions, PBC in the x direction and UPML boundaries in the y direction, and set spacing regions before the UPML regions. The reflectivity is calculated by numerical integration of the Poynting vector of the back-scattered field in a plane close to theUPML boundary.
The reflectivity spectrum as a function of the normalized frequency can be seen in Fig.
4.7, showing good agreement with Fig. 1 in [31] for both the amplitudes and positions of the peaks. The resonant reflection peak with nearly 100% reflectivity is reported at γ = 0.8714 in [31], using the FDFD method the result is γ = 0.8671, which means less than 0.5% difference, which is most probably caused by the finite approximation of the cylinders in the FDFD simulation. Another, smaller reflection peak can be seen at the doubled frequency. StrongerSHGis expected if the fundamental frequency is close to the resonance frequency. This structure can be considered as a Resonant Waveguide Grating (RWG), which is detailed in Sec.1.1.3. Specifically this geometry was examined there and analytical equations provided the same values for the positions of the resonant peaks as the FDFDmethod (see Fig. 1.5).
NL-FDFD simulation was executed to model SHG in the structure. The reflectivity of the second harmonic wave was calculated using field components of the 2ω mesh. Figs.
4.8a and4.8b show the magnitude of the Ez field at the resonance frequency for the FW and SHW meshes, respectively. The source code of the program resulting these figures is
Fig. 4.7. Reflectivity spectrum of one array of cylinders.
listed in Appendix A as PRG Fig. 4.8. The position of the cylinder is denoted by black circles. The graphs agree well with Fig. 3 in [31]. In Fig. 4.8a, the discontinuity of the field at y=−0.5 µm in the ω mesh is at the boundary of the TF-SF regions.
Fig. 4.8. Magnitude of Ez field components for FW (a) and SHW (b) at the resonant frequency for one array of cylinders.
Fig. 4.9 shows the normalized reflective power of the second harmonic wave defined by
R(2) =IRef2ω/I0. (4.17)
Close to the resonance frequency, an enhanced SHG is obtained, in good agreement with Fig. 2 in [31] showing the same magnitude of generated SHW. The difference of the amplitude can be related to the possible shift of the double resonance condition due to the discrete approximation of the cylinder boundaries, as mentioned above.
Fig. 4.9. Normalized power reflectivity ofSHW for one array of cylinders.
For five arrays of cylinders, a strongerSHGis expected. A double resonance condition, where reflectivity peaks are present at bothωand2ωfrequencies, can be found at an angle of incidenceθ = 7.6◦. Executing the same simulations as in the previous case the obtained spectrum is shown in Fig. 4.10, in good agreement with Fig. 4 in [31]. The plot of the normalized reflectivity power, presented in Fig. 4.11, also compares very well with Fig. 5 in [31]. Finally, perfect concordance is found for the magnitude of the Ez field shown at the resonance frequency in Fig. 4.12 compared with Fig. 6 in [31]. The Matlab program resulting the field distributions in Fig. 4.12 is listed in AppendixA as PRG Fig. 4.12.
Fig. 4.10. Reflectivity spectrum of five arrays of cylinders.
Fig. 4.11. Normalized power reflectivity of SHW for five arrays of cylinders.
Fig. 4.12. Magnitude of Ez field components for FW (a) and SHW (b) at the resonant frequency for five arrays of cylinders.
Chapter 5
Application of the pseudospectral technique for the FDFD method
In Chapter 4 a novel technique, called NonLinear Finite Difference Frequency Do-main (NL-FDFD) method, was presented for the simulation of Second Harmonic Gener-ation (SHG) processes, based on the linear Finite Difference Frequency Domain (FDFD) method. As mentioned in the introductory section in Chapter4, some linear methods have been already generalized forSHGsuch as the Beam Propagation Method (BPM), the Fi-nite Element Method (FEM) and the FiFi-nite Difference Time Domain (FDTD) method.
Each method has its advantage for certain purposes. For instance, simulation of SHGof a light pulse undoubtedly fits best to the FDTD method, while structures with narrow resonant frequency peaks promote the FDFDmethod (see Sec.4.5). Both theFDTDand FDFD methods require much finer grid resolution for SHG simulation than in the linear case. This comes from the fact that nonlinear optical processes are very sensitive to the phase of the interacting fields. Therefore, any numerical error caused by numerical disper-sion has a high impact on the simulation result. The need for more accurate simulation of the propagation emerged in the course of modeling second order nonlinear processes like SHG, sum frequency generation, difference frequency generation [32–34,105] and op-tical parametric amplification [42–44]. The pseudospectral method provides a solution to the problem. The generalization of the FDTD method has already been implemented, where the spatial finite differences are replaced by derivation in the pseudospectral sense, resulting in the PseudoSpectral Time Domain (PSTD) method [70,81,123,124]. Both the FDTD and the PSTDmethods have been applied for second order nonlinear simula-tions [81,115–117,125,126].
The calculation of spatial derivatives with higher accuracy is required in the case of theFDFDmethod as well. In Chapter4the extension of theFDFDmethod for simulation
in the nonlinear FDFD simulations, the need for the application of the pseudospectral method is even higher. There has been a follow-up paper of [127] introducing a method called “pseudospectral frequency domain method” [70]. However, the spatial discretization was performed with Chebyshev polynomials and not on the Yee mesh. This method has been applied to simulate light propagation in fibers [128], and in transverse anisotropic waveguides [129].
The aim of this chapter is to incorporate the pseudospectral method into the FDFD method defined on the Yee mesh, both in linear and nonlinear cases. In Sec. 5.1 the formulation of the PseudoSpectral Frequency Domain (PSFD) method will be carried out. In Sec. 5.2 the PSFDmethod will be examined and compared to theFDFD method by means of the numerical phase velocity and anisotropy. In Sec. 5.3 the PSFD method will be generalized for the nonlinear case in order to simulate SHG processes. Still in this section the NonLinear Finite Difference Frequency Domain (NL-FDFD) and the NonLinear PseudoSpectral Frequency Domain (NL-PSFD) methods will be compared by simulation of anSHGprocess in a homogeneous, nonlinear crystal. In Sec.5.4 the feature of oblique incidence from the FDFD method will be generalized for the PSFD method.
In Sec.5.5 this feature will be applied for the simulation of a tilted Quasi Phase Matched (QPM) grating in two ways: first theQPMgrating will be tilted and the angle of incidence for the Fundamental Wave (FW) will be 0◦, and second, theQPMgrating will be parallel with the simulation area and the FW will have oblique incidence with the same angle as the QPM grating in the previous case. Conversion efficiencies and propagation angles of the Second Harmonic Wave (SHW) obtained from the two cases will be compared to each other and to the analytical results as well.
Scaling properties of the FDFD and the PSFD methods will be analyzed roughly in Sec.5.6, and a new method called NonLinear PseudoSpectral Frequency Domain Sylvester equation based (NL-PSFD-SYLV) method will be introduced. It requires the linear dis-persion has to be homogeneous in at least one direction, but the nonlinear disdis-persion can still be two-dimensionally structured. The advantage of this method is that a very huge simulation area can be handled, and the memory consumption is quite low. In Sec. 5.7 the efficiency of theNL-PSFD-SYLV method will be demonstrated by simulation of a real two-dimensional nonlinear photonic crystal. A centered rectangular lattice is selected, and phase matching is realized in a non-collinear way with reciprocal vectorG01. In each grid point the nonlinear coefficient is inverted inside a circular area. The radius of the circles and thus the filling factor is optimized to achieve the highest possible effective nonlinear
coefficient. Conversion efficiency and angle of propagation will be compared to analytical results.