TamásSzarvas Simulationofelectromagneticwavepropagationinlinearandnonlinearstructureddielectrics Ph.D.Dissertation

(1)

Simulation of electromagnetic wave

propagation in linear and nonlinear structured dielectrics

Tamás Szarvas

Supervisor: Zsolt Kis, Ph.D.

Wigner Research Center for Physics Department of Quantum Optics and Quantum Informatics

Budapest University of Technology and Economics 2020

(2)

(3)

List of Figures

1.1 Grating structure showing different orders of diffracted modes, the polarization and thexcomponent of the wave vector. The original figure can be

found in [10]. . . 18

1.2 Categorization of gratings. . . 20

1.3 A typical binary grating illuminated with arbitrary angle of incidence and polarization. The original figure can be found in [11] as Fig. 1. . . 21

1.4 Illustration of a RWG. The original figure can be found in [30] as Fig. 2. . 22

1.5 RWG resonance on the θ_inc−λ₀/Λ map calculated with both methods. . . 24

1.6 Dependence of ordinary and extraordinary refractive indices on the wavelength. The original figure can be found in [46]. . . 29

1.7 Geometry of type I SHG in a negative uniaxial crystal. Polarization states of the fields are shown, and cdenotes the orientation of the crystal. . . 29

1.8 PPLN crystal for collinear quasi phase matching. Arrows indicate domain inversion. . . 31

1.9 Conversion efficiency for PPLN crystal and homogeneous crystal. . . 32

1.10 Intensity amplification for different order QPM structures and the ideal phase matching. The original figure can be found in [34] as Fig. 25. . . 33

1.11 Phase matching in tilted QPM grating. . . 33

2.1 Categorization of CEM. . . 36

2.2 Advantages and disadvantages of the four selected methods. . . 37

2.3 Route from FDTD to RCWA via FDFD and MoL methods. . . 38

2.4 The unit cell of Yee grid. (Axes are normalized by the grid constant.) . . . 39

2.5 Arrangement of field components in two dimensions for TMz mode. . . 41

2.6 Dependence of phase velocityv_p on propagation angleφand sampling density N_λ (S=0.5). . . 43

2.7 Dependence of phase velocityv_pon propagation angleφand stability factor S (N_λ=5). . . 44

(8)

2.8 Illustration of the numerical anisotropy. . . 44

2.9 Comparison of stable and unstable propagation. . . 45

2.10 TF-SF grid in two dimensions. . . 47

2.11 One-dimensional wave propagation using a TF-SF source.. . . 47

2.12 Scattering of a TE_z mode on a square-shaped metallic plate. . . 48

2.13 Scattering of a TM_z mode on a round-shaped dielectric (n = 3.5) cylinder. 49 2.14 Graphical interpretation of UPML boundary.. . . 50

2.15 Rearrangement of fields. . . 52

2.16 Realization of TF-SF source condition by masking procedure.f_s(x, y) is a plane wave,Q(x, y)is a mask matrix containing zeros in the TF region and ones in the SF region.<[b(x, y)] is the real part of the source term defined by Eq. (2.36) and hx(x, y) is the resulting field distribution.. . . 57

2.17 Scattering matrix approach. . . 61

3.1 Initial geometry of the SRG under consideration. . . 68

3.2 Reflectivity map of the SRG: the horizontal and vertical axes are the wavelength of the illuminating field and the angle of incidence, respectively. . . . 69

3.3 Simulation layout for the FDTD method. . . 70

3.4 The distribution of the electric and magnetic fields in the structure at 1100 nm (R ∼= 0%, left-hand panel) and at 1550 nm (R ∼= 100%, right-hand panel). The magnetic field strength is displayed as a density plot and the electric field as a vector field, the magnitude of the fields are normalized to the incident magnetic and electric fields, respectively. . . 71

3.5 Direct comparison of reflectivity plots resulting from independent numerical methods for TM polarized field (colored curves). Throughout the compu- tations the initial SRG geometry was used and the grating was illuminated at normal incidence. In this structure the reflection of TE polarized field (dashed curve) is much smaller. . . 72

3.6 Optimization steps for the widths of top and bottom layers. . . 74

3.7 Reflectivity map in terms of wavelength and width of top layer: R(λ₀,w_top). 76 3.8 Reflectivity comparison for the original and optimized structures. . . 76

3.9 Transmission map with the bottom layer width as scanning parameter. . . 78

3.10 Transmission of the bandpass filter. . . 78

3.11 Transmission of the longpass dichroic beamsplitter. . . 80

3.12 Transmission of the shortpass dichroic beamsplitter. . . 81

(9)

4.1 Conversion efficiency for homogeneous crystal. . . 90

4.2 Simulation layout for ω mesh. . . 91

4.3 Conversion efficiency for PPLN crystal - simulated and theoretical result. . 92

4.4 Convergence of the iteration process. . . 92

4.5 Presence of backward SHG in conversion efficiency: (a) l_c^F/l_c^B = 3, (b) l_c^F/l_c^B= 5. . . 94

4.6 Simulation layout for five cylinders to be infinitely repeated in thexdirection. 95 4.7 Reflectivity spectrum of one array of cylinders. . . 96

4.8 Magnitude ofE_zfield components for FW (a) and SHW (b) at the resonant frequency for one array of cylinders. . . 96

4.9 Normalized power reflectivity of SHW for one array of cylinders. . . 97

4.10 Reflectivity spectrum of five arrays of cylinders. . . 97

4.11 Normalized power reflectivity of SHW for five arrays of cylinders. . . 97

4.12 Magnitude ofE_zfield components for FW (a) and SHW (b) at the resonant frequency for five arrays of cylinders. . . 98

5.1 Simulation layout of circular wave propagation from a hard point source. The red and blue lines indicate the axial and diagonal directions, respectively.103 5.2 E_z field of a hard point source from FDFD simulation: theoretical result and simulation along grid axis and diagonal, with fitting curves. . . 104

5.3 E_z field of a hard point source from PSFD simulation: theoretical result and simulation along grid axis and diagonal, with fitting curves. . . 104

5.4 Normalized phase velocities: FDFD axial direction simulated results plus analytical curve, FDFD diagonal direction simulated results plus analytical curve, PSFD simulated results along axial and diagonal direction (which overlap each other). . . 105

5.5 Anisotropy: FDFD simulated results plus analytical curve, PSFD simulated result. . . 105

5.6 Simulation of SHG in a homogeneous crystal with the NL-FDFD method (N_λ = 100): the simulation result, the fitted curve coinciding with the simulation result, and the analytical curve are plotted. . . 107

