Chapter 6
1620 nm, and a longpass one with transmission range of 1520-1620 nm and reflection range of 1100-1320 nm. These mirrors can be applied to steer the propagation of bichromatic laser fields, e.g. with wavelengths 1310 and 1550 nm, especially in experiments where the higher thermal noise of traditional multilayer dielectric mirrors is inacceptable.
In Chapter4I have constructed a novel method for the simulation ofSHGby extension of the linear FDFD method to second order nonlinear media, the resulting method is named as NonLinear Finite Difference Frequency Domain (NL-FDFD) method. For the case of 3m symmetry of the underlying crystal lattice (e.g. LiNbO3), the structure of the nonlinear tensor was shown to allow for a two-dimensional simulation (see Sec. 4.1.1).
In order to simulate SHG in a two-dimensional structure, I introduced two Yee meshes for the fundamental and second harmonic fields, respectively. Maxwell’s equations for the Fundamental Wave (FW) and the Second Harmonic Wave (SHW) were discretized on the respective mesh. I have also worked out in Sec. 4.1.2 the general formulation for the three-dimensional case with arbitrary symmetry of the dijk tensor. Undepleted SHG was simulated in a homogeneous nonlinear crystal and the obtained conversion efficiency was shown to compare remarkably well with the analytical value. A similar result was obtained for SHG with pump depletion in a one-dimensional Periodically Poled Lithium Niobate (PPLN) crystal as well. I also carried out true two-dimensional simulations to demonstrate the capabilities of the method in higher dimensions - with the third dimension unstructured - using structures of one- and fivefold cylinder arrays assuming periodic variation in both the linear and the nonlinear susceptibilities. I have shown that the reflection spectra obtained by the NL-FDFD method agree with the theoretical values presented in Sec. 1.1.3 for these geometries. The spatial structures of the Ez field, both for ω and 2ω frequencies and both for one- and fivefold cylinder arrays obtained by the NL-FDFD method (Figs.4.8 and 4.12), exhibit very good agreement with the respective NonLinear Finite Difference Time Domain (NL-FDTD) method results in [31]. Moreover, while in [31] the resonant peaks calculated by a NL-FDTD simulation were found to be smaller than those derived by a semi-analytical method, in our NL-FDFD simulations the peak values and normalized reflectivity powers exhibit very good agreement with the semi-analytical results in [31]. The reason is that the FDFD method is better suited for structures having narrow resonance peaks in their spectrum than the FDTDmethod. So in structures where narrow resonance peaks play a role in the second order interaction, the NL-FDFDmethod seems to be a better choice than theNL-FDTDmethod. An additional benefit of the NL-FDFD method is the much easier implementation of oblique incidence
for the FW, and that of linear dispersion, compared to the NL-FDTDmethod.
In Chapter 5 I have further developed theFDFD method by replacing the first order discrete differences by derivation in pseudospectral sense. The resulting method called PseudoSpectral Frequency Domain (PSFD) method gives a much more accurate approxi-mation of the analytical derivatives compared to the standardFDFDmethod. In order to show this, the relative numerical phase velocity values obtained from simulation results were plotted together with the analytical curve for axial and diagonal directions of the simulation grid in Sec. 5.2. As for the FDFD method I have also implemented the sec-ond harmonic nonlinearity in the PSFD method yielding the NonLinear PseudoSpectral Frequency Domain (NL-PSFD) method. It was then compared to the NL-FDFD method by simulating a homogeneous nonlinear crystal. For comparison I have examined how the numerical coherence length and the maximum conversion efficiency value depend on the sampling density. NL-PSFD exhibits better accuracy in general, similar accuracy can be reached in NL-FDFD method by using a 6-8 times finer grid. Similarly to the FDFD method tilted angle of incidence can be realized in the PSFD method as well, which was executed in Sec. 5.5. The correctness of this procedure was demonstrated by simulating a tilted nonlinear QPM grating in two ways: first, the QPM grating was tilted and the wave vector of theFW was parallel to they axis of the simulation domain (see Fig.5.12), second, FW had oblique incidence and the grating was parallel to the y axis of the sim-ulation domain (see Fig. 5.13). I have calculated the conversion efficiency for both cases and compared them with each other and with the analytical solution as well: as shown in Fig. 5.14 the three curves practically coincide. The angle between the Poynting vector of SHW and the grating boundary is plotted in Fig. 5.15 for both cases together with the expected theoretical value, again, the simulation results provide identical solutions which agree with the theoretical one.
In Sec. 5.6 I have deduced a specific method for large volume SHG simulation when only the second order nonlinearity is structured spatially or the linear permittivity takes the form of ε(x, y) =ε(x) +ε(y). I have started from the NL-PSFDmethod but instead of scaling up the core matrix to obtain a large sparse matrix I kept it in the same size as the simulation area, which resulted in a Sylvester type equation. The method is there-fore calledNL-PSFD-SYLV, and the solution is obtained by iterations using the modified conjugate gradient method. In Sec. 5.7we illustrate the efficiency of theNL-PSFD-SYLV method in large scaleSHGsimulation on a real, two-dimensional nonlinear photonic crys-tal. Nonlinear interaction is shown to occur in a non-collinear way, via the G01 Fourier component of χ(2)(r). The phase matching criterion is calculated for the selected wave-length (1500 nm), determining also the required angle of incidence for the FW and the
10 µm of propagation, the wavefront takes the direction required by the phase matching condition. The conversion efficiency is also calculated from the Poynting vector, and found to be in good agreement with the analytical solution.
Thesis points
The major achievements of my Ph.D. research are summarized in the following thesis points:
[T1] I have determined the structure parameters for a number of basic optical devices derived from a T-shaped Surface Relief Grating (SRG) via numerical optimization.
Aimed first at a mirror, I have optimized the initial grating structure to achieve the widest high reflection plateau in the Near-InfraRed (NIR) range subject to the constraint for the narrowest widths of the bottom and top layers constituting the T-shaped ridges. Second, a bandpass filter was designed; its parameters were determined for a center wavelength of 1550 nm. Finally, two kinds of dichroic beam-splitters were realized: a short- and a longpass one. [P1]
[T2] I have extended the linear Finite Difference Frequency Domain (FDFD) method to second order nonlinear media, showing that Second Harmonic Generation (SHG) can be simulated with the NonLinear Finite Difference Frequency Domain (NL-FDFD) method. I have worked out the general formulation for the two-dimensional case with 3m symmetry of the dijk tensor, and also for the three-dimensional case with arbitrary symmetry of the dijk tensor. I have simulated a homogeneous non-linear crystal and a 1D Periodically Poled Lithium Niobate (PPLN) crystal and the obtained conversion efficiencies were shown to compare remarkably well with the analytical values. I have examined two-dimensional Resonant Waveguide Grating (RWG) structures of one- or fivefold cylinder arrays with periodic variation of both the linear and the nonlinear susceptibilities. I have shown that the reflection spectra, the normalized reflectivity power of the Second Harmonic Wave (SHW), and also the Ez field structure calculated at the resonance frequency, are in good agreement with results reported in [31] using a semi-analytical method.[P2]
[T3] I showed how the pseudospectral method can be applied to the FDFD method re-sulting in a new method called PseudoSpectral Frequency Domain (PSFD) method.
I have worked out the NonLinear PseudoSpectral Frequency Domain (NL-PSFD) method and compared it to theNL-FDFDmethod. I have introduced the
implemen-Phase Matched (QPM) grating in two variants: by normal incidence on a slanted structure and by oblique incidence on the same structure but oriented in the axial direction, and validated the results by direct comparison to each other and to the analytical solution. [P3]
[T4] I have introduced a specific method derived from the NL-PSFD method for large volume simulation, where the linear susceptibility tensor may show special spatial structure or be homogeneous. I have shown that in this method the core matrix, instead of being scaled up to a large sparse matrix (like in the FDFD and PSFD methods), can be kept in the same size as the simulation area. The resulting matrix equation in the obtainedNL-PSFD-SYLVmethod was shown to be of Sylvester type and has been solved by iteration using the modified conjugate gradient method. I have simulated a real, two-dimensional nonlinear photonic crystal of periodically poled cylinders in a centered rectangular lattice, where the nonlinear interaction occurs non-collinearly via the first Fourier component of the second order nonlinear tensor. I have calculated the angle of incidence of the Fundamental Wave (FW) for phase matching as well as the expected angle of the generated SHW. From the simulational result I have determined the direction of propagation of the SHWand also the conversion efficiency, both of which agree with the theoretical expectations.
[P3]
List of publications
Major publications related to thesis points
[P1] T. Szarvas and Zs. Kis. ”Optimization of a T-shaped optical grating for specific applications.” Optical Engineering,55(7): 077103 (2016).
[P2] T. Szarvas and Zs. Kis, ”Numerical simulation of nonlinear second harmonic wave generation by the finite difference frequency domain method.” Journal of the Optical Society of America B, 35(4): 731–740 (2018).
[P3] T. Szarvas and Zs. Kis, ”Application of pseudospectral method to the finite difference frequency domain method.”Journal of the Optical Society of America B,36(2): 333–
345 (2019).
Further publications and conferences
[P4] T. Szarvas and Zs. Kis, ”Simulation of wave propagation in periodically structured dielectrics.” In ”Kvantumelektronika 2014: VII. Szimpózium a hazai kvantumelek-tronikai kutatások eredményeiről,” (2014).
[P5] T. Szarvas and Zs. Kis, ”Optimization of surface relief grating for band filter appli-cation by numerical simulation.” In ”17th International Conference on Transparent Optical Networks (ICTON),” (2015).
[P6] T. Szarvas and Zs. Kis, ”Optical simulation of monolithic grating structure.” In
”Mesterpróba 2015” (2015).
[P7] T. Szarvas and Zs. Kis, ”Simulation of second harmonic wave generation by extended Finite Difference Frequency Domain method.” In ”International OSA Network of Students (IONS) Balvanyos 2017,” (2017).
(2018).
[P9] T. Szarvas and Zs. Kis, ”Másodharmonikus keltés numerikus modellezése az FDFD módszer kiterjesztésével.” In ”Kvantumelektronika 2018: Szimpózium a hazai kvan-tumelektronikai kutatások eredményeiről,” (2018).
[P10] T. Szarvas and Zs. Kis, ”Simulation of nonlinear photonic crystals by modified finite difference frequency domain method.” In ”Frontiers in Optics / Laser Science,”
(2018).
Appendices
Appendix A: List of source codes
PRG Fig. 1.5 RWG.m
PRG Fig. 3.2 T_SRG_RCWA_R_lambda_theta.m
PRG Fig. 3.4 T_SRG_FDFD_FieldDist_1100nm.m
PRG Fig. 3.4 T_SRG_FDFD_FieldDist_1550nm.m
PRG Fig. 3.5 FDTD T_SRG_FDTD.m
PRG Fig. 3.5 FDFD T_SRG_FDFD.m
PRG Fig. 3.5 MoL T_SRG_MoL.m
PRG Fig. 3.5 RCWA T_SRG_RCWA.m
PRG Fig. 3.7 T_SRG_R_lambda_wtop.m
PRG Fig. 3.9 T_SRG_BF_R_lambda_wbottom.m
PRG Fig. 3.10 MoL T_SRG_BF_MoL_R_lambda.m PRG Fig. 3.10 RCWA T_SRG_BF_RCWA_R_lambda.m
PRG Fig. 3.11 T_SRG_DFLP_RCWA_R_lambda.m
PRG Fig. 3.12 T_SRG_DFSP_RCWA_R_lambda.m
PRG Fig. 4.1 NLFDFD_HOMCRYS.m
PRG Fig. 4.3 NLFDFD_PPLNCRYS.m
PRG Fig. 4.4 NLFDFD_ITER.m
PRG Fig. 4.8 NLFDFD_onearraycyl_FieldDIst.m PRG Fig. 4.12 NLFDFD_fivearraycyl_FieldDIst.m PRG Figs. 5.1 5.2 5.3 PSFD_CircularWaves_PointSource.m
PRG Fig. 5.12 NL_PSFD_TILTEDPPLN_NormFW.m
PRG Fig. 5.13 NL_PSFD_PPLN_TiltedFW.m
PRG Figs. 5.18 5.19 5.20 NL_PSFD_NLPC.m
Appendix B: TF-SF formulation in the FDTD method
TF-SF formulation in one dimension
Fig.6.1 shows theTF-SF domain definition for one dimension.
Fig. 6.1. TF-SF domain definition in one dimension. (This figure is a part of Fig. 5.8b in [81].)
We cannot apply theFDTDalgorithm as in the case of a normal grid, instead, we have to correct the fields at the boundary between the total field and the scattered field. At the boundary the FDTD equations contain different type of fields (scattered and total), but these equations can be corrected if we know the incoming field at the boundary.
ConsideringEz,tot|iL and theFDTD equation at grid point iL: Ez,tot|n+1i
L =Ez,tot|n+1i
L + ∆t
ε0∆
Hy,tot|n+1/2i
L+1/2−Hy,scat|n+1/2i
L−1/2
, (6.1)
obviously we have to apply the following correction:
Ez,tot|n+1i
L ={Ez,tot|n+1i
L }Y ee− ∆t
ε0∆Hy,inc|n+1/2i
L−1/2, (6.2)
and the incoming field is known.
ForHy,scat|n+1/2i
L−1/2 the FDTD equation and the correction are:
Hy,scat|n+1/2i
L−1/2 =Hy,scat|n−1/2i
L−1/2+ ∆t
µ0∆ Ez,tot|ni
L−Ez,scat|ni
L−1
, (6.3a)
Hy,scat|n+1/2i
L−1/2 ={Hy,scat|n+1/2i
L−1/2}Y ee− ∆t
µ0∆Ez,inc|ni
L. (6.3b)
Ez|iR ={Ez|iR }Y ee+
ε0∆Hy,inc|i
R+1/2 (6.4a)
Hy|n+1/2i
R+1/2 ={Hy|n+1/2i
R+1/2}Y ee+ ∆t
µ0∆Ez,inc|ni
R. (6.4b)
TF-SF formulation in two dimensions
We will repeat the procedure like in the one-dimensional case, namely using theFDTD algorithm for the whole grid to obtain field values in the next time point, but now we have to match the fields at the four TF-SF boundaries. We provide the formulas for TEz, for TMz mode they can be derived analogously. Like in one dimension we have to modify the FDTD equations at the TF-SF boundary, as an example at the bottom boundary line1:
Ez|n+1i,j
0 ={Ez|n+1i,j
0 }Y ee+ ∆t
ε0∆Hx,inc|n+1/2i,j0−1/2 (i=i0, . . . , i1) (6.5a) Hx|n+1/2i,j
0−1/2 ={Hx|n+1/2i,j
0−1/2}Y ee+ ∆t
µ0∆Ez,inc|ni,j0 (i=i0, . . . , i1) (6.5b) In order to execute the correction at TF-SF boundaries we need the values of the incident field there. To determine those values we are using an obliquely oriented one-dimensional grid with a plane wave propagating on it, and take the projections of its field values to the TF-SF boundary. This procedure is illustrated in Fig. 6.2, the point m0 of the one-dimensional grid is situated at the bottom left corner of the two-dimensional TF-SF boundary.
For the determination of the incoming field components we introduce the normalized incoming wave vector: ˆkinc = cosφˆx+ sinφˆy, and the position vector:rcomp = (icomp− i0)ˆx+ (jcomp−j0)ˆy and define the distance
d=ˆkinc·ˆrcomp, (6.6)
which corresponds to the projection of the position vector of a point at theTF-SF bound-ary to the one-dimensional grid laid through the origin (i0, j0) (see Fig. 6.2).
On the one-dimensional grid (on which∆and∆tagree with the ones used on the two-dimensional grid) wave propagation is started at the initial point of the one-two-dimensional grid using a hard source. When the incoming wave reaches the two-dimensional grid at point m0, the TF-SF boundary on the two-dimensional grid is excited via the cou-pling equations. The equations for the one-dimensional propagation compensated by the
1The derivation of formulas and the corrections on the other three boundaries can be found on page 193–196 in [81]
Fig. 6.2. Determination of the incoming field values. (The original figure can be found in [81] as Fig. 5.11a.)
anisotropy are
Einc|n+1m =Einc|nm+ 1 hv˜p(φ=0◦)
˜ vp(φ)
i
∆t ε0∆
Hinc|n+1/2m−1/2−Hinc|n+1/2m+1/2
, (6.7a)
Hinc|n+1/2m+1/2 =Hinc|n−1/2m+1/2+ 1 h˜v
p(φ=0◦)
˜ vp(φ)
i
∆t
µ0∆ Einc|nm−Einc|nm+1
. (6.7b)
The reason of the compensation is that on the one-dimensional grid axial propagation is executed, while on the two-dimensional grid it happens under an angle φ. Due to the numerical anisotropy (see in Sec. 2.1.2) the velocity of the propagation is different in the two cases (note that the same ∆and ∆t are used). Hence we have to adjust the velocity of the propagation on the one-dimensional grid to make it equal to the velocity of the propagation on the two-dimensional grid.
TF-SFcorrectionsEz,inc|nd,Hx,inc|n+1/2d ,Hy,inc|n+1/2d are calculated using the field values from the one-dimensional grid. Index d indicates the field value which is situated at a distance d from(i0, j0).
The E component of the incident field (d0 =d−[d]) is given by
Einc|nd = (1−d0)Einc|nm0+[d]+d0Einc|nm0+[d]+1, (6.8a)
Ez,inc|nd =Einc|nd. (6.8b)
Hinc|n+1/2d = (1−d0)Hinc|n+1/2m
0−1/2+[d00]+d0Hinc|n+1/2m
0+1/2+[d00], (6.9a)
Hx,inc|n+1/2d =Hinc|n+1/2d sinφ, (6.9b)
Hy,inc|n+1/2d =−Hinc|n+1/2d cosφ. (6.9c)
Appendix C: Periodic boundary condition in the FDTD method
For simulation of periodic structures we use Periodic Boundary Condition (PBC).
The sketch of the procedure is shown in Fig. 6.3. Mathematically, for a selected mode (k = 2πfc ):
Hx(x, y =yp+ ∆y/2, t) = Hx(x, y = ∆y/2, t−ypksinφ1) (6.10a) Ez(x, y = 0, t) = Ez(x, y =yp, t+ypksinφ1) (6.10b)
Fig. 6.3. PBC. The original figure can be found in [81] as Fig. 13.1.
If the angle of incidence differs from 0◦, a problem occurs due to the time shift. On the one hand a previously calculated value should be used (Hx), which can be executed, but on the other hand a future value, which is not yet computed. This is impossible (at least for general waveform, several modes). In case of normal incidence there is no such problem (φ1 = 0), and Eqs. (6.10a) and (6.10b) are simplified to the forms (which can be applied for a general waveform):
Hx(x, y =yp+ ∆y/2, t) = Hx(x, y = ∆y/2, t) (6.11a) Ez(x, y = 0, t) = Ez(x, y =yp, t) (6.11b) Although simulation of oblique incidence is not possible for general waveform (as we explained above), considering a single frequency plane wave it still can be simulated with the so-called Sine-Cosine method. Converting PBC equations to the Fourier space they take the following form:
H˘x(x, y =yp+ ∆y/2) = ˘Hx(x, y = ∆y/2)e−jkyyp (6.12a) E˘z(x, y = 0) = ˘Ez(x, y =yp)e+jkyyp, (6.12b) where ky =ksinφ1. We consider two grids and apply PBC on both ones and select the
Ez(C) =<
(Ez(A) +jEz(B))ejkyyp (6.13a) Ez(D) ==
(Ez(A) +jEz(B))ejkyyp (6.13b) Hx(A) =<
(Hx(C) +jHx(D))e−jkyyp (6.13c) Hx(B) ==
(Hx(C) +jHx(D))e−jkyyp (6.13d)
Fig. 6.4. Sine-Cosine method. The original figure can be found in [81] as Fig. 13.9.
Appendix D: Deduction of the UPML boundary condi-tion
We examine the following situation: two domains are separated atx= 0. The boundary is a plane perpendicular to the xaxis, let us select an isotropic material for the left-hand domain (x < 0), and an anisotropic one for the right-hand domain (x > 0) with the following permittivity and permeability tensors:
ε2 =ε2
a 0 0 0 b 0 0 0 b
, µ2 =µ2
c 0 0 0 d 0 0 0 d
. (6.14)
The source is an incoming plane wave from the left-hand domainH˘inc =H0e−jβ1xx−jβ1yy. In the right-hand domain Maxwell’s equations take the form:
β2×E˘ =ωµ2H,˘ (6.15a)
β2×H˘ =−ωε2E,˘ (6.15b)
where β2 =β2xˆx+β2yy, and the wave equation is obtained:ˆ
β2×(ε−12 β2)×H˘ +ω2µ2H˘ = 0. (6.16)
By substituting the material properties from Eq. (6.14) into the wave equation Eq.
(6.16) the following formulas can be deduced:
k22−β22xb−1d−1−β22ya−1d−1 = 0, (T Ez : ˘Hx,H˘y = 0), (6.17a) k22−β22
xb−1d−1−β22
yb−1c−1 = 0, (T Mz : ˘Hz = 0), (6.17b) where k22 =ω2µ22.
Taking into account the reflection, the electric and magnetic fields in the left-hand domain are
H˘1 = ˆzH0 1 + Γe2jβ1xx
e−jβ1xx−jβ1yy, (6.18a)
E˘1 =
−ˆxβ1y
ω1 1 + Γe2jβ1xx
+ ˆyβ1x
ω1 1−Γe2jβ1xx
H0e−jβ1xx−jβ1yy, (6.18b)
H2 = ˆzH0τ e , (6.19a) E˘2 =
−ˆx β2y
ω2a + ˆy β2x ω2b
H0τ e−jβ2xx−jβ2yy, (6.19b) where Γ is the reflection coefficient and τ is the transmission coefficient, they are de-termined by the continuous transition of the tangential components, which provides us:
Γ = β1x −β2xb−1
β1x +β2xb−1, τ = 1 + Γ. (6.20)
Due to phase matching at thex= 0interface (which is equivalent with the Snell-Descartes refraction law):
β2y =β1y. (6.21)
Substitution of β2y from Eq. (6.21) into Eq. (6.17a) (TEz case) results in β2x =q
k22bd−β12ya−1b. (6.22)
By selection of ε1 = ε2, and µ1 =µ2 we have k2 = k1. Furthermore, if we choose b =d and a−1 = b, it causes in Eq. (6.22): β2x = bβ1x, from which we can see Γ = 0 in Eq.
(6.20) independently from the angle of incidence and wavelength. So the latter choices correspond to Perfectly Matched Layer (PML). If we substitute β2y from Eq. (6.21) into Eq. (6.17b) (TMz case), we can conclude after similar steps that Γ = 0 is satisfied if b =d and c−1 =d. In summary the selection of b =d, a−1 =b, and c−1 =d is the PML condition.
Appendix E: Implementation of the UPML boundary condition in the FDTD method
For implementation of theUPMLboundary condition in theFDTDmethod we define the following quantities:
D˘x =εsz
sxE˘x, D˘y =εsx
syE˘y, D˘z =εsy
szE˘z. (6.23)
Using these definitions in Eq. (2.15a) we get:
∂H˘z
∂y − ∂H˘y
∂z
∂H˘x
∂z −∂H˘z
∂x
∂H˘y
∂x − ∂H˘x
∂y
=jω
sy 0 0 0 sz 0 0 0 sx
D˘x D˘y D˘z
, (6.24)
using the definition of sx, sy, sz (see Eq. (2.17)) and converting the problem to the time domain we obtain:
∂Hz
∂y − ∂Hy
∂z
∂Hx
∂z −∂Hz
∂Hy ∂x
∂x − ∂Hx
∂y
= ∂
∂t
κy 0 0 0 κz 0 0 0 κx
Dx Dy Dz
+1
ε
σy 0 0 0 σz 0 0 0 σx
Dx Dy Dz
. (6.25)
Using Eq. (6.25) we get for the time stepping ofDx (applying a semi-implicit interpo-lation like in Eq. (2.3)):
Dx|n+1i+1/2,j,k =
2κy−σy∆t 2κy+σy∆t
Dx|ni+1/2,j,k+
2∆t 2κy+σy∆t
· Hz|n+1/2i+1/2,j+1/2,k−Hz|n+1/2i+1/2,j−1/2,k
∆y − Hy|n+1/2i+1/2,j,k+1/2−Hz|n+1/2i+1/2,j,k−1/2
∆z
! .
Consider nowD˘x from Eq. (6.23) and insert the definition of {si}i=x,y,z:
κx+ σx jω
D˘x =ε
κz + σz jω
E˘x, (6.26)
multiplying by jω and converting to the time domain:
∂
∂t(κxDx) + σx
ε Dx =ε ∂
∂t(κzEx) + σz ε Ex
, (6.27)
Ex|i+1/2,j,k =
2κz+σz∆t Ex|i+1/2,j,k+
+ 1
(2κz+σz∆t)ε ·h
(2κx+σx∆t)Dx|n+1i+1/2,j,k−(2κx−σx∆t)Dx|ni+1/2,j,ki . Formulas for the rests of the components can be derived similarly.
Appendix F: Numerical phase velocity and anisotropy in the FDFD method
In TEz polarization the discretized Maxwell’s equations read
−jωµ0Hx = 1
∆y
Ez
I,J+1−Ez I,J
, (6.28a)
−jωµ0Hy =− 1
∆x
Ez
I+1,J −Ez I,J
, (6.28b)
jωε0εrEz = 1
∆x
Hy
I,J−Hy I−1,J
− 1
∆y
Hx
I,J −Hx I,J−1
. (6.28c)
By substitution of a plane wave solution Ez
I,J =Ez,0e−j(˜kx∆xI+˜ky∆yJ), (6.29)
where ˜kx and ˜ky are the components of the numerical wavenumber vector, into the Eq.
(6.28a) one finds
−jωµ0Hx,0 = 1
∆yEz,0
e−j˜ky∆y/2−ej˜ky∆y/2
, (6.30)
so Hx,0 and similarly Hy,0 are obtained Hx,0 = 2
∆yωµ0Ez,0sin (˜ky∆y/2), Hy,0 =− 2
∆xωµ0Ez,0sin (˜kx∆x/2), (6.31) and for Ez,0 we have
Ez,0 =− 2
∆xωεHy,0sin (˜kx∆x/2) + 2
∆yωεHx,0sin (˜ky∆y/2). (6.32) SubstitutingHx,0 and Hy,0 from Eq. (6.31) into Eq. (6.32) the result is
εµ= 2
∆xω 2
sin2(˜kx∆x/2) + 2
∆yω 2
sin2(˜ky∆y/2). (6.33) Using the dispersion relation ω=kc, we get an expression for k
k2 = 2
∆x 2
sin2(˜kx∆x/2) + 2
∆y 2
sin2(˜ky∆y/2) (6.34) Defining the numerical sampling density: Nλ = ∆xλ0 and assuming an equidistant grid
∆x/2 =
Nλλ0/2 =
Nλk. (6.35)
Expressing k˜ in polar coordinate system (˜kx = cosφ˜k and k˜y = sinφk), Eq. (6.34) takes˜ the following form
k2 =k2 Nλ
π 2
sin2 π Nλ
cosφ k˜ k
! +k2
Nλ π
2
sin2 π Nλ
sinφ
˜k k
! ,
and after simplification by k2 π
Nλ 2
= sin2 π
Nλ cosφk˜ k
!
+ sin2 π
Nλ sinφk˜ k
!
. (6.36)
For φ= 0◦ (major grid axis):
k˜=kNλ
π sin−1 π
Nλ. (6.37)
For φ= 45◦ (grid diagonal axis):
k˜=k
√2Nλ
π sin−1 π
√2Nλ. (6.38)
Appendix G: Code snippet for the NL-PSFD-SYLV method
A Matlab [9] code snapshot is presented for the NL-PSFD-SYLVmethod based on the modified conjugate gradient method:
C2 = -d33*mask.*Ez.*Ez;
R = C2 - A2*Ez2w - Ez2w*B2;
Q_ = A2’*R + R*B2’;
Z_ = Q_;
for i1 = 1:Niter
gamma = trace(R’*R)/trace(Z_’*Z_);
Ez2w = Ez2w + gamma*Z_;
Rp = R;
R = C2 - A2*Ez2w - Ez2w*B2;
Q_ = A2’*R + R*B2’;
eta = trace(R’*R)/trace(Rp’*Rp);
Z_ = Q_ + eta*Z_;
end
component for a centered rectangular lattice
The Fourier component ofχ(2)(r) belonging to the reciprocal vector G01 is χ(2)(G01) = 1
V Z
e−jG01rχ(2)(r)dr. (6.39)
The basis of the NonLinear Photonic Crystal (NLPC) is a circle, the nonlinear coefficient is inverted outside the circles (−χ0) and and remains unchanged (χ0) inside them:
χ(2)(r) = χ0 2Θ r2−R2
−1
, (6.40)
where the radius of the circles is denoted by R. Substituting χ(2)(r) from Eq. (6.40) into Eq. (6.39) results in:
χ(2)(G01) = 2χ0 V
Z
e−jG01rΘ r2−R2
dr− χ0 V
Z
e−jG01rdr. (6.41) The reciprocal vector is selected to be
G01= 2π a
1, −1
tanα
, G01 =|G01|= 2π
asinα. (6.42)
The second term in Eq. (6.41) has no contribution to the result, as it corresponds to the Fourier component of a homogeneous nonlinear crystal. We carry out the integration for one primitive cell:
χ0 ab
Z a 0
dxe−j2πax Z b
0
dyejatan2παy = χ0 ab
Z b 0
dyejatan2παyja 2π
e−j2πaa−1
= 0. (6.43)
The integral in the first term in Eq. (6.41) will be evaluated for the Wigner-Seitz cell.
Due to the step function in the integrand, the boundary is a circle with radius R:
2χ0 ab/2
Z
|r|<R
e−jG01rdr= 4χ0 ab
Z R 0
drr Z 2π
0
dϕe−jasinα2π rcos(ϕ−ϕ0). (6.44) Thenϕis substituted byϕ0 =ϕ−πand a trigonometric identity is applied:−cos(ϕ−ϕ0) = sin(ϕ0−ϕ0+π/2)and ϕ0 =π/2 is selected:
χ(2)(G01) = 4χ0 ab
Z R 0
drr Z π
−π
dϕ0ejasinα2π rsinϕ0. (6.45)
Application of the integral representation of the Bessel function J0(x) = 2π1 Rπ
−πdθejxsinθ yields
χ(2)(G01) = 4χ0 ab
Z R 0
drr2πJ0
2π asinαr
= 4χ0 ab 2π
asinα 2π
2Z asinα2π R 0
dxxJ0(x) (6.46) Finally, by using the Bessel function identity Rx
0 dx0x0J0(x0) =xJ1(x) one finds χ(2)(G01) = 4χ0asinα
ab RJ1
2π asinαR
. (6.47)
Acknowledgement
I would like to thank my supervisor, Zsolt Kis for his always supporting and devoted mentoring since the beginning of my MSc thesis. I have delved into the topic of optical simulation methods for his influence. Beyond giving me professional advice he also pro-moted my participation in conferences and involved me in exciting projects beyond my Ph.D. work. My thanks are also due to my internal consultant, Pál Maák, who not only performed administrative tasks but I also acquired deep knowledge in laser physics and design of laser systems from him. I wish to express my thanks to Gábor Erdei for the knowledge I have learned in lessons of general optics and optical system design and for his valuable lecture notes.
I deeply appreciate the patience and support from my family, from my wife Zsuzsanna and my children Flóra and Dániel. They made me possible to carry out all these things.
I still want to thank my parents to start on my way, the time devoted by my mother to enhance my math knowledge after school; and the time devoted by my father for the preparation for physics competitions and carrying out physical experiments together; and for my sister, who always supported me at the competitions.
I am grateful to my workplace, Semilab Co. Ltd., to Tibor Pavelka, Áron Pap and György Nádudvari that they enabled me to attend the courses of the Ph.D. school next to my work.
At last but not least I am also indebted to my physics teacher in primary school, to Győző Soleczki to have engaged my attention for physics.
Bibliography
[1] J. C. Maxwell. A treatise on electricity and magnetism. Cambridge University Press (1873).
[2] J. D. Jackson. Classical Electrodynamics, 3rd Edition. Wiley (1998).
[3] D. J. Griffiths. Introduction to Electrodynamics, 4th Edition. Pearson (2013).
[4] K. Simonyi. Elméleti Villamosságtan. Tankönyvkiadó (1988).
[5] R. C. Booton. Computational Methods for Electromagnetics and microwaves. Wiley (1992).
[6] D. B. Davidson. Computational Electromagnetics for RF and Microwave Engineer-ing, Second Edition. Cambridge University Press (2010).
[7] A. F. Peterson, S. L. Ray, and R. Mittra. Computational Methods for Electromag-netics. Wiley-IEEE Press (1997).
[8] T. Rylander, P. Ingelström, and A. Bondeson. Computational Electromagnetics.
Springer (2005).
[9] MATLAB, The MathWorks, Inc., Natick, Massachusetts, United States. http:
//www.mathworks.com.
[10] N. Bonod and J. Neauport. “Diffraction gratings: from principles to applications in high-intensity lasers.” Advances in Optics and Photonics,8(1): 156–199 (2016).
[11] M. Moharam, E. B. Grann, D. A. Pommet, et al. “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings.”
Journal of the Optical Society of America A, 12: 1068 – 1076 (1995).
[12] V. Braginsky, M. Gorodetsky, and S. Vyatchanin. “Thermodynamical fluctuations and photo-thermal shot noise in gravitational wave antennae.” Physics Letters A, 264: 1–10 (1999).
[13] G. M. Harry, A. M. Gretarsson, P. R. Saulson, et al. “Thermal noise in interfer-ometric gravitational wave detectors due to dielectric optical coatings.” Classical and Quantum Gravity,19: 897–917 (2002).
[14] E. Popov, L. Mashev, and D. Maystre. “Theoretical Study of the Anomalies of Coated Dielectric Gratings.” Optica Acta, 33: 607 – 619 (1986).
[15] I. Avrutsky and V. Sychugov. “Reflection of a Beam of Finite Size from a Corrugated Waveguide.” Journal of Modern Optics, 36: 1527 – 1539 (1989).
[16] D. Rosenblatt. “Resonant Grating Waveguide Structures.” IEEE Journal of Quan-tum Electronics, 33: 2038 – 2059 (1997).
[17] A. Bunkowski, O. Burmeister, D. Friedrich, et al. “High reflectivity grating waveg-uide coatings for1064 nm.” Classical and Quantum Gravity, 23: 7297–7303 (2006).
[18] F. Brückner, D. Friedrich, T. Clausnitzer, et al. “Demonstration of a cavity coupler based on a resonant waveguide grating.” Optics Express, 17: 163–169 (2009).
[19] F. Brückner, T. Clausnitzer, O. Burmeister, et al. “Monolithic dielectric surfaces as new low-loss lightmatter interfaces.” Optics Letters, 33: 264–266 (2008).
[20] F. Brückner, D. Friedrich, T. Clausnitzer, et al. “Realization of a monolithic high-reflectivity cavity mirror from a single silicon crystal.” Physical Review Letters, 104: 163903, 1–4 (2010).
[21] C. Lee, K. Hane, and S. Lee. “The optimization of sawtooth gratings using RCWA and its fabrication on a slanted silicon substrate by fast atom beam etching.” Journal of Micromechanics and Microengineering,18: 045014–045021 (2008).
[22] S.-D. Wu, T. K. Gaylord, and E. N. Glytsis. “Optimization of sawtooth surface-relief gratings: effects of substrate refractive index and polarization.” Applied Optics, 45: 3420 – 3424 (2006).
[23] S.-D. Wu, T. K. Gaylord, J. S. Maikisch, et al. “Optimization of anisotropically etched silicon surface-relief gratings for substrate-mode, optical interconnects.” Ap-plied Optics, 45: 15 – 21 (2006).
[24] J. Maikisch and T. Gaylord. “Optimum parallel-face slanted surface-relief gratings.”
Applied Optics, 46: 3674 – 3681 (2007).
[25] S. Kroker, T. Käsebier, F. Brückner, et al. “Reflective cavity couplers based on resonant waveguide gratings.” Optics Express, 19: 16466 – 16479 (2011).
[26] S. Kroker, F. Brückner, E.-B. Kley, et al. “Enhanced angular tolerance of resonant waveguide grating reflectors.” Optics Letters, 36: 537 – 539 (2011).
[27] S. Kroker, T. Kasebier, E.-B. Kley, et al. “Coupled grating reflectors with highly angular tolerant reflectance.” Optics Letters, 38: 3336 – 3339 (2013).
[28] H. Zhao and D. Yuan. “Design of fused-silica rectangular transmission gratings for polarizing beam splitter based on modal method.” Applied Optics, 49: 759 – 763 (2010).
[29] M. Mutlu, A. Akosman, G. Kurt, et al. “Experimental realization of a high-contrast grating based broadband quarter-wave plate.” Optics Express, 20: 27966–27973 (2012).
[30] G. Quaranta, G. Basset, O. J. F. Martin, et al. “Recent Advances in Resonant Waveguide Gratings.” Laser and Photonics Reviews,12(9): 1800017 (2018).
[31] L. Yuan and Y. Y. Lu. “Analyzing second harmonic generation from arrays of cylin-ders using Dirichlet-to-Neumann maps.” Journal of the Optical Society of America B,26(4): 587–594 (2009).
[32] A. Yariv. Quantum Electronics 3rd Edition. John Wiley & Sons (1989).
[33] R. Boyd. Nonlinear Optics 3rd Edition. Elsevier, Academic Press (2008).
[34] R. L. Sutherland, D. G. McLean, and S. Kirkpatrick. Handbook of Nonlinear Optics 2nd Edition. Marcel Dekker Inc., New York (2003).
[35] ELI-ALPS. https://www.eli-alps.hu/.
[36] J. A. Armstrong, N. Bloembergen, J. Ducuing, et al. “Interactions between Light Waves in a Nonlinear Dielectric.” Physical Review,127: 1918 (1962).
[37] M. Manzo, F. Laurell, V. Pasiskevicius, et al. “Two-dimensional domain engineering in LiNbO3 via a hybrid patterning technique.” Optical Materials Express,1(3): 365–
371 (2011).
[38] M. L. Ren and Z. Y. Li. “Giant enhancement of second harmonic generation in non-linear photonic crystals with distributed Bragg reflector mirrors.” Optics Express, 17(17): 14502–14510 (2009).
[39] M. L. Ren, D. L. Ma, and Z. Y. Li. “Experimental demonstration of super quasi-phase matching in nonlinear photonic crystal.” Optics Letters,36: 3696–3698 (2011).
[40] J. Yuan and J. Yang. “Computational design for efficient second-harmonic genera-tion in nonlinear photonic crystals.” Journal of the Optical Society of America B, 30: 205 (2013).
[41] V. Berger. “Nonlinear Photonic Crystals.” Physical Review Letters, 81(19): 4136–
4139 (1998).
[42] H. C. Liu and A. H. Kung. “Substantial gain enhancement for optical parametric amplification and oscillation in two-dimensional nonlinear photonic crystals.” Optics Express, 16(13): 9714–9725 (2008).
[43] M. Lazoul, A. Boudrioua, L. M. Simohamed, et al. “Multi-resonant optical para-metric oscillator based on 2D-PPLT nonlinear photonic crystal.” Optics Letters, 40(8): 1861–1864 (2015).
[44] M. Lazoul, A. Boudrioua, L. M. Simohamed, et al. “Experimental study of multi-wavelength parametric generation in a two-dimensional periodically poled lithium tantalate crystal.” Optics Letters,38(19): 3892–3894 (2013).
[45] I. V. Soboleva, S. A. Seregin, A. A. Fedyanin, et al. “Efficient bidirectional optical harmonics generation in three-dimensional photonic crystals.” Journal of the Optical Society of America B, 28: 1680–1684 (2011).
[46] RP Photonics Encyclopedia - Birefringence. https://www.rp-photonics.com/
birefringence.html.
[47] P. Sumithra and D. Thiripurasundari. “A review on Computational Electromagnet-ics Methods.” Advanced Electromagnetics, 6(1): 42–55 (2017).
[48] ZEMAX OpticStudio. https://www.zemax.com/.
[49] S. G. Johnson and J. D. Joannopoulos. “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis.” Optics Express, 8: 173–190 (2001).
[50] R. D. Meade, A. M. Rappe, K. D. Brommer, et al. “Accurate theoretical analysis of photonic band-gap materials.” Physical Review B, 48: 8434–8437 (1993).
[51] H. S. Sozuer and J. P. Dowling. “Photonic band calculations for woodpile struc-tures.” Journal of Modern Optics,41: 231–239 (1994).
[52] S. Shouyuan, C. Chen, and D. W. Prather. “Plane-wave expansion method for calculating band structure of photonic crystal slabs with perfectly matched layers.”
Journal of the Optical Society of America A, 21: 1769–1775 (2004).
[53] K. M. Leung and Y. F. Liu. “Photon band structures: The plane-wave method.”
Physical Review B,41: 10188–10190 (1990).
[54] S. Guo and S. Albin. “Simple plane wave implementation of photonic crystal calcu-lations.” Optics Express, 11: 167–175 (2003).
[55] J. B. Pendry. “Photonic band structures.” Journal of Modern Optics, 41: 209–229 (1994).
[56] H. Oraizi and M. Afsahi. “Analysis of planar dielectric multilayers as FSS by trans-mission line transfer matrix method (TLTMM).” Progress In Electromagnetics Re-search, 74: 217–240 (2007).
[57] Z. Li and L. Lin. “Photonic band structures solved by a plane-wave-based transfer-matrix method.” Physical Review E, 67: 046607 (2003).
[58] A. J. Ward. Transfer Matrices, Photonic Bands and Related Quantities. Ph.D.
thesis, Imperial College of Science, Technology and Medicine, London (1996).
[59] E. X. Pérez. Design, fabrication and characterization of porous silicon multilayer optical devices. Ph.D. thesis, Universitat Rovira I Virgili (2007).
[60] A. Alejo-Molina, J. J. Sánchez-Mondragón, C. Velésquez-Ordonez, et al. “Transfer Matrix and Reflexion in a Metallo-Dielectric Photonic Crystal.” AIP Conference Proceedings,992: 501 (2008).
[61] P. Bell, J. Pendry, L. M. Moreno, et al. “A program for calculating photonic band structures and transmission coefficients of complex structures.” Computer Physics Communications, 85: 306–322 (1995).
[62] R. F. Harrington. Field Computation by Moment Methods. IEEE Press (1993).
[63] W. C. Gibson. The Method of Moments in Electromagnetics. Chapmann & Hall CRC Press (2008).
[64] X.-Q. Sheng and W. Song. Essentials of Computational Electromagnetics. Wiley (2012).
[65] N. Engheta, W. D. Murphy, V. Rokhlin, et al. “The Fast Multipole Method for Electromagnetic Scattering Computation.” IEEE Transactions on Antennas and Propagation,40: 634–641 (1992).
[66] L. G. Ozgur Ergul. The Multilevel Fast Multipole Algorithm (MLFMA) for Solving Large-Scale Computational Electromagnetics Problems. Wiley-IEEE Press (2014).
[67] J. Jin. The Finite Element Method in Electromagnetics 2nd ed. Wiley (2002).
[68] J. P. Webb. “The Finite-Element Method for Finding Modes of Dielectric-Loaded Cavities.” IEEE Transactions on Microwave Theory and Techniques, 33: 635–639 (1985).
[69] M. G. Larson and F. Bengzon. The Finite Element Method: Theory, Implementa-tion, and Applications. Springer (2013).
[70] Q. H. Liu. “The PSTD algorithm: A time-domain method requiring only two cells per wavelength.” Microwave and Optical Technology Letters, 15: 158–165 (1997).
[71] Q. H. Liu. “Large-scale simulations of electromagnetic and acoustic measurements using the pseudospectral time-domain (PSTD) algorithm.” IEEE Transactions on Geoscience and Remote Sensing,37: 917–926 (1999).
[72] Q. H. Liu and G.-X. Fan. “A frequency-dependent PSTD algorithm for general dispersive media.” IEEE Microwave and Guided Wave Letters, 9: 51–53 (1999).
[73] W. Sun, G. Videen, Q. Fu, et al. “Scattered-field FDTD and PSTD algorithms with CPML absorbing boundary conditions for light scattering by aerosols.” Journal of Quantitative Spectroscopy and Radiative Transfer, 131: 166–174 (2013).
[74] G. L. Pedrola. Beam Propagation Method for Design of Optical Waveguide Devices.
Wiley (2015).
[75] J. M. Liu and L. Gomelsky. “Vectorial beam propagation method.” Journal of the Optical Society of America A,9: 1574–1585 (1992).
[76] X. Ma, S. Zhu, L. Li, et al. “Finite-difference beam-propagation method for anisotropic waveguides with torsional birefringence.” Optics Express,26: 3995–4003 (2018).
[77] H. El-Refaei, D. Yevick, and I. Betty. “Stable and Noniterative Bidirectional Beam Propagation Method.” IEEE Photonics Technology Letters, 12(4): 389–391 (2000).
[78] P. L. Ho and Y. Y. Lu. “A Stable Bidirectional Propagation Method Based on Scat-tering Operators.” IEEE Photonics Technology Letters,13(12): 1316–1318 (2001).
[79] Computational Electromagnetics - EMPossible. https://empossible.net/
academics/emp5337/.
[80] K. Yee. “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media.” IEEE Transactions on Antennas and Propagation, 14 (3): 26274 – 26284 (1966).
[81] A. Taflove. Computational Electrodynamics The Finite-Difference Time-Domain Method 3rd ed. Artech House (2005).
[82] J.-P. Berenger. “A Perfectly Matched Layer for the Absorption of Electromagnetic Waves.” Journal of Computational Physics,114: 185–200 (1994).
[83] R. Rumpf. “Simple implementation of arbitrarily shaped total-field/scattered-field regions in finite-difference frequency-domain.” Progress In Electromagnetics Re-search B,36: 221 – 248 (2012).
[84] R. Rumpf. Design and optimization of nano-optical elements by coupling fabrication to optical behavior. Ph.D. dissertation, University of Central Florida (2006).
[85] A. Farjadpour, D. Roundy, A. Rodriguez, et al. “Improving accuracy by subpixel smoothing in the finite-difference time domain.” Optics Letters, 31(20): 2972–2974 (2006).
[86] M. N. O. Sadiku. Numerical Techniques in Electromagnetics 2nd Edition. Boca Raton, FL: CRC Press (2001).
[87] R. Pregla and W. Pascher. The method of lines, Numerical Techniques for Mi-crowave and Millimeter Wave Passive Structures. T. Itoh (ed.), 381-446, J. Wiley Publ. (1989).
[88] R. Pregla. “The method of lines for the analysis of dielectric waveguide bends.”
Journal of Lightwave Technology, 14(4): 634–639 (1996).
[89] S. F. Helfert. “Applying Oblique Coordinates in the Method of Lines.” Progress In Electromagnetics Research, 61: 271–278 (2006).
[90] R. Rumpf. “Improved formulation of scattering matrices for semi-analytical meth-ods that is consistent with convention.” Progress In Electromagnetics Research B, 35: 241 – 261 (2011).
[91] E. F. Beckenbach and R. Redheffer. Modern mathematics for the engineer. Part 3. Difference equations and functional equations in transmission-line theory. New York, McGraw-Hill (1961).
[92] L. Li. “Formulation and comparison of two recursive matrix algorithms for mod-eling layered diffraction gratings.” Journal of the Optical Society of America A, 13(5): 1024–1035 (1996).
[93] M. G. Moharam and T. K. Gaylord. “Rigorous coupled-wave analysis of planar-grating diffraction.” Journal of the Optical Society of America, 71(7): 811–818 (1981).
[94] M. G. Moharam and T. K. Gaylord. “Rigorous coupled-wave analysis of grating diffraction - E-mode polarization and losses.” Journal of the Optical Society of America, 73(4): 451–455 (1983).
[95] M. G. Moharam and T. K. Gaylord. “Three-dimensional vector coupled-wave analysis of planar-grating diffraction.” Journal of the Optical Society of America, 73(9): 1105–1112 (1983).
[96] M. G. Moharam and T. K. Gaylord. “Rigorous coupled-wave analysis of metallic surface-relief gratings.” Journal of the Optical Society of America A, 3(11): 1780–
1787 (1986).
[97] P. Lalanne and G. M. Morris. “Highly improved convergence of the coupled-wave method for TM polarization.” Journal of the Optical Society of America A, 13(4): 779–784 (1996).
[98] M. Moharam, D. Pommet, and E. Grann. “Stable implementation of the rigorous coupled-wave analysis for surface-relief gratings: enhanced transmittance matrix approach.” Journal of the Optical Society of America A,12: 1077 – 1086 (1995).
[99] M. A. Green and M. J. Keevers. “Optical properties of intrinsic silicon at 300 K.”
Progress in Photovoltaics, 3: 182 – 192 (1995).
[100] H. H. Li. “Refractive index of silicon and germanium and its wavelength and tem-perature derivatives.” Journal of Physical and Chemical Reference Data, 9: 561 – 658 (1993).
[101] R. H. Byrd, J. C. Gilbert, and J. Nocedal. “A Trust Region Method Based on Interior Point Techniques for Nonlinear Programming.” Mathematical Programming,89: 149 – 185 (2000).
[102] R. A. Waltz, J. L. Morales, J. Nocedal, et al. “An interior algorithm for nonlin-ear optimization that combines line snonlin-earch and trust region steps.” Mathematical Programming, 107: 391 – 408 (2006).
[103] J. D. McMullen. “Optical parametric interactions in isotropic materials using a phase-corrected stack of nonlinear dielectric plates.” Journal of Applied Physics, 46: 3076–3081 (1975).
[104] M. M. Fejer, G. A. Magel, D. H. Jundt, et al. “Quasi-Phase-Matched Second Har-monic Generation: Tuning and Tolerances.” IEEE Journal of Quantum Electronics, 28(11): 2631–2654 (1992).
[105] V. G. Dmitriev, G. G. Gurzadyan, and D. N. Nikogosyan. Handbook of Nonlinear Optical Crystals 3rd Edition. Springer (1999).
[106] S. Enoch and H. Akhouayri. “Second-harmonic generation in multilayered devices:
theoretical tools.” Journal of the Optical Society of America B, 15(3): 1030–1041 (1998).
[107] M. Bertolotti. “Wave interactions in photonic band structures: an overview.” Journal of Optics A, 8(4): S9 (2006).
[108] M. L. Ren and Z. Y. Li. “Exact iterative solution of second harmonic generation in quasi-phase-matched structures.” Optics Express, 18(9): 7288–7299 (2010).
[109] J. V. Roey, J. van der Donk, and P. E. Lagasse. “Beam-propagation method: analysis and assessment.” Journal of the Optical Society of America,71(7): 803–810 (1981).
[110] H.-P. Nolting and R. Marz. “Results of Benchmark Tests for Different Numerical BPM Algorithms.” Journal of Lightwave Technology,13(2): 216–224 (1995).
[111] P. S. Weitzman and U. Österberg. “A Modified Beam Propagation Method to Model Second Harmonic Generation in Optical Fibers.” IEEE Journal of Quantum Electronics,29(5): 1437–1443 (1993).
[112] A. Locatelli, F. M. Pigozzo, F. Baronio, et al. “Bidirectional beam propagation method for secondharmonic generation in engineered multilayered waveguides.” Op-tical and Quantum Electronics, 35: 429–452 (2003).
[113] A. Locatelli, F. M. Pigozzo, D. Modotto, et al. “Nonlinear bidirectional beam prop-agation method based on scattering operators for periodic microstructured waveg-uides.” Journal of the Optical Society of America B, 20(8): 1724–1731 (2003).
[114] H. M. Al-Mudhaffar, M. A. Alsunaidi, and H. M. Masoudi. “Full-wave Solution of the Second Harmonic Generation Problem Using a Nonlinear FDTD Algorithm.”
In “Progress In Electromagnetic Research Symposium,” (2007).
[115] A. Bourgeade and E. Freysz. “Computational modeling of second-harmonic gen-eration by solution of full-wave vector Maxwell equations.” Journal of the Optical Society of America B, 17(2): 226–234 (2000).
[116] T. W. Lee and S. C. Hagness. “Pseudospectral time-domain methods for model-ing optical wave propagation in second-order nonlinear materials.” Journal of the Optical Society of America B, 21(2): 330–342 (2004).
[117] T. W. Lee and S. C. Hagness. “Pseudospectral time-domain analysis of second-harmonic generation in 2D nonlinear photonic crystals.” In “The 16th Annual Meeting of the IEEE Lasers and Electro-Optics Society (LEOS),” (2003).
[118] F. L. Teixeira. “Time-Domain Finite-Difference and Finite-Element Methods for Maxwell Equations in Complex Media.” IEEE Transactions on Antennas and Prop-agation, 56(8): 2150–2166 (2008).
[119] S. V. Makarov, M. I. Petrov, U. Zywietz, et al. “Efficient Second-Harmonic Gen-eration in Nanocrystalline Silicon Nanoparticles.” Nano Letters, 17(5): 3047–3053 (2017).
[120] M. Luo and Q. H. Liu. “Enhancement of second-harmonic generation in an air-bridge photonic crystal slab: simulation by spectral element method.” Journal of the Optical Society of America B, 28(12): 2879–2887 (2011).
[121] C. Forestiere, A. Capretti, and G. Miano. “Surface integral method for second har-monic generation in metal nanoparticles including both local-surface and nonlocal-bulk sources.” Journal of the Optical Society of America B, 30: 2355–2364 (2013).
[122] D. C. Marinica, A. G. Borisov, and S. V. Shabanov. “Second harmonic generation from arrays of subwavelength cylinders.” Physical Review B, 76: 085311 (2007).
[123] Y. Ohmura and Y. Okamura. “Staggered Grid Pseudo-spectral Time-domain Method for Light Scattering Analysis.” PIERS Online, 6(7): 632–635 (2010).
[124] P. G. Dinesen, J. S. Hesthaven, J. P. Lynov, et al. “Pseudospectral method for the analysis of diffractive optical elements.” Journal of the Optical Society of America A, 16(5): 1124–1130 (1999).