Materials and Methods
CHAPTER 4 MATERIALS AND METHODS
4.1.2 Photonic Crystal Simulations
The diffraction of light through PhCs cannot be calculated by the ray tracing used for the optical MC simulations. To evaluate the impact of PhCs on the light propagation, a Maxwell-solver is required that calculates the diffraction by the PhC. For this purpose this work uses GD-Calc software (KJ Innovation, Santa Clara, CA, USA) which is based on the rigorous coupled wave analysis [147]. The following sections are mainly based on Ref. [148].
4.1.2.1 Representation of Photonic Crystals in GD-Calc
In GD-Calc, a PhC slab is represented through a layer (called stratum) positioned between the superstrate from which light is incident (with RI nSup) and the substrate (nSub) as depicted in Fig. 4.3a. Within the PhC stratum, the grating geometry is modeled through cuboids, all of which have a constant RI. If the PhC contains round features, these have to be approximated in a staircase manner as depicted in Fig. 4.3b,c. A more complex geometry in which the PhC structures also vary with the z-position, such as conical-shaped cavities, requires the use of several layers. However, this work only considers PhC structures that do not depend on the z-position like cylinders. GD-Calc is only applicable to structures that are periodic within thex-y-plane and every grating is defined through a set of parameters as illustrated in Fig. 4.3 including
• primitive lattice vectors (⃗a1,⃗a2)
• lattice pitch a=|⃗a1|=|⃗a2| 54
4.1 SIMULATIONS
Superstrate
Substrate PhC stratum
x z
E ~i E ~r
E ~t a
y x y x
b c
a1
a2 ~
~ ~a2 ~a1
a h s
' '
Figure 4.3: a) Drawing of the simulated domain in GD-Calc. The arrows illustrate the directions of the EM plain waves corresponding to the electric fields E. The superscripts⃗ denote incident,reflected, andtransmitted waves. Perspective view of the cuboid represen-tation of an orthogonal (b) and a hexagonal (c) PhC grating within GD-Calc. Note that in (b) and (c), the number of cuboids approximating the rounded edges has been reduced and the filling material of the cavities as well as the superstrate are not shown for illustration purposes.
• cavity size s
• grating height h
• grating bulk RI nBulk
• cavity RI nCav
The grating symmetry is defined by the angle ϕ between ⃗a1 and ⃗a2. A value of ϕ = 90◦ yields an orthogonal grating and ϕ = 60◦ corresponds to a hexagonal grating.
4.1.2.2 Expansion of Electromagnetic Field and Permittivity
The electric field in a configuration as depicted in Fig. 4.3 can be expressed through a sum of the incident field E⃗i in the superstrate, the reflected field E⃗r in the superstrate, and the transmitted field E⃗t in the substrate,
E⃗ =E⃗i+E⃗r+E⃗t . (4.2)
In the following, the superscripts "i", "r", and "t" always connote quantities corresponding to incidence, reflection, and transmission, respectively. The incident field is given through a plane wave
E⃗i(⃗r) = E⃗0iexp(i⃗ki·⃗r), (4.3) with E⃗0i being a constant vector and⃗ki being the wave vector of the incident field. As shown in Sec. 3.2.2, thex-y-component of the total electric field is a periodic function in the presence of a PhC slab. Hence, E⃗ can be represented through a Fourier expansion in x and y,
E(⃗⃗ r) =
m1,m2
F2DE(⃗⃗ r, m1, m2) expi⃗kxy(m1, m2)·⃗rxy , (4.4)
CHAPTER 4 MATERIALS AND METHODS
where the 2D Fourier harmonics (m1, m2) correspond to the 2D diffraction orders of the PhC grating,⃗kxy and⃗rxy are thex-y-plane projections of⃗k and ⃗r, andF2DE(⃗⃗ r, m1, m2) is the 2D Fourier coefficient given through
F2DE(⃗⃗ r, m1, m2) =E⃗0(m1, m2) expi⃗k(m1, m2)·⃗r . (4.5) Generally, the diffraction orders m1 and m2 in the series of Eq. 4.4 range from −∞ to ∞.
For numerical computations, they must be limited to a finite number. In analogy to E, the⃗ total magnetic field H⃗ can be expressed through
H(⃗⃗ r) =
m1,m2
F2DH(⃗⃗ r, m1, m2) expi⃗kxy(m1, m2)·⃗rxy . (4.6)
Since the PhC structure in GD-Calc is modeled through a stack of layers that consist of cuboids (in our case a single layer), the permittivity within a stratum does not depend onz,
ε(⃗r) = ε(⃗rxy). (4.7)
As for the electric and magnetic fields, the periodicity of the PhC allows representing the permittivity within the grating through a Fourier series,
ε(⃗rxy) =
l1,l2
F2Dε(l1, l2) expi (l1⃗b1+l2⃗b2)·⃗rxy , (4.8) using the Fourier order indicesl1 and l2 as well as the primitive reciprocal lattice vectors⃗b1
and⃗b2. The in-plane wave vectors of the diffraction orders are correlated to the projection of the incident wave vector⃗kixy through
⃗kxy(m1, m2) =⃗kxyi +m1⃗b1+m2⃗b2 , (4.9) Based on Eq. 4.9, the third component kz(m1, m2) can be derived through the correlations between the wave vector lengths,
|⃗ki|=|⃗kr|, (4.10)
|⃗ki|= nSup
nSub|⃗kt|, (4.11)
where the first equation applies to reflected orders and the second applies to transmitted orders.
4.1.2.3 Calculation of Scattering Matrices
The propagation of the EM wave in the presence of the PhC slab is calculated by GD-Calc using scattering matrices (S-matrices), that represent a linear mapping function correlating the amplitudes and phases of an EM wave incident on one side of a stratum and outgoing from the opposite side. If the grating consists of more than one stratum, the individual S-matrices are stacked to obtain the total S-matrix.
56
4.1 SIMULATIONS
The S-matrix of a stratum is derived from the macroscopic Maxwell Eqs. 3.10-3.13 using the series representations forE,⃗ H, and⃗ εgiven in Eqs. 4.4, 4.6 and 4.8. The number of diffraction orders retained in the Fourier expansions (NOrders) have a significant impact on the accuracy and computation time of the algorithm. On the one hand, the stair-case approximation of the geometry caused by the block-wise partitioning requires sufficiently large values for NOrders to avoid spikes of the EM fields at the block corners. On the other hand, the calculation time scales approximately with NOrders3 . Thus, an initial convergence test was performed to evaluate a value for NOrders that balances accuracy and computation time.
4.1.2.4 Correlation between Incident and Diffracted Fields
Consider the incidence of a plane EM wave with arbitrary polarization vector ⃗q on the PhC slab. As introduced in Sec. 2.5.1, any polarization can be represented in the base (⃗s, ⃗p) of the linear S- and P-polarization states with complex-valued coefficients A and B,
⃗
q =A ⃗s+B ⃗p . (4.12)
The electric field of the incident wave is then given through
E⃗A,Bi (⃗r) =E0i(A ⃗s+B ⃗p) expi⃗ki·⃗r , (4.13) wherein E0i is a constant. For each diffraction order (m1, m2) retained in the Fourier series, the GD-Calc algorithm yields two S-matrices, one for reflection (R(m1, m2)) and one for transmission (T(m1, m2)). These are 2×2 matrices with complex-valued entries,
R(m1, m2) =
R and T correlate the incident and diffracted electric fields through linear transformations,
Here, the subscripts S and P connote the projections of the electric field onto the ⃗s and ⃗p polarization states. For the transformations throughRandT, the electric fields are evaluated at points ⃗r1 = (x1, y1, z1) and⃗r2 = (x2, y2, z2), which have identical in-plane coordinates,
x1 =x2 , (4.18)
y1 =y2 . (4.19)
CHAPTER 4 MATERIALS AND METHODS
Superstrate
Substrate PhC stratum
x z
ki
~
z2 z1
x1=x2
kr(m1,m2)
~
kt(m1,m2)
~
Figure 4.4: Illustration of wave vectors of incident, reflected, and transmitted field for a GD-Calc calculation showing the coordinates used for the transformations in Eqs. 4.16 and 4.17.
The vertical coordinatesz1 andz2indicate the superstrate- and substrate-oriented boundaries of the PhC slab (see Fig. 4.4), respectively.
The power of the incident and diffracted EM waves are proportional to the square of the electric field amplitude. Hence, Eqs. 4.16 and 4.17 can be used to calculate the fraction of incident power that has been diffracted into a certain order (m1, m2), i.e. the diffraction efficiency D(m1, m2). For the polarization state ⃗q defined through (A, B), the efficiencies of the reflected and transmitted orders are given through
DrA,B(m1, m2) = −|A RSS+B RSP|2+|A RPS+B RPP|2
|A|2+|B|2 · Re (krz)
Re (kiz) , (4.20) DtA,B(m1, m2) = −|A TSS+B TSP|2+|A TPS+B TPP|2
|A|2+|B|2 · Re (kzt)
Re (kzi) , (4.21) wherein Re denotes the real part of a complex number and kz are the z-components of the wave vectors which can be derived through Eqs. 4.10 and 4.11. Note that all matrix elements ofRandT as well askrz andktz implicitly depend on the diffraction order (m1, m2), which has been omitted for brevity. These correlations indicate that only propagating orders have non-zero diffraction efficiences, since evanescent orders have purely imaginary krz orkzt components. The overall reflection and transmission coefficients R and T of a PhC can be calculated by summing over all diffraction efficiencies,
RA,B=
m1,m2
DA,Br (m1, m2), (4.22)
TA,B=
m1,m2
DA,Bt (m1, m2). (4.23)
58
4.1 SIMULATIONS