5.7 Same as Fig. 5.6 but using the NL-PSFD method with N_λ = 30. . . 108

5.8 Fitted amplitude (peak value) of the conversion efficiency normalized to the analytical amplitude as a function of grid sampling density for NL-FDFD and NL-PSFD methods. . . 109

(10)

5.9 Fitted numerical coherence length of the conversion efficiency normalized to the analytical coherence length as a function of grid sampling density for NL-FDFD and NL-PSFD methods. . . 109 5.10 Rearrangement of the block matrixPS^e/h_x/y for the case of three periods. . . 111 5.11 Simulation layout for tilted PPLN grating indicating phase matching con-

dition. . . 111 5.12 Normal incidence with tilted PPLN grating. Image on the left shows the

normally incident FW, in the middle the tilted PPLN grating, and on the right the generated SHW. . . 112 5.13 Oblique incidence with parallel PPLN grating. Image on the left shows the

obliquely incident FW, in the middle the paralel PPLN grating, and on the right the generated SHW. . . 113 5.14 Conversion efficiencies for both types of simulations together with the the-

oretical curve. . . 114 5.15 Propagation angles of generated SHWs for both types of simulations to-

gether with the calculated angle by phase matching criterion. . . 114 5.16 Phase matching in 2D NLPC structure. . . 117 5.17 Normalized SHG nonlinear coefficient as the function of the filling factor. . 119 5.18 Simulation result at the end face of the crystal (the distance is measured in

units of mesh ∆x). Image on the left shows the obliquely incident FW, in the middle the two-dimensionally structured NLPC, and on the right the generated SHW. . . 120 5.19 Propagation angle of generated SHW in 2D NLPC. . . 121 5.20 Conversion efficiencies for 2D NLPC. . . 121 6.1 TF-SF domain definition in one dimension. (This figure is a part of Fig.

5.8b in [81].) . . . 133 6.2 Determination of the incoming field values. (The original figure can be

found in [81] as Fig. 5.11a.) . . . 135 6.3 PBC. The original figure can be found in [81] as Fig. 13.1. . . 137 6.4 Sine-Cosine method. The original figure can be found in [81] as Fig. 13.9. . 138

(11)

List of Tables

3.1 Summary chart of simulational parameters for the methods. . . 72 3.2 Initial and optimized parameters for broadband mirror including the con-

strains for the optimization . . . 75 3.3 Initial and optimized parameters for bandpass filter including the con-

strains for the optimization . . . 77 3.4 Initial and optimized parameters for the long- and shortpass dichroic filter

including the constrains for the optimization . . . 79 5.1 Normalized phase velocities for axial and diagonal directions and numerical

anisotropy for the FDFD method. In case of the PSFD method, all values are 100.0%. . . 106 5.2 Fitted amplitude (peak value) and numerical coherence length of the con-

version efficiency normalized to analytical values for NL-FDFD and NL- PSFD methods. . . 110

(12)

(13)

Abbreviations

ABC Absorbing Boundary Condition

BPM Beam Propagation Method

CEM Computational ElectroMagnetics

CPML stretched-Coordinate Perfectly Matched Layer

EM ElectroMagnetic

FDFD Finite Difference Frequency Domain

FDTD Finite Difference Time Domain

FEM Finite Element Method

FW Fundamental Wave

LN Lithium Niobate

MoL Method of Lines

NIR Near-InfraRed

NL-FDFD NonLinear Finite Difference Frequency Domain NL-FDTD NonLinear Finite Difference Time Domain NL-PSFD NonLinear PseudoSpectral Frequency Domain

NL-PSFD-SYLV NonLinear PseudoSpectral Frequency Domain Sylvester equation based

NLPC NonLinear Photonic Crystal

PBC Periodic Boundary Condition

PEC Perfect Electric Conductor

PMC Perfect Magnetic Conductor

PML Perfectly Matched Layer

(14)

PSTD PseudoSpectral Time Domain

QPM Quasi Phase Matching

RCWA Rigorous Coupled Wave Analysis

RWG Resonant Waveguide Grating

SF Scattered Field

SHG Second Harmonic Generation

SHW Second Harmonic Wave

SRG Surface Relief Grating

TF Total Field

TF-SF Total Field - Scattered Field

TE Transverse Electric

TM Transverse Magnetic

TMM Transfer Matrix Method

UPML Uniaxially Perfectly Matched Layer

(15)

Notations

σ Scalar quantity

E 3D vector

E_x Component of a 3D vector, funtion of r

E_x Component of a discretized, multiple-index quantity

ε 3D 3x3 tensor

E˘ Complex 3D vector (phasor)

ˆ

x 3D unit vector

D_e Normal (not block) matrix: N_x×N_y,N_x×N_x,N_y×N_y, . . . e_z Column vector with size N_x, N_y, . . .

A Block matrix: N_xN_y×N_xN_y

e_z Block column vector with sizeN_xN_y

I_N Identitity matrix: N ×N

· Pointwise matrix multiplication (Hadamard product)

⊗ Kronecker product

~ Redheffer Star product

d_ijk Three-index tensor

E² E² tensor, a column vector with six specific components e Block column vector with 3 components, where each com-

ponent is a a N_xN_yN_z column vector

C Block matrix 3×3, where each component is a NxNyNz × N_xN_yN_z block matrix

d Three-index tensor, where each component is a N_xN_yN_z × N_xN_yN_z block matrix

(16)

(17)

Chapter 1 Introduction

One of the oldest and most complete parts of physics is the description of electromagnetism by Maxwell’s equations [1]. This is a system of partial differential equations, which can describe static and dynamic problems as well. The solution of a specific problem strongly depends on the boundary conditions. There are several analytical solutions for specific problems [2–4], but when the geometry is complex as in structured dielectrics, numerical simulations become unavoidable [5–8]. For numerical calculation discretization has to be carried out, the question is at what stage of the calculation. Semi-analytical methods do the job analytically until a point where usually integration of a formula should be executed numerically. In another approach discretization is carried out from the very beginning at the level of Maxwell’s equations.

In this Thesis I am focussing on specific fields of computational electromagnetism.

On the one hand, I am going to study the extension of a monolithic double-layer (T- shaped) Surface Relief Grating (SRG) to specific applications such as a broadband mirror, a bandpass filter and short- and longpass dichroic beamsplitters. As an introduction, in Sec. 1.1.1 an overview concerning the principle of gratings will be given, including their classification. Next in Sec. 1.1.2 the so-called monolithic type is detailed including the binary and higher order SRG. Finally in Sec. 1.1.3 a special type of grating, the so-called Resonant Waveguide Grating (RWG) will be discussed thoroughly. The following four methods are used for the simulation and the optimization of the geometry of the SRG:

the Finite Difference Time Domain (FDTD) method, the Finite Difference Frequency Domain (FDFD) method, the Method of Lines (MoL), and the Rigorously Coupled Wave Analysis (RCWA) method; they are classified among other methods of Computational ElectroMagnetics (CEM) in Fig.2.1, and detailed in Chapter2. The initialSRGstructure, the steps of the optimization, and the optimized structures are presented as my first contribution in Chapter 3.

(18)

On the other hand, a more challenging and up-to-date topic inCEM is the numerical simulation of nonlinear processes in structured nonlinear dielectrics. So a further moti- vation or goal of this work is to develop novel methods for the simulation of nonlinear processes in second harmonic crystals. After a presentation of the basic principles of the field in Sec. 1.2, we work out a novel method as an extension of the linearFDFD method for the simulation of Second Harmonic Generation (SHG) processes. The formulation of this method and detailed simulational examples are included in Chapter 4.

In former works the generalization of the FDTD method, leading to the PseudoSpec- tral Time Domain (PSTD) method has already been proposed by replacing numerical derivation via finite differences by a pseudospectral treatment. However, this generalization of the FDFD method has not been executed yet in the literature. This is presented in Chapter 5 and the resulting method is called the PseudoSpectral Frequency Domain (PSFD) method. The method is extended for nonlinear simulation as well and several numerical examples are presented to show its effectiveness.

In this Thesis for all numerical calculations and simulations I used computer programs implemented by myself in Matlab [9]. In Appendix A the source codes of the programs are listed generating the figures for major achievements. The reference of the source code for those figures is always given in the description of the figure. The files of the programs are enclosed for this Thesis (either sent in email or attached in a memory stick).

1.1 Linear structured dielectrics

The basic equations of electrodynamics were derived by Maxwell; this is a classical field theory suitable for the description of any classical electromagnetic phenomena. Their differential and integral forms are given by

∂B

∂t =−∇ ×E−M, ∂

∂t Z

F

BdF=− I

∂F

Edl− Z

F

MdF, (1.1a)

∂D

∂t = +∇ ×H−J, ∂

∂t Z

F

DdF= + I

∂F

Hdl− Z

F

JdF, (1.1b)

∇D = 0,

I

F

DdF= 0, (1.1c)

∇B = 0,

I

F

BdF= 0, (1.1d)

(19)

where B is the magnetic inductivity, H is the magnetic field,M is the magnetism of the medium,Eis the electric field,Dis the electric displacement, andJis the electric current.

The symbol F denotes a surface; ∂F denotes the perimeter of the surface. In Eq. (1.1c) we neglect the presence of free space charge. Eqs. (1.1a)–(1.1d) have to be completed with the material equations. In case of linear media:

D =εE, B =µH; (1.2a)

J =J^s+σE, M =M^s+σ^∗H; (1.2b)

where ε is the electric permittivity, µ is the magnetic permeability, J^s and M^s are the external electric and magnetic current densities, respectively. Finally, σ and σ^∗ are the electric conductivity and the equivalent magnetic loss.

In two dimensions, assuming translation invariance in the z-direction, the fields vary only in thex−yplane. Therefore, the equations are decoupled into two groups: Transverse Electric (TE) and Transverse Magnetic (TM) modes.

Equations for theTE_z mode:

∂H_x

∂t = 1 µ

−∂E_z

∂y −(M_x^s+σ^∗Hx)

, (1.3a)

∂H_y

∂t = 1 µ

∂E_z

∂x −(M_y^s+σ^∗Hy)

, (1.3b)

∂E_z

∂t = 1 ε

∂H_y

∂x −∂H_x

∂y −(J_z^s+σEz)

. (1.3c)

Equations for theTM_z mode:

∂Ex

∂t = 1 ε

∂Hz

∂y −(J_x^s+σE_x)

, (1.4a)

∂Ey

∂t = 1 ε

−∂Hz

∂x −(J_y^s+σE_y)

, (1.4b)

∂Hz

∂t = 1 µ

∂Ex

∂y −∂Ey

∂x −(M_z^s+σ^∗H_z)

. (1.4c)

In one dimension, Eqs. (1.3a)–(1.3c) and Eqs. (1.4a)–(1.4c) simplify even further.

Assuming field propagation in the x direction one obtains:

∂H_y

∂t = 1 µ

∂E_z

∂x −(M_y^s+σ^∗H_y)

∂E_z

∂t = 1 ε

∂H_y

∂x −(J_z^s+σE_z)

, (1.5a)

∂E_y

∂t = 1 ε

−∂H_z

∂x −(J_y^s+σE_y)

∂H_z

∂t = 1 µ

−∂E_y

∂x −(M_z^s+σ^∗H_z)

. (1.5b)

(20)

1.1.1 Diffraction gratings

Diffraction gratings are special optical components having a periodic modulation at the wavelength scale on an interface between two or more materials. A very good summary on gratings can be found in Ref. [10]. Light impinging on the periodically modulated surface is reflected or transmitted in specific angles only, which is not the case if the modulation is aperiodic. By exploiting the periodicity and determining the condition of the constructive interference one can deduce the grating equations:

n_refsinθ_m =n_incsinθ_inc−mλ₀

Λ, (1.6a)

n_trnsinθ_m =n_incsinθ_inc−mλ₀

Λ, (1.6b)

whereλ₀is the wavelength in vacuum,Λis the grating period,θ_incis the angle of incidence, m is the grating diffraction order,θ_m is the propagation angle of the diffracted order, and n_inc,n_ref, n_trn are the incident, reflective, and transmissive refractive indices respectively.

Fig.1.1 shows an example for a grating indicating different orders of diffraction. BothTE and TM polarizations are illustrated.

Fig. 1.1. Grating structure showing different orders of diffracted modes, the polarization and the xcomponent of the wave vector. The original figure can be found in [10].

Mathematically a grating generates an infinite number of orders of diffraction. But in practice only some of them are propagating while the others are evanescent. The incident field at the grating surface is E_inc(x, y) = E₀e^j(kînc^x ^x+kînc^y ^y), where k_xînc = k₀n_incsinθ_inc, k_yînc =k₀n_inccosθ_inc, and k₀ = 2π/λ₀ is the wavenumber in vacuum. The phase matching

(21)

condition for the grating provides k_x,m^grat=k_x^inc−m2π

Λ , (1.7)

so k_x,m is always real. Together with the dispersion relation the y component of the wave vector can be calculated too:

k₀n_grat

2

=

k_x,m^grat 2

+

k_y,m^grat 2

, (1.8a)

k^grat_y,m = r

k₀n_grat2

− k^gratx,m

2

. (1.8b)

On the right-hand side of Fig.1.1 the propagating modes are presented; for the given choice of parameters only the 1^st, -1^st, and -2^nd orders are propagating. If we consider the x axis projection of the second order mode (see Eq. (1.7)) it would be longer than the wave vector itself (so not possible geometrically); in this case the sign of the terms under the root becomes negative in Eq. (1.6b) and so k_y becomes complex, i.e. this mode is evanescent. The longer the period length is, the more diffracted modes are propagating, so gratings can be classified as follows:

1. Subwavelength gratings: Λ < λ₀

n_inc, only specular order is propagating (as if the grating would not be present).

2. Low order gratings: Λ> λ₀

n_inc, only few reflected diffracted modes are propagating.

3. High order gratings : Λ λ₀

n_inc, lots of reflected diffractive modes are propagating.

There are several types of diffraction gratings. They can mainly be categorized as reflective or transmissive, or as volume or surface relief gratings. Fig.1.2 shows also some further categorization of diffraction gratings. In this work we are dealing with reflective type, surface relief, monolithic gratings as marked in the figure.

1.1.2 Surface relief gratings

One important category of diffraction gratings are Surface Relief Gratings (SRG), where the diffractive layer is a periodically modulated groove structure on top of the homogeneous substrate. SRGs can be made up as multilayer sandwich structures where high and low contrast layers alternate, the top layer includes the grooves. If the refractive index contrast is realized only by microfabrication from the same material (substrate) the grating is called monolithic. The simplest type of monolithic gratings is the binary grating,

(22)

Fig. 1.2. Categorization of gratings.

where the structure contains three regions: incident region (usually air, n=1), grating with certain filling factor, depending on the width of the grooves; and the substrate (n_s). The height of the grooves is denoted by d, if d is of the same size as the period or higher (d ≥ Λ), the grating is deep, if d Λ the grating is shallow. Fig. 1.3 depicts a typical binary grating for arbitrary angle of incidence and polarization.

So far we have concentrated on the diffraction properties of aSRG, but other aspects are also of interest. If we calculate average diffractive indices for different layers in the monolithic structure the grating can be considered as a multilayer high-low index device like an optical coating. Nowadays the application of antireflection coatings for various optical elements is a standard procedure. Devices for many purposes require maximum transmission and minimum reflection. For some products the coating is the active component, which provides the desired feature, for instance a broadband dielectric mirror, a dichroic beamsplitter or a bandpass filter. There is no lateral change in the thickness and refractive index inside a given layer of the coating, which therefore can be simulated as a one-dimensional structure. Despite the fact that very good optical properties can be achieved using these multilayer dielectric structures, they have some drawbacks.

Thermally driven motion of the optical elements gives rise to noise in high precision in- terferometers [12], and it turns out that mirrors with multilayer dielectric coating exhibit higher thermal noise than monolithic ones [13]. This thermal noise of the dielectrically coated mirrors mainly comes from the mechanical loss of the multilayer coating structure.

(23)

Fig. 1.3. A typical binary grating illuminated with arbitrary angle of incidence and polarization. The original figure can be found in [11] as Fig. 1.

High-reflectivity mirrors can be obtained from monolithic mediums by preparing aSRGon their surface, which has been demonstrated both theoretically and experimentally. Some early works on anomalous features of dielectric gratings have been studied in Refs. [14,15]

and a review has been published in Ref. [16].

A binary grating mirror has been proposed in Ref. [17], and subsequently demonstrated experimentally in Ref. [18], where the authors reported 99% reflectivity. A two layer SRG structure has been proposed in Ref. [19] with 100% theoretical reflectivity, in the experimental realization [20] a 99.79% reflection has been measured. In this latter case the grating ridges have a T-shape, which results in a high-low-high effective refractive index sequence. The light interference between the layers yields a highly reflective surface in a wide spectral range.SRG structures have been deeply examined for expedience in optical systems. Numerous proposals and successful experimental realizations of specific optical elements have been described based on SRGs. The scattering properties of sawtooth [21]

and right-angle-face trapezoidal SRGs have been compared in Ref. [22]. Substrate-mode optical interconnects are obtained from a V-groove siliconSRG[23], or from a parallel-face slanted SRG [24]. Cavity couplers with high angular tolerance have been demonstrated in Refs. [25–27]. Polarizing beamsplitter [28] and quarter waveplate [29] have also been developed based on SRGs.

(24)

1.1.3 Resonant waveguide gratings

The diffracted modes are present not only in the reflection region but also in the transmission region (so all regions contain them). An extremely interesting phenomenon occurs when a waveguide compatible structure is placed close to the grating. Under certain circumstances the diffracted modes can be coupled into the waveguide which leads to the presence of guided modes. When guided modes are generated reflection (or transmission) peaks appear in the spectra, which can be really narrow and high (maximal reflection or transmission). A very good summary and review on these phenomena can be found in Ref. [30]. Fig.1.4 illustrates a standard Resonant Waveguide Grating (RWG).

Fig. 1.4. Illustration of a RWG. The original figure can be found in [30] as Fig. 2.

Considering shallow grating depth d, the waveguide mode is weakly perturbed by the grating (so scattering by grating ridges and grooves can be ignored), and can be approx- imated by a propagating mode of a pure slab waveguide. Under these circumstances, the equation describing the modes in a slab waveguide can be coupled with the diffraction grating equation, by setting the propagation wavevector of the mode in the slab waveguide equal to the wave vector of the light diffracted by the grating. This yields the following equation in TE mode (see Ref. [30]):

tanκ_mt= κ_m(γ_m+δ_m)

κ²_m−γ_mδ_m , (1.9)

where κ_m =

q

n²_wk²₀−β_m², (1.10a)

γ_m = q

β_m² −n²_supk₀², (1.10b)

δ_m = q

β_m² −n²_subk₀², (1.10c)

β_m =k₀

n_supsinθ−mλ₀ Λ

, (1.10d)

(25)

whereθ is the polar angle of illumination, tis the waveguide thickness, and n_sup,n_w,n_sub are the refractive indices of the superstrate, waveguide, and substrate, respectively.

The above equations provide a rigorous description of RWGs. There is another approach providing the regions on the θ_inc − λ₀/Λ map where resonance may exist. To derive the respective formulas we review the physical properties of waveguides. A guided mode in a waveguide can only exist when:

1. Total reflection is satisfied: θ > θ_c^sup and θ > θ^sub_c , where θ^sup_c = sin⁻¹(n_sup/n_w)and θ_c^sub = sin⁻¹(n_sub/n_w). This also involves the condition: n_w >max{n_sup, n_sub}, 2. Transverse resonance condition is satisfied i.e. the repeatedly reflected wave has

constructive interference with itself. In this case only discrete angles: θ_m will be allowed in the range determined in the previous point. The longitudinal component of the wave vector equals β =k₀n_wsinθ_m.

The waveguide phase velocity is defined byv_p =ω/β, so that the effective refractive index can be defined as

n_eff =c/v_p =cβ/ω=β/k=n_wsinθ_m. (1.11)

The grating equation Eq. (1.6b) provides n_eff =n_incsinθ_inc−mλ₀

Λ. (1.12)

The conditions for n_eff to represent a guided mode are:

max{n_sup, n_sub} ≤n_eff < n_w, (1.13a) max{n_sup, n_sub} ≤

n_incsinθ_inc−mλ₀ d

< n_w. (1.13b)

As an example we analyze an infinite array of dielectric cylinders, which can be considered as an RWG structure. This structure will later be used in this Thesis to demon- strate the simulation of second harmonic generation by the NonLinear Finite Difference Frequency Domain (NL-FDFD) method in Sec.4so it is worth to examine itsRWGprop- erties here. The normalized wavelength is selected to be λ₀/Λ = 0.5. . .2. The period of the structure is Λ, and the radius of the cylinders is r = 0.1Λ, the refractive index of the cylinders is n=√

2. First we determine the effective refractive index:

n_eff = nr²π+ (Λ² −r²π)

/Λ² = 0.01πn+ (1−0.01π) = 1.0651. (1.14)

(26)

We determine the RWG map in the λ₀/Λ −θ plane indicating the contribution of the different modes. We apply two methods: the rigorous method (Eq. (1.9)), which provides

“exact” curves and the previous approximation (Eq. (1.13b)) which provides valid regions for the resonant process.

Fig. 1.5. RWG resonance on theθ_inc−λ₀/Λ map calculated with both methods.

Fig.1.5 showsRWGresonance on the θ_inc−λ₀/Λmap calculated with both methods.

The source code of the program resulting this figure is listed in Appendix A as PRG Fig. 1.5. The different diffracted orders are marked by distinctive colors, the results of the rigorous method (Eq. (1.9)) are denoted by solid lines and the regions resulting from Eq. (1.13b) by dashed lines. Determination of the dashed lines from Eq. (1.13b) is direct but the solid lines (Eq. (1.9)) require the application of a nonlinear solver. The simulation result shows that for the selected angle of incidence (θ_inc = 8.25^◦) only three modes are active in the given λ₀/Λ region. The exact values of the normalized wavelength of the first and second order diffraction orders at the selected angle are 0.593 and 1.184, their ratio is1.997. If the cylinders are made from second order nonlinear crystals the structure is resonant for both the fundamental and generated second harmonic waves at the same time. That is why such a structure provides high efficient second harmonic generation and serves as a good tool for the demonstration of numerical methods simulating second harmonic generation. Reflection peaks of Fig. 4.7 show good correlation with resonance positions at the selected angle of incidence in Fig. 1.5 (red line) taking into account the

(27)

fact that the definition of the x axis is different (λ₀/Λ vs a/λ₀, where a is the grating period (= Λ) but with different notation.). The reason of using a different notation in Fig. 4.7 is to correlate the results with former ones in Ref. [31].

1.2 Nonlinear structured dielectrics

The nonlinear response of a dielectric medium to an incoming electromagnetic wave is included into Maxwell’s equations via the electric polarization vector (the physical background of the effect and the causes of the nonlinear response can be found in Refs.

[32–34]), the series expansion of the polarization P in terms of the powers of the electric field E is given by

P(t) = ε₀

χ⁽¹⁾E(t) +χ⁽²⁾E²(t) +χ⁽³⁾E³(t) +. . .

, (1.15)

here field values are considered to be scalar for simplicity (as in Eq. (1.1.2) in [33]). The term χ⁽¹⁾ describes the linear response, the rest of the terms are the nonlinear effects.

Nonlinear optics can be divided into three groups according to the terms above. The higher harmonics are represented byχ⁽²⁾,χ⁽³⁾, and higher terms gaining large importance recently in the physical processes induced by ultra short and high energy pulsed lasers (ELI [35]). In most third order nonlinear materials (termχ⁽³⁾), the tensor character of the third order susceptibility is not manifested and reduces to a scalar, i.e. a constant times the unit tensor (isotropy of the medium being a necessary but not sufficient condition thereof).

Nonlinearity results in an intensity dependent refractive index: the associated nonlinear process is the optical Kerr effect [32–34]. In this Thesis we focus on the simulation of the processes associated with the second order term (χ⁽²⁾), and so only this term will be considered in the nonlinear polarization from now on.

Materials with second order nonlinearity cannot be centrosymmetric, more precisely, they have to be polar [32–34]. An important process resulting from second order nonlinearity is Second Harmonic Generation (SHG) in nonlinear optical crystals. Efficient generation of a Second Harmonic Wave (SHW) requires phase matching between the fundamental and generated waves [32–34]. In the bulk of unstructured crystals phase matching can be achieved by appropriately choosing the geometry (direction of propagation and polarization with respect to the crystal optical axes) for the Fundamental Wave (FW). Another way of realizing phase matching is modulating spatially the second order susceptibility of the nonlinear crystal; this leads to Quasi Phase Matching (QPM) that can be achieved e.g. by periodic poling of the second order susceptibility [33,34,36,37].

The interaction between the fundamental and generated waves is further enhanced in 1D

(28)

multilayer (i.e. 1D NonLinear Photonic Crystal (NLPC)) structures [38–40]. Based on novel technological developments, periodic poling can be done in a 2D pattern, which is exploited in various nonlinear diffraction processes [41–45].

Maxwell’s equations in the time domain read

∇ ×H=∂_tε₀E+∂_tP, (1.16a)

∇ ×E=−∂_tµ₀H. (1.16b)

The polarization vector P consists of a linear and a nonlinear part: P=P_L+P_{N L}. The linear part PL = ε0χE can be merged with the ε0E term in Eq. (1.16a), which yields ε0(I+χ)E, where χ is the two-index linear susceptibility tensor. All vector and tensor quantities have to be expressed in the principal optical axis coordinate system of the crystal. In this coordinate system the electric displacement vector, which describes the linear response is given by

D =ε0(I+χ)E=ε0







εxx 0 0 0 ε_yy 0 0 0 ε_zz











 Ex

E_y E_z







, (1.17)

where the relative susceptibility tensor is introduced with diagonal components ε_ii. The nonlinear polarization serves as a source term in the frequency domain Maxwell’s equations, from which the NL-FDFDmethod is derived. Assuming that the incoming electric field consists of a superposition of narrow band fields E(ω_n), the i-th component of the second order nonlinear polarization can be written as (see Eq. (1.5.21) in [33])

P_i^{N L}(ω_n+ω_m) = ε₀X

nm

d_ijkE_j^ωⁿE_k^ω^m, (1.18)

where d_ijk is the nonlinear susceptibility tensor (see Eq. (1.5.20) and Eq. (1.1.2) in [33]), and the Einstein convention is used for the summation of spatial indices. There is a factor of 2 difference between Eq. (1.18) and Eq. (1.5.21) in [33], which is caused by a different definition of complex field components. The definition conform to linear FDFD and adopted by us reads:

E(t) = X

n

< E^ωⁿe^jωⁿ^t

=X

n

1

2 E^ωⁿe^jωⁿ^t+c.c.

, (1.19)

for more details see Sec. 1.2.2 in [33] and the comment/footnote after Eq. (1.2.8) at the bottom of page 7 in that Ref..

(29)

In the cases of second harmonic generation (ω, ω → 2ω) and difference frequency generation (2ω, ω →ω) Eq. (1.18) leads to the following formulas, respectively (see Eqs.

(1.5.26) and (1.5.27) in [33] applying the field definition above):

P_i^{N L,2ω} =ε₀d_ijkE_j^ωE_k^ω, (1.20a)

P_i^{N L,ω} = 2ε0dijkE_j^2ω(E_k^ω)^∗ . (1.20b)

A difference of factor 2 between Eqs. (1.20a) and (1.20b) is due to the degeneracy (see Eq. (55) in [34]).

The nonlinear polarization in Eq. (1.20a) is the origin of the SHG process, which indicates the appearing and building up of the 2ω field. The nonlinear polarization in Eq. (1.20b) describes the back action of the generated SHW field on the fundamental (FW) field, which is manifested as difference frequency generation between the2ω and ω fields. In the undepleted pump approximation the ω field is considered constant, and the generation process (Eq. (1.20a)) is regarded. If pump depletion is taken into account, this is no longer the case and both equations, Eqs. (1.20a) and (1.20b) are considered.

1.2.1 Homogeneous crystal

As a first example we examine a homogeneous crystal unstructured both in its linear susceptibility (no micromachining) and in its nonlinear susceptibility (no poling technique). There is no phase matching, the process of SHG is examined in this condition.

For second harmonic generation the conversion efficiency is defined by η_ω = I_ω

I₀, η_2ω = I_2ω

I₀ , (1.21)

where I₀ = 1

2ε₀n_ωcE₀², (1.22)

is the intensity of incoming light, and in general

I_ω =|S_ω|, S_ω = ¹₂<(E_ω×H^∗_ω), (1.23) I_2ω =|S_2ω|, S_2ω = ¹₂<(E_2ω×H^∗_2ω), (1.24) where Z₀ =q

µ0

ε0 is the vacuum impedance.

For the case of two-dimensional TE_z, if the generated SHW is not collinear with the

(30)

FW, the components of the Poynting vector are:

S_x^2ω = 1

2Z₀<(E_z^ω(jH˜_y^ω)^∗), S_y^2ω = 1

2Z₀<(E_z^ω(jH˜_x^ω)^∗), (1.25) the intensity and the angle of energy flow are:

I_2ω = q

(S_x^2ω)²+ S_y^2ω2

, ](S^2ω) = arctan S_y^2ω

S_x^2ω

. (1.26)

Second harmonic generation at low conversion efficiency can be easily calculated analytically. Using Table 3 on page 44 of [34] the analytical solution for the conversion efficiency is:

η_2ω = y

L_{N L} 2

sinc²(∆ky/2) =

2

∆kL_{N L} 2

sin²(∆ky/2), (1.27) where L_{N L} is the nonlinear length:

L_{N L} = 1 4πdeff

s

2ε₀n²_ωn_2ωcλ²₀ I0

= λ₀√ n_ωn_2ω 2πdeffE0

, (1.28)

where I₀ is substituted from Eq. (1.22) and λ₀ is the wavelength of the FW in vacuum.

Using our definition for SHG, see Eqs. (1.20a) and (1.20b): d_eff = ¹₂χ⁽²⁾. The phase mismatch, ∆k is detailed in the next section. Neither pump depletion nor backward SHW generation is taken into account in this formula.

The amplitude of the alternation of the conversion efficiency depends on the phase mismatch and the intensity ofFW. For real materials and laser fields the phase mismatch is typically in the range of low conversion efficiencies so Eq. (1.27) is valid. This formula describes perfectly the SHG process in a non-phase-matched homogeneous crystal.

1.2.2 Phase Matching

Due to the linear dispersion of the underlying dielectric medium, the propagation velocities of the FW and SHW are different. This results in a continuous phase slip between FW and SHW in the course of propagation. Accordingly the newly generated SHWwill periodically constructively or destructively interfere with the pre-existingSHW.

The period of the alternation is twice the coherence length. The coherence length of the SHG process is given by [34], Eq. (41) page 44:

l_c = π

|∆k|, ∆k =k_2ω−2k_ω, k_2ω = 2π

λ0/2n_2ω, 2k_ω = 22π λ0

n_ω. (1.29)

(31)

The standard way to obtain phase matching for collinear beams is using birefringent optical crystals. There are several ways to achieve phase matching by exploiting the difference in the refractive indices for ordinary and extraordinary polarizations. For more detail see Tab. 2.3.2 in [33]. Here we consider as an example type I phase matching forSHG in a negative uniaxial crystal. Fig. 1.6 shows the dependence of ordinary and extraordinary refractive indices on the wavelength in such a crystal. Fig. 1.7 illustrates the process.

Fig. 1.6. Dependence of ordinary and extraordinary refractive indices on the wavelength.

The original figure can be found in [46].

Fig. 1.7. Geometry of type I SHG in a negative uniaxial crystal. Polarization states of the fields are shown, and cdenotes the orientation of the crystal.

The criterion for phase matching is

n_e(2ω, θ) =n_o(ω), (1.30)

where ne is the extraordinary refractive index as a function of the frequency and the angle of the crystal axis, while the ordinary refractive index n_o is only a function of frequency. Whether a wave propagates as an ordinary or extraordinary wave in a crystal is determined by its polarization state. In a birefringent crystaln_e(θ)can be obtained (see

(32)

Eq. (2.3.8) in [33]) by 1

ne(θ)² = sin²θ

n²_e +cos²θ

n²_o . (1.31)

Using Eqs. (1.30) and (1.31) the criterion for phase matching is sin²θ

n_e(2ω)² + cos²θ

n_o(2ω)² = 1

n_o(ω)². (1.32)

Considering dispersion relations (see Fig.1.6) for ordinary and extraordinary refractive indices one can find the angle satisfying the phase matching condition. There are some remarks for this type of phase matching:

1. It is very sensitive for the angle, so the orientation should be kept stable and in long crystals only parallel, slightly focused beams can be applied.

2. In a birefringent crystal the Poynting vectorSand the wave vectorkare not parallel, so ordinary and extraordinary rays with parallel wave vectors quickly diverge from one another. This is the so-called spatial walk-off effect, which limits the spatial overlapping of the generating fields.

3. It only works with certain components of the d_ijk tensor, which are sometimes not the highest ones. (This limitation will be circumvented by quasi phase matching.) As a final remark, the phase matching described above is usually called “angle tuning”, another option is “temperature tuning”, which eliminates the walk-off effect of the angle- tuned version.

1.2.3 Quasi Phase Matching

Another way to overcome the phase mismatch problem is the so-called Quasi Phase Matching (QPM) based on the periodical modulation of the sign of a selecteddijk component. In practice this can be achieved by the poling technique [33,34,36,37]. The method works in the following way: as seen in homogeneous crystals there is a constructive interference within the coherence length between the newly generated SHWand the pre-existing SHW which becomes destructive beyond the coherence length, leading to smaller SHW amplitude. By switching the sign of the d_ijk component at the coherence boundary the interference continues to be constructive.

As a basic example we consider the Periodically Poled Lithium Niobate (PPLN) crystal, where phase matching takes place collinearly. In an ideal PPLN material the periodicity of the poling is Λ = 2l_c, with 50% duty cycle, i.e. the length of each stack is equal

(33)

to the coherence length, which is in this case:

lc = π

|∆k| = λ₀

4(n^e_2ω−n^e_ω). (1.33)

Although this formula is derived for the negligible pump depletion case, it is also valid for pump depletion (see top of page 89 in [34]). Fig.1.8 shows the layout of such a structure.

Fig. 1.8. PPLN crystal for collinear quasi phase matching. Arrows indicate domain inversion.

The conversion efficiency is shown simultaneously in Fig. 1.9 in the homogeneous Lithium Niobate (LN) slab and in the PPLN slab, for the same input FW field strength.

Constructive interference of SHG can be seen as periodic poling is applied in accordance with the coherence length. In the homogeneous case (green curve) after one coherence length of propagation, the second harmonic field starts decreasing until it vanishes totally, and then the build up and drop down stages are periodically repeated. In the PPLNcase the sign of the nonlinear coefficient d33 gets inverted when the drop down is about to start, thus the second harmonic fields continue building up (blue curve).

If phase matching is applied the negligible pump depletion approximation is not valid anymore, so Eq. (1.27) cannot be applied for high conversion efficiencies (pump depletion).

In other words for negligible pump depletion only Eq. (1.20a) has to be considered, while in the case of pump depletion Eq. (1.20b) as well. The conversion efficiency for pump depletion can be calculated analytically (see [34], table 4, page 47):

η_2ω = tanh²(y/L_{N L}), (1.34)

In the case of QPM the effective nonlinear coefficient (dQ) is smaller than for perfect phase matching, the exact value is calculated as (see [33], Eq. (2.4.7), page 87):

d_Q = 2

πd_eff. (1.35)

(34)

Fig. 1.9. Conversion efficiency for PPLN crystal and homogeneous crystal.

So in Eq. (1.28)d_eff has to be replaced byd_Q = _π²d_eff = _π¹χ⁽²⁾, remarking the fact that the resulting “theoretical” curve will be just a smoothed approximation of the real conversion efficiency in the PPLNcase. In Fig.1.9 the red curve indicates the “theoretical” result for perfect phase matching when applying the effective value for d33.

The PPLN structure described above is the so-called first order QPMstructure since the change of the sign of the nonlinear coefficient occurs each time after a coherence length. Less efficient phase matching can be obtained when sign change occurs after an odd multiple of the coherence length, these types are called third order, etc. phase matching. Intensity amplification for different order QPMstructures and the ideal phase matching can be seen on Fig.1.10. This figure also indicates that even if we forced keeping constructive interference by sign change, the quasi phase matching cannot be as effective as “normal” phase matching if we use the same d_ijk component. But the big advantage of QPM is the freedom of selection of the d_ijk component, which can be then chosen to be the maximal one. Therefore, higher amplification can be reached for certain crystals by QPM than by normal phase matching.

If the second order dispersion is spatially structured like a grating, the phase matching can be realized in a non-collinear way. The phase matching condition is k_g =k_2ω−2k_ω. In Fig. 1.11, the dashed line is parallel to the interface between adjacent oppositely poled regions, k_g is perpendicular to them. The wave vector of the FW, k_ω makes an angle

(35)

Fig. 1.10. Intensity amplification for different orderQPMstructures and the ideal phase matching. The original figure can be found in [34] as Fig. 25.

φ with the grating, while the angle between k_2ω and k_ω is θ. For simplicity we neglect dispersion. In this case the wave vectors of the FW and SHW will point to symmetrical directions with respect to the grating, i. e. θ = 2φ, and the phase matching condition for the grating is

d= λ₀

4nsinφ, (1.36)

where λ₀ is the wavelength of the FW in vacuum, n is the refractive index of the crystal, and d is the period of modulation.

Fig. 1.11. Phase matching in tilted QPMgrating.

Even if (quasi) phase matching is applied, but either the intensity of incoming FW is small or the propagation length is short, the process can be considered to be in the low conversion regime, so pump depletion can be ignored. Thus the analytical conversion effi-

(36)

ciency can be obtained by considering the first term in the Taylor series of the conversion efficiency in Eq. (1.34) (see in [33]):

η_2ω = y

L_{N L} 2

, (1.37)

or it can be obtained simply by taking the limit of the sinc function in Eq. (1.27) when

∆k goes to zero.

(37)

Chapter 2 Simulational methods

Computational ElectroMagnetic (CEM) methods can be classified in several ways [47]. A summary chart can be seen in Fig. 2.1, which categorizes most types of CEM methods. High frequency methods are implemented by most commercial optical design softwares (e.g. Zemax [48]): standard ray tracing and diffraction simulation using either the Fourier or Huygens integral method are well-known examples. These methods can only be applied if structure parameters are much larger than the wavelength (d λ).

When structure parameters are comparable to the wavelength (d ∼λ), “true” numerical methods are needed, where the polarization and the vector nature of the ElectroMagnetic (EM) fields become important. Two methods are usually distinguished. The plane wave expansion method [49–54] describes the band structure of the photonic crystal, and the Transfer Matrix Method (TMM) [55–61] calculates the transmission and reflection of a one-dimensional layered structure.

Numerical methods can be classified further based on the criterion whether the differential or integral form of Maxwell’s equations are discretized. The most popular methods derived from the integral form are the Method of Moments [62–64], and the Fast Multipole Method [65,66]. Methods derived from the differential form require a grid on which theEM field values are defined. Depending on whether the grid is structured or non-structured the methods can be further classified. If the geometry of the structure is complex and contains many curved surfaces the usage of a non-structured grid seems to be the right choice, this method is called the Finite Element Method (FEM) [67–69], which is widely used in commercial softwares.

The methods, which will be used in this work belong to the structured grid category. Two main classes are included depending on whether the method is defined in the time or in the frequency domain. The Finite Difference Time Domain (FDTD) method (see Sec. 2.1) is very popular recently, due to its versatile application for wide range of

(38)

problems, like active devices, and much more. The PseudoSpectral Time Domain (PSTD) method [70–73] is an advanced form ofFDTDwhere the discrete differential operators are replaced by derivation using Fourier transformations. In the frequency domain the methods can be further divided into real space and Fourier space type. Among the real space methods we consider first the Finite Difference Frequency Domain (FDFD) method (see Sec.2.2). It starts from Maxwell’s equations in the frequency space, but only one frequency is regarded (monochromatic simulation). By rearranging the two- or three-dimensional field matrices into a vector, the final field distribution is obtained via a matrix inversion.

The Beam Propagation Method (BPM) [74–76] is basically a technique for simulating the propagation of light in slowly varying approximation, it is a one-way model, multiple reflections (like in fiber Bragg grating structures) are eliminated, however, a stable bidi- rectional beam propagation method has already been worked out [77,78]. The Method of Lines (MoL) (see Sec. 2.3) leaves one axis analytical and discretization is present only in the transversal plane. Accordingly the structure is divided into “lines”, fields are matched between the layers, reflection and transmission are calculated by connecting the first layer to the last one. Finally the Rigorously Coupled Wave Analysis (RCWA) (see Sec. 2.4)

(39)

is mentioned, which operates on the Fourier components of the fields. This method is widely used for optimization of subwavelength structures both in academic and industrial applications (optical critical dimension calculations).

The various methods are optimal for different problems. Fig. 2.2 details the strength and weaknesses of the four selected methods (FDTD, FDFD, MoL, RCWA) indicating their optimal applications. This summary table is based on lecture notes in [79]. The FDTD method is an early discretized representation of Maxwell’s equations, hence its presentation embeds the concepts for the other methods. Necessary tools for numerical electromagnetic simulations such as the implementation of the source and the absorbing boundary condition are introduced by the FDTD method. The FDFD method is in the main focus of this Thesis since this method will be extended for second harmonic nonlinear simulations. The FDFD method uses the same formulation for spatial derivatives as FDTD, the main difference is that the FDFD method operates in the frequency domain.

MoL plays the role of a natural bridge between the FDFD and RCWA methods. It will be detailed not only for the sake of completeness but MoL also fits better for metals or structures complex in longitudinal dimension than the RCWA method. The RCWA method is used for the optimization of structures in Chapter 3. It is very fast and used in commercial environments as well.

Fig. 2.2. Advantages and disadvantages of the four selected methods.

FDTD and RCWA are the most used and widespread methods on structured grids in CEM. Fig. 2.3 shows the step-by-step evolution of the FDTD method leading to the RCWA method, giving a comprehensive picture of numerical simulational methods. One of the reasons of dealing with four different methods is the possibility of direct comparison

(40)

of the results, another is finding - by scrutinizing their implementation - the method best fitted for generalization in second order nonlinear simulations.

Fig. 2.3. Route fromFDTD to RCWA via FDFDand MoL methods.

In the rest of this chapter we elaborate the details of these methods. All of them will be used for simulation of linear dielectric structures in Chapter 3, where the RCWA method will also be used for the optimization of the structures. The FDFD method will be developed further for second harmonic nonlinearity in Chapter4. In Chapter5a novel method called PSFDwill be detailed, which is based on the FDFD method.

2.1 FDTD method

In the FDTD method Maxwell’s equations are directly converted from analytical formulas to numerical ones: from a system of differential equations to a system of finite difference equations. The discretization can be carried out in several ways, we choose one preferred way to do it using a very stable and popular approach, the Yee algorithm [80,81].

2.1.1 Discretization of Maxwell’s equations, the Yee algorithm

As mentioned previously, simulation of Maxwell’s equations are performed in the FDTD method by the approximation of the infinitesimal partial derivation by discrete differentiation. As a first approach we consider a rectangular grid and all six EM field components (Ex, Ey, Ez, Hx, Hy, Hz) are defined in each grid point. By taking into account the structure and symmetry of Maxwell’s equations a much more efficient field arrangement can be achieved. The Yee grid (Yee mesh) realizes the optimal arrangement of fields, which fits best to the symmetry of Maxwell’s equations [80,81].

TamásSzarvas Simulationofelectromagneticwavepropagationinlinearandnonlinearstructureddielectrics Ph.D.Dissertation

Simulation of electromagnetic wave

propagation in linear and nonlinear structured dielectrics

Tamás Szarvas

Supervisor: Zsolt Kis, Ph.D.

Wigner Research Center for Physics Department of Quantum Optics and Quantum Informatics

Contents

List of Figures

List of Tables

Abbreviations

Notations

Chapter 1 Introduction

1.1 Linear structured dielectrics

1.1.1 Diffraction gratings

1.1.2 Surface relief gratings

1.1.3 Resonant waveguide gratings

1.2 Nonlinear structured dielectrics

1.2.1 Homogeneous crystal

1.2.2 Phase Matching

1.2.3 Quasi Phase Matching

Chapter 2

Simulational methods

2.1 FDTD method

2.1.1 Discretization of Maxwell’s equations, the Yee algorithm