Modeling of second-harmonic generation in periodic nanostructures by the Fourier modal method with matched coordinates

J. Defrance; M. Schäferling; T. Weiss

doi:10.1364/OE.26.013746

1. Introduction

In recent years, a major work has been devoted to utilizing resonant phenomena in nanophotonic systems [1–6]. This includes the control of far-field emission patterns [7–9] and the ability to achieve high intensities in the resonant near field [10–13], resulting in investigations of resonantly enhanced nonlinear effects in nanophotonics [14–19].

In order to understand the optical properties of nanophotonic systems and particularly for enhancing nonlinear effects, numerical methods play a crucial role. Different approaches have been developed to account for the nonlinear contributions to Maxwell’s equations in either time or frequency domain [20, 21]. In our paper, we focus on how to derive the second-harmonic emission from nanostructures with the Fourier modal method [22,23]. This numerical technique belongs to the class of modal methods, which are frequency domain solvers that are very efficient for calculating the optical properties of layered structures. The reason is that modal methods do not require a numerical discretization in the stacking direction of the layered structures. Instead, one searches for eigenmodes that can propagate or decay in the stacking direction. These eigenmodes serve as a basis for deriving the propagation of fields through a layer, which is usually carried out by the scattering matrix approach [24]. In the case of the Fourier modal method, the eigenmodes are calculated by decomposing Maxwell’s equations in a finite Fourier basis. This implies periodic boundary conditions, so that the Fourier modal method is ideal for layered periodic structures, i.e., for structures that are typically fabricated by lithographic methods.

Possible formulations for the calculation of harmonic generation in the Fourier modal method are described in [25,26]. They are based on the so-called undepleted pump approximation, where the energy depletion of the pump field is neglected in a first-order approximation, because the energy transfer from the pump field to any harmonic is rather small. As a consequence, the undepleted pump approximation allows for using a two-step process, in which we first derive the pump field from Maxwell’s equations for a given incident field. The resulting near fields are potential sources for fields at the harmonic. The emission intensity of these sources at the harmonic depends on the nonlinear susceptibility tensor as well as on the linear optical properties of the given system at the harmonic. Therefore, we need to solve only two linear problems – a homogeneous and an inhomogeneous differential equation for the pump field and the field generated at the harmonic, respectively. For the sake of simplicity, we focus here on second-harmonic generation, but the approach can be generalized for higher-harmonic generation.

While previous approaches for solving the inhomogeneous equations at the second harmonic rely on deriving the solutions for entire layers [25,26], we propose here to consider nonlinear sources in planes normal to the stacking direction and to coherently integrate over the emission of theses sources. Thus, we treat the second-harmonic emission in the same way as the emission from currents in certain layers [8,27], with the currents originating now in the nonlinear polarization. Moreover, we combine the derivation of the second-harmonic emission with two advanced formulations of the Fourier modal method: Factorization rules [22,28] and matched coordinates with adaptive spatial resolution [23, 29]. These formulations have been introduced in order to prevent inaccuracies in the Fourier modal method related to the Gibbs phenomenon that occur due to the finite Fourier basis. Regarding the factorization rules, we find that their formulation for the second-harmonic susceptibility described in previous works [25, 26] can be further improved. In particular, we show that it is necessary to account for the discontinuities of the fields at the pump energy and the second harmonic simultaneously. This new formulation of the Fourier factorization of the second-harmonic susceptibility tensor allows for including matched coordinates and adaptive spatial resolution. Thus, even nontrivial lateral geometries with large contrast in the dielectric function can be calculated efficiently by using a coordinate system in which all material interfaces are described by surfaces of constant coordinates with an increased spatial resolution in the vicinity of these interfaces.

The paper is structured in three main sections. The first one describes how to derive the second-harmonic emission of two-dimensional periodic nanostructures in a modal method. The second section addresses the Fourier modal method and contains our extended formulation of the factorization rules and the matched coordinates applied to the nonlinear susceptibility tensor. In the last part, we present numerical results.

2. Second-harmonic generation in modal methods

Modal methods are well-suited for multilayer systems that are composed of layers with a common translational symmetry. Modal methods benefit from this symmetry, because they are based on solving Maxwell’s equations separately in each layer. The solutions are then combined by considering the boundary conditions between adjacent layers.

For our derivations, we use the covariant formulation of Maxwell’s equations:

∊^{α β γ} \partial_{β} E_{γ} = i k_{0} μ^{α β} H_{β},

∊^{α β γ} \partial_{β} H_{γ} = - i k_{0} (ε^{α β} E_{β} + 4 π P^{α}) + \frac{4 π}{c} j^{α} .

Equations (1) and (2) hold for higher harmonics and for the pump field, with P^α and j^α denoting the α component of the nonlinear polarization and free currents, respectively. Furthermore, ∊^αβγ is the Levi-Civita symbol. The time dependence of the fields E and H is assumed to be exp(−iωt); the vacuum wave number is denoted as k₀ = ω/c. In the undepleted pump approximation for second-harmonic generation, P^α (ω) = 0 at the pump energy, while the simplest form at the second harmonic is given by [26]

P^{α} (2 ω) = χ^{(2), α β γ} (2 ω; ω, ω) E_{β} (ω) E_{γ} (ω),

with χ^(2),αβγ being the components of the second-harmonic susceptibility tensor. Therefore, in the undepleted pump approximation, the equations describing the pump field are independent of the fields at the second harmonic, whereas the pump field enters as an inhomogeneity in the equations at the second harmonic.

Let us first address the pump frequency. Without the loss of generality, we define x³ as the direction of translational symmetry (see Fig. 1). By eliminating the E₃ and H₃ components from Maxwell’s equations, we obtain the following system of equations [22]:

- i \partial_{3} F (x^{1}, x^{2}, x^{3}) = M (x^{1}, x^{2}) F (x^{1}, x^{2}, x^{3}) .

Here, F = (E₁, E₂, H₁, H₂)^T denotes a supervector of transversal field components. The Fourier transformed expression of the matrix operator M is given in Appendix A. Most importantly, it is independent of x³, so that we can make the ansatz F(x¹, x², x³) = F_m(x¹, x²) exp(iγ_mx³), which transforms Eq. (4) into the eigenvalue problem

γ_{m} F_{m} = M F_{m} .

with eigenmodes F_m that can propagate or decay with propagation constants γ_m in the direction of translational symmetry. The eigenmodes form a basis, in which we can expand arbitrary solutions of Maxwell’s equations inside that layer. In the following, we split the full set of eigenmodes {F_m} into subsets {

F_{m}^{\pm}

} with eigenvalues

γ_{m}^{\pm}

. The superscript + or − denotes the sign of the real part of

γ_{m}^{\pm}

(i.e., forward or backward propagation along the x³ direction) in the case that the imaginary part of

γ_{m}^{\pm}

vanishes. Otherwise,

sign (Im γ_{m}^{\pm}) = \pm 1

(i.e., forward or backward decay in the x³ direction). Thus, an arbitrary field can be expanded as

F (x^{1}, x^{2}, x^{3}) = \sum_{m} A_{m}^{+} F_{m}^{+} (x^{1}, x^{2}) e^{i γ_{m}^{+} (x^{3} - x_{+}^{3})} + A_{m}^{-} F_{m}^{-} (x^{1}, x^{2}) e^{i γ_{m}^{-} (x^{3} - x_{-}^{3})},

where

A_{m}^{\pm}

is the expansion coefficient calculated at positions

x_{\pm}^{3}

. The expansion coefficients

A_{m}^{\pm}

can be summarized as elements of supervectors

A_{l}^{\pm}

as the amplitude vectors for each layer l. The connection of the amplitude vectors at positions

x_{l}^{3}

and

x_{l^{'}}^{3}

(

x_{l}^{3} < x_{l^{'}}^{3}

) is given by the scattering matrix S_l′,l [24]:

\underset{\equiv A_{l^{'}, l}^{out}}{\underset{︸}{(\begin{matrix} A_{l'}^{+} \\ A_{l}^{-} \end{matrix})}} = \underset{\equiv S_{l^{'}, l}}{\underset{︸}{(\begin{matrix} S_{l^{'}, l}^{+ +} & S_{l^{'}, l}^{+ -} \\ S_{l^{'}, l}^{+ -} & S_{l^{'}, l}^{- -} \end{matrix})}} \underset{\equiv A_{l^{'}, l}^{in}}{\underset{︸}{(\begin{matrix} A_{l}^{+} \\ A_{l^{'}}^{-} \end{matrix})}} .

In the case that there is a free current located at the x¹x² plane through position $x_{0}^{3}$ ( $x_{l}^{3} < x_{0}^{3} < x_{l^{'}}^{3}$ ), such that $J (x^{3}) = J_{0} δ (x^{3} - x_{0}^{3})$ , it is possible to construct a similar matrix $Σ_{l^{'}, l} (x_{0}^{3})$ for the emission of that current based on [8,27]:

A_{l^{'}, l}^{out} = Σ_{l^{'}, l} (x_{0}^{3}) J_{0} .

The supervector J(x³) as well as Σ_l′,l(x³) are specified in Appendix B. Owing to the linearity of this equation, we may consider an arbitrary volumetric source J(x³) as a superposition of sources at positions x³, so that

A_{l^{'}, l}^{out} = \int_{x_{l}^{3}}^{x_{l^{'}}^{3}} d x^{3} Σ_{l^{'}, l} (x^{3}) J (x^{3}) .

From Eq. (2), it can be seen that we can treat the nonlinear polarization P^α in the same way as the current j^α, i.e., we define j^α = −iωP^α in the absence of free currents, and use Eq. (9) in order to calculate the emission of the nonlinear source to the far field.

Fig. 1 Lateral view of a layer of nonlinear material with circular air holes. This structure is one of the examples that will be treated in this paper. As illustrated, the linear field induces a nonlinear polarization that results in a volumetric source at the second harmonic, which we discretize in a set of planar emitting layers. For each discrete source at a position $x_{0}^{3}$ , we calculate the propagation of the field through the structure using the scattering matrix formalism. Then, we integrate coherently over all contributions to obtain the total far field. The inset in the right part depicts a schematic of the scattering matrix formalism, with the amplitude vectors A⁺ and A⁻ of downward and upward propagating or decaying eigenmodes in a certain layer [cf. equation (5)], respectively. The boxes with S_b and S_t indicate the scattering matrices of the sub-structures above and below a layer.

Download Full Size | PDF

3. Second-harmonic generation in the Fourier modal method

3.1. Fourier factorization rules

In the Fourier modal method, equation (5) is solved in a truncated Fourier space. However, it is known that the expansion of discontinuous functions in a truncated Fourier space exhibits over- and undershoots close to the jump discontinuity that are shifted to the discontinuity when increasing the dimension of the truncated Fourier space, with a pointwise convergence. This behavior, known as Gibbs phenomenon, hampers the performance of the Fourier modal method, because the dielectric function is discontinuous in layers that contain different materials. Even worse, whenever we multiply two discontinuous functions with concurrent jump discontinuities such as permittivity and normal component of the electric field at an interface between two materials, Lifeng Li has shown that the corresponding convolution in Fourier space does not converge at the concurrent jump discontinuity [28]. He has formulated so-called factorization rules, which allow to overcome this issue by reformulating Maxwell’s equations such that there are no products of functions with concurrent jump discontinuities [22]. These factorization rules have to be obeyed when transforming the material tensors to reciprocal space, including the second-order susceptibility tensor.

Let us consider an interface with a normal vector that equals the contravariant basis vector e^α. This means that the normal component of the dielectric displacement with respect to the interface is D^α, while E_β with β ≠ α gives the two tangential electric field components. The constitutive equation for D^α is as follows:

D^{α} = ε^{α α} E_{α} + ε^{α β} E_{β} + 4 π \underset{P^{α}}{\underset{︸}{χ^{(2), α σ ρ} {\hat{E}}_{σ} {\hat{E}}_{ρ}}}, β \neq α .

Note that the condition β ≠ α implies that there is no summation over indices α and the summation over β runs only over indices β ≠ α. For all other indices, we maintain the standard sum convention over all three components. The hat denotes quantities at the pump frequency.

In a first step, we express the discontinuous component of the electric field as a function of the other field components:

E_{α} = {(ε^{α α})}^{- 1} D^{α} - {(ε^{α α})}^{- 1} ε^{α β} E_{β} - {(ε^{α α})}^{- 1} 4 π χ^{(2), α σ ρ} {\hat{E}}_{σ} {\hat{E}}_{ρ}, β \neq α .

In this equation, the product of the linear permittivity and the field component E^α with a concurrent jump discontinuity has been removed, so that we can Fourier transform Eq. (11) in the x^α direction according to Li’s factorization rules. We denote this Fourier transformation by square brackets with subscript α. After reorganization of the terms, we obtain the Fourier transform of the electric displacement:

{[D^{α}]}_{α} = {[{(ε^{α α})}^{- 1}]}_{α}^{- 1} ({[E_{α}]}_{α} + {[{(ε^{α α})}^{- 1} ε^{α β}]}_{α} {[E_{β}]}_{α} + 4 π {[{(ε^{α α})}^{- 1} χ^{(2), α σ ρ} {\hat{E}}_{σ} {\hat{E}}_{ρ}]}_{α}), β \neq α .

This simple derivation provides the accurate Fourier transformation of D^α for e^α being normal to the interface. The Fourier transformation of the other two components D^β in direction x^α can be found by substituting E_α in the constitutive equation for D^β by Eq. (11):

\begin{array}{l} {[D^{β}]}_{α} & = & {[ε^{β α} {(ε^{α α})}^{- 1}]}_{α} {[D_{α}]}_{α} + {[ε^{β γ} - ε^{β α} {(ε^{α α})}^{- 1} ε^{α γ}]}_{α} {[E_{γ}]}_{α} β \neq α & γ \neq α \\ + 4 π {[χ^{(2), β σ ρ} - ε^{β α} {(ε^{α α})}^{- 1} χ^{(2), α σ ρ} {\hat{E}}_{σ} {\hat{E}}_{ρ}]}_{α} . \end{array}

Inserting Eq. (12) into Eq. (13) and repeating these steps for the second direction of periodicity yields the correctly factorized Fourier transform of the permittivity tensor. Following the notation of L. Li [22], we can summarize the above steps by defining an operator

l_{τ}^{+}

and its inverse

l_{τ}^{-}

, with

B \equiv l_{τ}^{\pm} A, B^{ρ σ} = {\begin{array}{l} {(A^{τ τ})}^{- 1}, & ρ = τ, σ = τ; \\ {(A^{τ τ})}^{- 1} A^{τ σ}, & ρ = τ, σ \neq τ; \\ A^{ρ τ} {(A^{τ τ})}^{- 1}, & ρ \neq τ, σ = τ; \\ A^{ρ σ} \pm A^{ρ τ} {(A^{τ τ})}^{- 1} A^{τ σ}, & ρ \neq τ, σ \neq τ, \end{array}

for arbitrary rank two tensors A and B. Thus, the correct application of the Fourier factorization rules to the permittivity tensor yields

\tilde{ε} = l_{2}^{+} F_{2} l_{2}^{-} l_{1}^{+} F_{l} l_{1}^{-} ε,

where F_τ is the Fourier transformation in the direction τ. The order of the Fourier transformation is arbitrary, but it should be noted that it has a small influence [22].

Similarly, it is possible to construct from Eqs. (12) and (13) an operator $m_{τ}^{\pm}$ that should be applied to the second-harmonic susceptibility tensor and accounts for the discontinuities at the second harmonic:

G \equiv m_{τ}^{\pm} (B) C, G^{α ρ σ} = {\begin{array}{l} B^{τ τ} C^{τ ρ σ}, & α = τ; \\ C^{α ρ σ} \pm B^{α τ} C^{τ ρ σ}, & α \neq τ . \end{array}

Here, the tensors B and C are results of applying specific operators to the permittivity and second-harmonic susceptibility, respectively. For the latter, we have to additionally take into account the discontinuities of the pump field. Hence, we express the discontinuous components of the electric field of the pump as a function of the continuous components of the electric field and displacement as well as the permittivity tensor evaluated at the pump energy:

{\hat{E}}_{α} = {({\hat{ε}}^{α α})}^{- 1} ({\hat{D}}^{α} - {\hat{ε}}^{α β} {\hat{E}}_{β}), β \neq α .

Inserting Eq. (17) into Eqs. (12) and (13) results in lengthy expressions that contain no products of functions with concurrent jump discontinuities, so that they can be Fourier transformed. They consist of products of $m_{τ}^{-} (l_{τ}^{-} ε)$ , ε̂ and χ⁽²⁾, which can be written in the compact form $m_{τ}^{-} (l_{τ}^{-} ε) j_{τ}^{-} (l_{τ}^{-} \hat{ε}) χ^{(2)}$ , with the new operator $j_{τ}^{\pm}$ being defined by

Γ \equiv j_{τ}^{\pm} (B) G, Γ^{α ρ σ} = {\begin{array}{l} G^{α τ τ} B^{τ τ} B^{τ τ}, & ρ = τ, σ = τ; \\ G^{α τ σ} B^{τ τ} \pm G^{α τ τ} B^{τ τ} B^{τ σ}, & ρ = τ, σ \neq τ; \\ G^{α ρ τ} B^{τ τ} \pm G^{α τ τ} B^{τ ρ} B^{τ τ}, & ρ \neq τ, σ = τ; \\ G^{α ρ σ} \pm G^{α τ σ} B^{τ ρ} \pm G^{α ρ τ} B^{τ σ} + G^{α τ τ} B^{τ ρ} B^{τ σ}, & ρ \neq τ, σ \neq τ . \end{array}

Note that

j_{τ}^{\pm} j_{τ}^{\mp} = m_{τ}^{\pm} m_{τ}^{\mp} = 1

, which allows us to recast Eqs. (12) and (13) into the form of Eq. (10) after carrying out the Fourier transformation. For brevity of notations, we define a new operator T_τ with

T_{τ} (ε, \hat{ε}) χ^{(2)} = j_{τ}^{+} (F_{τ} l_{τ}^{-} \hat{ε}) m_{τ}^{+} (F_{τ} l_{τ}^{-} ε) F_{τ} m_{τ}^{-} (l_{τ}^{-} ε) j_{τ}^{-} (l_{τ}^{-} \hat{ε}) χ^{(2)},

which gives the correctly factorized Fourier transform χ⁽²⁾ in direction x^τ. Thus, the final result for the correctly factorized Fourier transform of the second harmonic susceptibility yields

{\tilde{χ}}^{(2)} = T_{2} (l_{1}^{+} F_{1} l_{1}^{-} ε, l_{1}^{+} F_{1} l_{1}^{-} \hat{ε}) T_{1} (ε, \hat{ε}) χ^{(2)} .

This equation holds for the second-harmonic susceptibility tensor χ⁽²⁾, but it can be generalized to higher harmonics as well.

3.2. Adaptive spatial resolution and matched coordinates

While factorization rules allow for accurately describing products of functions with concurrent jump discontinuities in a finite Fourier space, they cannot overcome the Gibbs phenomenon in general. In particular, the Fourier expansion of functions with discontinuities exhibits over- and undershoots at the discontinuity that are approximately 18 % of the jump height [30]. This means that the dimension of the finite Fourier basis in the Fourier modal method needs to be large in order to achieve a convincing convergence behavior in the case that the index contrast in the dielectric function is high. This occurs, e.g., at interfaces between semiconductors and air, as well as at interfaces between dielectrics and metals. Therefore, Granet developed a new formulation of the Fourier modal method [29], in which he implemented coordinate transformations such that the spatial resolution near material interfaces is increased. The implementation of this adaptive spatial resolution requires replacing the original permittivity tensor by a redefined tensor that contains the discontinuous dielectric function times the metric components of the coordinate transformation. Thus, the influence of Gibbs phenomenon is drastically reduced, because the metric components are close to zero near the jump discontinuities, providing a smaller effective jump height in the adaptive coordinates.

However, the implementation of adaptive spatial resolution and factorization rules is only straight-forward in cases, where all material interfaces are aligned along the directions of periodicity [31,32]. The application of factorization rules in layers with complex cross sections requires to find at least a decomposition of the fields into normal and tangential components [33,34]. In the covariant formulation of Maxwell’s equations, normal and tangential field components are directly provided by the contra- and covariant field components in the case that all material interfaces are aligned along surfaces of constant coordinates [23]. This means that one needs to find matched coordinate systems that are adapted to the geometry of interest. Using the covariant formulation of Maxwell’s equations [Eqs. (1) and (2)], it is straight-forward to combine this approach with adaptive spatial resolution and factorization rules in order to achieve a fast convergence behavior even for metallo-dielectric systems [23,35,36]. An example of coordinates with adaptive spatial resolution matching a circular geometry is displayed in Fig. 2.

Fig. 2 Example of matched coordinates for the structure depicted in Fig. 1. Note that the matched coordinates include the material interfaces as surfaces of constant coordinates with an increased spatial resolution in their vicinity. The coordinate transformation is included in the numerical calculations in the form of redefined permittivity, permeability, and nonlinear susceptibility tensors.

Download Full Size | PDF

In the following, x̄^α and x^α denote the coordinates of uniform (Cartesian) and nonuniform (curvilinear) coordinate systems, respectively. The electric and magnetic field are covariant vectors of weight zero, whereas the electric displacement, the magnetic induction, and the nonlinear polarization are contravariant vectors of weight one. The permittivity and permeability are tensors of rank two and weight one and the second-harmonic susceptibility is a tensor of rank three and weight one. The transformation from uniform to curvilinear coordinates for two arbitrary tensors M and N of rank two and three, respectively, and weight one is given by

M^{ρ σ} = \frac{\partial (x^{1}, x^{2}, x^{3})}{\partial ({\bar{x}}^{1}, {\bar{x}}^{2}, {\bar{x}}^{3})} {\bar{M}}^{α β} \frac{\partial x^{ρ}}{\partial {\bar{x}}^{α}} \frac{\partial x^{σ}}{\partial {\bar{x}}^{β}}, N^{ρ σ τ} = \frac{\partial (x^{1}, x^{2}, x^{3})}{\partial ({\bar{x}}^{1}, {\bar{x}}^{2}, {\bar{x}}^{3})} {\bar{N}}^{α β γ} \frac{\partial x^{ρ}}{\partial {\bar{x}}^{α}} \frac{\partial x^{σ}}{\partial {\bar{x}}^{β}} \frac{\partial x^{τ}}{\partial {\bar{x}}^{γ}} .

In the Fourier modal method with matched coordinates, we now simply replace all quantities of the uniform coordinate system by those of the nonuniform coordinate system. Most importantly, ε, μ, and χ⁽²⁾ should be replaced by the redefined tensors given through Eq. (21) before carrying out their Fourier transformation according to the correct Fourier factorization rules given in section 3.1. The downside of the Fourier modal method with matched coordinates is that it requires solving homogeneous and isotropic layers numerically, even though knowing their analytical solutions in uniform coordinates. This additional effort is, however, usually compensated by the fast convergence of the results in layers with high index contrast, and it can be reduced by additional means [37].

4. Numerical results

In this section, we illustrate the performance of our method for two example geometries. The first geometry is an asymmetric one-dimensional grating made of gallium arsenide (GaAs). As this structure possesses only layers with trivial material distributions, we can easily compare our new formulation of the factorization rules for the second-harmonic susceptibility tensor with previous formulations of the factorization rules [25,26]. In the following, we refer to our formulation of the factorization rules as full factorization rules, whereas we denote previous formulations as simple factorization rules. In addition, we illustrate the impact of adaptive spatial resolution on the performance of the derivation of the second-harmonic emission for this simple example geometry. The second example geometry is a two-dimensional grating made of lithium niobate (LiNbO₃) with a circular cross section (see Fig. 1). This structure exhibits a high contrast of refractive index at its interfaces, and it possesses a nontrivial lateral cross section, so that previous formulations of the Fourier modal method cannot be used in order to derive the second-harmonic emission from this system accurately.

4.1. One-dimensional grating

A schematic of the asymmetric one-dimensional grating is depicted in Fig. 3(a). The structure consists of two adjacent periodic layers with GaAs surrounded by air. Each layer is 100 nm thick. The width of the GaAs region in the upper and lower layer is 0.66 μm and 0.75 μm, respectively. The period is 1 μm. Without the loss of generality, we consider a constant refractive index of n = 3.85 for GaAs.

Fig. 3 (a) One-dimensional periodic asymmetric grating made of GaAs and surrounded by air. This structure consists of two layers, where each layer has a of thickness of 100 nm and a period of 1 μm. The upper part of the grating has a width of 0.66 μm; the lower part has a width of 0.75 μm. The second-harmonic emission originates in the GaAs region. We consider normal incidence with the incident electric field aligned along the x¹ direction and an incoming power of 1W per unit cell area. (b) Linear transmission (red solid line), second-harmonic unpolarized intensity emitted above (black dashed line) and below (gray solid line) the structure as well as resonance energies (black arrows). (c) Relative error of the second-harmonic unpolarized intensity below the system at a pump energy of 670 meV for three different implementations: Simple factorization rules (dark blue line), full factorization rules with (light blue line) and without (blue line) adaptive spatial resolution.

Download Full Size | PDF

The crystallographic lattice of GaAs is face centered cubic; in Gaussian units, the reduced representation [38] of the second-harmonic susceptibility is

χ^{(2)} = 10^{- 7} (\begin{matrix} 0 & 0 & 0 & 5.73 & 0 & 0 \\ 0 & 0 & 0 & 0 & 5.73 & 0 \\ 0 & 0 & 0 & 0 & 0 & 5.73 \end{matrix}) \frac{cm}{statV} .

We consider an incident pump field that is normally incident with 1 W of incoming power per unit cell area. The incident electric field is linearly polarized and aligned along the direction of periodicity. The unpolarized second-harmonic intensity is calculated by summing up the s- and p-polarized intensities radiated in the fundamental diffraction orders.

In the numerical calculation of the second-harmonic emission, we use a spatial grid of 512 points per unit cell. Figure 3(b) shows the linear transmittance (red solid line) and the second-harmonic unpolarized intensity emitted above (black dashed line) and below (gray solid line) the structure calculated for a reciprocal basis (i.e., number of plane waves) of size 167. In order to illustrate the impact of resonances on the second-harmonic intensity, we denote the resonance energies of that structure by arrows. The resonance energies have been derived from the scattering matrix by the methods described in [37,39,40]. The fundamental resonance at 842 meV possesses a large linewidth of 84 meV, the higher-order resonances are narrower (from 4.4 meV to 60 meV linewidth). Evidently, we have a peak of emission whenever there is a resonance at the second-harmonic energy.

In order to illustrate the convergence of our calculations, we calculated the relative error of the second-harmonic unpolarized intensity for a successively increasing number of plane waves in the reciprocal basis. The relative error is calculated as the relative difference of the second-harmonic intensity for finite and infinite Fourier basis, where the latter has been derived from an exponential fit to the second-harmonic intensity as a function of the number of plane waves in the Fourier basis. The results derived by the Fourier modal method with simple (dark blue line) and full factorization rules with (light blue line) and without (blue line) adaptive spatial resolution are depicted in Fig. 3(c) at a pump energy of 670 meV. The formulation with full factorization rules converges significantly faster than the approach with simple factorization rules. The best convergence behavior can be obtained by the formulation with full factorization rules and adaptive spatial resolution.

4.2. Two-dimensional grating

The two-dimensional grating consists of a LiNbO₃ film with circular air holes (see Fig. 1). The period is 1.400 μm; the radius of the air holes is 600 nm; the thickness of the film is 80 nm. The incoming pump field is normally incident and linearly polarized with 1 W power per unit cell area. The incident electric field is aligned along the x¹ direction. The permittivity of LiNbO₃ is a diagonal tensor with ε¹¹ = ε²² = 2.29² and ε³³ = 2.21². The second-harmonic susceptibility is given by

χ^{(2)} = 10^{- 9} (\begin{matrix} 0 & 0 & 0 & 0 & - 13 & - 6.6 \\ - 6.6 & 6.6 & 0 & - 13 & 0 & 0 \\ - 13 & - 13 & - 75.9 & 0 & 0 & 0 \end{matrix}) \frac{cm}{statV} .

Figure 4(a) displays the linear transmittance and the second-harmonic unpolarized intensity of this structure calculated for 19 × 19 plane waves in the Fourier basis (quadratic truncation) and 512 × 512 spatial points. The structure exhibits a large number of resonances in the given energy range, but only few resonances can be excited by the pump or the second-harmonic energy, with resonance energies of 842 − 8.8i meV, 1152 − 25i meV, 1183 − 40i meV, 1700 − 46i meV and 2345 − 10i meV (the absolute value of the imaginary part is half of the resonance linewidth). Figure panels 4(b–c) contain the convergence behavior of the relative error calculated at pump energies of (b) 580 meV (dark blue line) and (c) 842 meV (blue line).

Fig. 4 (a) Second-harmonic unpolarized intensity above (black dashed line) and below (gray solid line) the structure, and linear transmittance spectra (red line) of the structure presented in Fig. 1. The thickness of the structure is 80 nm, the period is 1.400 μm along the x and y directions and the radius of the holes is 600 nm. (b–c) Convergence behavior of the calculated second-harmonic unpolarized intensity with the pump field at (a) 580 meV (dark blue line) and (b) 842 meV (blue line). The triangles in (a) mark the energies at which the convergence curves have been calculated.

Download Full Size | PDF

5. Conclusion

We have developed a formulation of the Fourier modal method for the calculation of the second-harmonic intensity in crossed gratings, in which we can combine factorization rules with matched coordinates and adaptive spatial resolution. In particular, our formulation of the factorization rules accounts for the permittivity tensor at the pump and the second-harmonic energy. Thus, we can calculate accurately the second-harmonic emission of periodic nanostructure arrays with complex cross sections, as demonstrated by our numerical results.

A. Basic principles of the Fourier modal method

In the Fourier modal method, we solve Maxwell’s equations in reciprocal space, with the electric and magnetic fields being decomposed in Fourier-Floquet series:

F_{α} (x^{1}, x^{2}, x^{3}) = \sum_{m, n} F_{α, m n} (x^{3}) e^{i k_{1, m} x^{1} + k_{2, n} x^{2}} .

Here, F_α denotes the vector components of either E_α or H_α for α ∈ {1, 2, 3} and m, n ∈ ℤ. Furthermore,

k_{α, m} = k_{α, 0} + G_{α, m}, G_{α, m} = m \frac{2 π}{P_{α}}, k_{1, 0} = N sin θ cos ϕ, k_{2, 0} = N sin θ sin ϕ,

with P_α as the period in direction α; N = 1 for the pump field and N = 2 at the second-harmonic energy [25]. The electric displacement and the second-harmonic polarization are given by

D_{m n}^{ρ} = \sum_{p, q} {\tilde{ε}}_{m n, p q}^{ρ σ} E_{σ, p q} + 4 π P_{m n}^{ρ},

P_{m n}^{ρ} = \sum_{p^{'}, q^{'}, p^{″}, q^{″}} {\tilde{χ}}_{m n, p^{'}, q^{'}, p^{″}, q^{″}}^{(2), ρ α β} {\hat{E}}_{α, p^{'}, q^{'}} {\hat{E}}_{β, p^{″}, q^{″}} .

Here, ε̃ and χ̃⁽²⁾ refer to the Fourier transform of the permittivity and the second-harmonic susceptibility carried out according to the factorization rules as defined in Eqs. (15) and (20).

The matrix M in Eq. (5) is given by:

\begin{array}{l} M & = & [\begin{array}{l} - {\overset{ˇ}{μ}}^{23} K_{2} - K_{1} {\overset{ˇ}{ε}}^{31} & {\overset{ˇ}{μ}}^{23} K_{1} - K_{1} {\overset{ˇ}{ε}}^{32} \\ {\overset{ˇ}{μ}}^{13} K_{2} - K_{2} {\overset{ˇ}{ε}}^{31} & - {\overset{ˇ}{μ}}^{13} K_{1} - K_{2} {\overset{ˇ}{ε}}^{32} \\ - k_{0} {\overset{ˇ}{ε}}^{21} - \frac{1}{k_{0}} K_{1} {\overset{ˇ}{μ}}^{33} K_{2} & - k_{0} {\overset{ˇ}{ε}}^{22} + \frac{1}{k_{0}} K_{1} {\overset{ˇ}{μ}}^{33} K_{1} \\ k_{0} {\overset{ˇ}{ε}}^{11} - \frac{1}{k_{0}} K_{2} {\overset{ˇ}{μ}}^{33} K_{2} & k_{0} {\overset{ˇ}{ε}}^{12} + \frac{1}{k_{0}} K_{2} {\overset{ˇ}{μ}}^{33} K_{1} \end{array} \\ \begin{array}{l} k_{0} {\overset{ˇ}{μ}}^{21} + \frac{1}{k_{0}} K_{1} {\overset{ˇ}{ε}}^{33} K_{2} & k_{0} {\overset{ˇ}{μ}}^{22} - \frac{1}{k_{0}} K_{1} {\overset{ˇ}{ε}}^{33} K_{1} \\ - k_{0} {\overset{ˇ}{μ}}^{11} + \frac{1}{k_{0}} K_{2} {\overset{ˇ}{ε}}^{33} K_{2} & - k_{0} {\overset{ˇ}{μ}}^{12} - \frac{1}{k_{0}} K_{2} {\overset{ˇ}{ε}}^{33} K_{1} \\ - {\overset{ˇ}{ε}}^{23} K_{2} - K_{1} {\overset{ˇ}{μ}}^{31} & {\overset{ˇ}{ε}}^{23} K_{1} - K_{1} {\overset{ˇ}{μ}}^{32} \\ {\overset{ˇ}{ε}}^{13} K_{2} - K_{2} {\overset{ˇ}{μ}}^{31} & - {\overset{ˇ}{ε}}^{13} K_{1} - K_{2} {\overset{ˇ}{μ}}^{32} \end{array}], \end{array}

where

\overset{ˇ}{ε} = l_{3}^{-} \tilde{ε}

and

\overset{ˇ}{μ} = l_{3}^{-} \tilde{μ}

[22]. Furthermore, K₁ and K₂ are diagonal matrices with diagonal elements k_1,m and k_2,m, respectively.

B. The scattering matrix formalism

The scattering matrix approach is a general formalism that provides the relations between the incoming and outgoing fields of a system. In a modal method, this system can be an interface, a layer or a system containing several layers and interfaces. In the case of the latter, the propagation through the layers and interfaces is calculated separately, and the results are combined recursively to obtain the scattering matrix of the multilayer. This procedure is performed in three steps [24].

In the first step, we build the propagation matrix P_l(L) in layer l using the eigenvalues of this layer:

(\begin{matrix} A_{l}^{+} (x_{0}^{3} + L) \\ A_{l}^{-} (x_{0}^{3} + L) \end{matrix}) = \underset{P_{l} (L)}{\underset{︸}{(\begin{matrix} e^{i Γ_{l}^{+} L} & 0 \\ 0 & e^{i Γ_{l}^{-} L} \end{matrix})}} [\begin{matrix} A_{l}^{+} (x_{0}^{3}) \\ A_{l}^{-} (x_{0}^{3}) \end{matrix}] .

Here, L is the thickness of the layer, and P_l(L) is constituted by blocks of sub-matrices defined by the matrix exponentials of

Γ_{l}^{\pm}

, i.e., diagonal matrices containing the eigenvalues of layer l.

In the second step, we build the transfer matrix I_l,l−1 of an interface between the layer l − 1 and layer l using the boundary conditions of the tangential field components:

(\begin{matrix} A_{l}^{+} \\ A_{l}^{-} \end{matrix}) = \underset{I_{l, l - 1}}{\underset{︸}{F_{l}^{- 1} F_{l - 1}}} (\begin{matrix} A_{l - 1}^{+} \\ A_{l - 1}^{-} \end{matrix}) .

Here F_l and F_l−1 are the material matrices of layers l and l − 1, respectively. The columns of these matrices are the eigenvectors of the corresponding layers.

In the last step, we use an iterative procedure to combine subsequently the matrices of the interfaces between layers and the propagation inside each layer. The recursion starts with a unit scattering matrix. At each recursion, the current scattering matrix is combined with a matrix X that is either an interface matrix I_l,l−1 or an inverse propagation matrix P_l(−L), with

X \equiv (\begin{matrix} X^{+ +} & X^{+ -} \\ X^{- +} & X^{- -} \end{matrix}) .

Formally, we denote the combination of the scattering matrix S and matrix X by the operator *, which results in the scattering matrix for the next recursion

X * S = [\begin{matrix} - {(W^{+ +})}^{- 1} S^{+ +} & - {(W^{+ +})}^{- 1} W^{+ -} \\ S^{- +} - S^{- -} X^{- +} {(W^{+ +})}^{- 1} S^{+ +} & S^{- -} X^{- -} - S^{- -} X^{- +} {(W^{+ +})}^{- 1} W^{+ -} \end{matrix}] .

Here,

W^{+ \pm} = S^{+ -} X^{- \pm} - X^{+ \pm} .

Let us address now the second-harmonic emission. The corresponding source term in Eq. (8) can be derived from Eqs. (10) and (13) in [8], which yields

J (x^{3}) = - 4 π i F_{0}^{- 1} {\begin{matrix} - K_{1} {({\tilde{ε}}_{0}^{33})}^{- 1} P^{3} \\ - K_{2} {({\tilde{ε}}_{0}^{33})}^{- 1} P^{3} \\ k_{0} [{\tilde{ε}}_{0}^{23} {({\tilde{ε}}_{0}^{33})}^{- 1} P^{3} - P^{2}] \\ k_{0} [P^{1} - {\tilde{ε}}_{0}^{13} {({\tilde{ε}}_{0}^{33})}^{- 1} P^{3}] \end{matrix}},

with the eigenvector matrix F₀ and the Fourier transformed permittivity tensor ε̃₀ being that of the layer at x³. P^α is a supervector containing the amplitudes

P_{m n}^{α}

of the correctly Fourier-factorized second-harmonic polarization. Combining Eqs. (16) to (19) of [8], the scattering matrix for the emission of that source can be written as

Σ_{l^{'}, l} (x_{0}^{3}) = [\begin{matrix} S_{b}^{+ +} V e^{i Γ_{0}^{+} (L - x_{0}^{3})} & - S_{b}^{+ +} V e^{i Γ_{0}^{+} L} S_{t}^{+ -} e^{- i Γ_{0}^{-} x_{0}^{3}} \\ S_{t}^{- -} U e^{- i Γ_{0}^{-} L} S_{b}^{- +} e^{i Γ_{0}^{+} (L - x_{0}^{3})} & - S_{t}^{- -} U e^{- i Γ_{0}^{-} (x_{0}^{3})} \end{matrix}] .

Here, L is the thickness of the layer generating the second harmonic, S_t and S_b denote the matrices of the groups of layers above and below this region, respectively, and the exponential terms are diagonal matrices that account for the propagation inside the emitting layer. U and V are matrices that account for internal reflections:

U = {(𝟙 - e^{- i Γ_{0}^{-} L} S_{b}^{- +} e^{i Γ_{0}^{+} L} S_{t}^{+ -})}^{- 1}, V = {(𝟙 - e^{i Γ_{0}^{+} L} S_{t}^{+ -} e^{- i Γ_{0}^{-} L} S_{b}^{- +})}^{- 1} .

Funding

DFG SPP 1839; MWK Baden-Württemberg; VW foundation.

Acknowledgement

We acknowledge discussions with M. L. Nesterov, N. A. Gippius and S. G Tikhodeev on the proposed method for calculating the second harmonic emission by the Fourier modal method.

References and links

1. W. L. Barnes, “Surface plasmon-polariton length scales: a route to sub-wavelength optics,” J. Opt. A-Pure Appl. Op. 8, S87–S93 (2006). [CrossRef]

2. W. A. Murray and W. L. Barnes, “Plasmonic Materials,” Adv. Mater. 19, 3771–3782 (2007). [CrossRef]

3. M. Hentschel, M. Schäferling, X. Duan, H. Giessen, and N. Liu, “Chiral plasmonics,” Sci. Adv. 3, e1602735 (2017). [CrossRef] [PubMed]

4. A. Moreau, G. Granet, F. I. Baida, and D. V. Labeke, “Light transmission by subwavelength square coaxial aperture arrays in metallic films,” Opt. Express 11, 1131–1136 (2003). [CrossRef] [PubMed]

5. P. Gay-Balmaz and O. J. F. Martin, “Electromagnetic resonances in individual and coupled split-ring resonators,” Appl. Phys. 92, 2929–2936 (2002). [CrossRef]

6. N. P. de Leon, B. J. Shields, C. L. Yu, D. E. Englund, A. V. Akimov, M. D. Lukin, and H. Park, “Tailoring light-matter interaction with a nanoscale plasmon resonator,” Phys. Rev. Lett. 108, 226803 (2012). [CrossRef] [PubMed]

7. A. G. Curto, G. Volpe, T. H. Taminiau, M. P. Kreuzer, R. Quidant, and N. F. van Hulst, “Unidirectional emission of a quantum dot coupled to a nanoantenna,” Science 329, 930–933 (2010). [CrossRef] [PubMed]

8. S. V. Lobanov, T. Weiss, D. Dregely, H. Giessen, N. A. Gippius, and S. G. Tikhodeev, “Emission properties of an oscillating point dipole from a gold Yagi-Uda nanoantenna array,” Phys. Rev. B 85, 155137 (2012). [CrossRef]

9. D. Dregely, K. Lindfors, M. Lippitz, N. Engheta, M. Totzeck, and H. Giessen, “Imaging and steering an optical wireless nanoantenna link,” Nat. Commun. 2, 4354 (2014).

10. F. Neubrech, A. Pucci, T. W. Cornelius, S. Karim, A. García-Etxarri, and J. Aizpurua, “Resonant plasmonic and vibrational coupling in a tailored nanoantenna for infrared detection,” Phys. Rev. Lett. 101, 157403 (2008). [CrossRef] [PubMed]

11. M. L. Nesterov, X. Yin, M. Schäferling, H. Giessen, and T. Weiss, “The role of plasmon-generated near fields for enhanced circular dichroism spectroscopy,” ACS Photonics 4, 578–583 (2016). [CrossRef]

12. T. Schumacher, M. Brandstetter, D. Wolf, K. Kratzer, M. Hentschel, H. Giessen, and M. Lippitz, “The optimal antenna for nonlinear spectroscopy of weakly and strongly scattering nanoobjects,” Appl. Phys. B 122, 91 (2016). [CrossRef]

13. K. Weber, M. L. Nesterov, T. Weiss, M. Scherer, M. Hentschel, J. Vogt, C. Huck, W. Li, M. Dressel, H. Giessen, and F. Neubrech, “Wavelength dcaling in antenna-enhanced infrared spectroscopy: Towards the far-IR and THz region,” ACS Photonics 4, 45–51 (2017). [CrossRef]

14. B. Metzger, M. Hentschel, and H. Giessen, “Ultrafast nonlinear plasmonic spectroscopy: From dipole nanoantennas to complex hybrid plasmonic structures,” ACS Photonics 3, 1336–1350 (2016). [CrossRef]

15. G. Purvinis, P. S. Priambodo, M. Pomerantz, M. Zhou, T. A. Maldonado, and R. Magnusson, “Second-harmonic generation in resonant waveguide gratings incorporating ionic self-assembled monolayer polymer films,” Opt. Lett. 29, 1108–1110 (2004). [CrossRef] [PubMed]

16. J. A. H. van Nieuwstadt, M. Sandtke, R. H. Harmsen, F. B. Segerink, J. C. Prangsma, S. Enoch, and L. Kuipers, “Strong modification of the nonlinear optical response of metallic subwavelength hole arrays,” Phys. Rev. Lett. 97, 146102 (2006). [CrossRef] [PubMed]

17. M. Kauranen and A. V. Zayats, “Nonlinear plasmonics,” Nat. Photonics 6(11), 737–748 (2012). [CrossRef]

18. K. Thyagarajan, S. Rivier, A. Lovera, and O. J. F. Martin, “Enhanced second-harmonic generation from double resonant plasmonic antennae,” Opt. Express 20, 12860–12865 (2012). [CrossRef] [PubMed]

19. H. Aouani, M. Rahmani, M. Navarro-Cía, and S. A. Maier, “Third-harmonic-upconversion enhancement from a single semiconductor nanoparticle coupled to a plasmonic antenna,” Nat. Nanotechnol. 9, 290–294 (2014). [CrossRef] [PubMed]

20. D. M. Sullivan, “Nonlinear FDTD formulations using Z transforms,” IEEE T. Microw. Theory 43, 676–682 (1995). [CrossRef]

21. K. Hayata, M. Nagai, and M. Koshiba, “Finite-element formalism for nonlinear slab-guided waves,” IEEE T. Microw. Theory 36, 1207–1215 (1988). [CrossRef]

22. L. Li, “Fourier modal method for crossed anisotropic gratings with arbitrary permittivity and permeability tensors,” J. Opt. A-Pure Appl. Op. 5(4), 345 (2003). [CrossRef]

23. T. Weiss, G. Granet, N. A. Gippius, S. G. Tikhodeev, and H. Giessen, “Matched coordinates and adaptive spatial resolution in the Fourier modal method,” Opt. Express 17, 8051–8061 (2009). [CrossRef] [PubMed]

24. S. G. Tikhodeev, A. L. Yablonskii, E. A. Muljarov, N. A. Gippius, and T. Ishihara, “Quasiguided modes and optical properties of photonic crystal slabs,” Phys. Rev. B 66, 045102 (2002). [CrossRef]

25. B. Bai and J. Turunen, “Fourier modal method for the analysis of second-harmonic generation in two-dimensionally periodic structures containing anisotropic materials,” J. Opt. Soc. Am. B 24, 1105–1112 (2007). [CrossRef]

26. T. Paul, C. Rockstuhl, and F. Lederer, “A numerical approach for analyzing higher harmonic generation in multilayer nanostructures,” J. Opt. Soc. Am. B 27, 1118–1130 (2010). [CrossRef]

27. D. M. Whittaker and I. S. Culshaw, “Scattering-matrix treatment of patterned multilayer photonic structures,” Phys. Rev. B 60, 2610–2618 (1999). [CrossRef]

28. L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures,” J. Opt. Soc. Am. A 13, 1870–1876 (1996). [CrossRef]

29. G. Granet, “Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution,” J. Opt. Soc. Am. A 16, 2510–2516 (1999). [CrossRef]

30. A. Kufner and J. Kadlec, Fourier Series, 1st ed. (Iliffe Books, 1971).

31. G. Granet and J. P. Plumey, “Parametric formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. A-Pure Appl. Op. 4, S145–S149 (2002). [CrossRef]

32. T. Weiss, N. A. Gippius, S. G. Tikhodeev, G. Granet, and H. Giessen, “Efficient calculation of the optical properties of stacked metamaterials with a Fourier modal method,” J. Opt. A-Pure Appl. Op. 11, 114019 (2009). [CrossRef]

33. T. Schuster, J. Ruoff, N. Kerwien, S. Rafler, and W. Osten, “Normal vector method for convergence improvement using the RCWA for crossed gratings,” J. Opt. Soc. Am. A 24, 2880–2890 (2007). [CrossRef]

34. R. Antos, “Fourier factorization with complex polarization bases in modeling optics of discontinuous bi-periodic structures,” Opt. Express 17, 7269–7274 (2009). [CrossRef] [PubMed]

35. S. Essig and K. Busch, “Generation of adaptive coordinates and their use in the Fourier modal method,” Opt. Express 18, 23258–23274 (2010). [CrossRef] [PubMed]

36. J. Küchenmeister, T. Zebrowski, and K. Busch, “A construction guide to analytically generated meshes for the Fourier modal method,” Opt. Express 20, 17319–17347 (2012). [CrossRef] [PubMed]

37. T. Weiss, N. A. Gippius, S. G. Tikhodeev, G. Granet, and H. Giessen, “Derivation of plasmonic resonances in the Fourier modal method with adaptive spatial resolution and matched coordinates,” J. Opt. Soc. Am. A 28, 238–244 (2011). [CrossRef]

38. R. W. Boyd, Nonlinear optics, 3rd ed. (Academic Press, 2008).

39. T. Weiss, M. Mesch, M. Schäferling, H. Giessen, W. Langbein, and E. A. Muljarov, “From dark to bright: First-order perturbation theory with analytical mode normalization for plasmonic nanoantenna arrays applied to refractive index sensing,” Phys. Rev. Lett. 116, 237401 (2016). [CrossRef] [PubMed]

40. D. A. Bykov and L. L. Doskolovich, “Numerical Methods for Calculating Poles of the Scattering Matrix With Applications in Grating Theory,” J. Lightwave Technol. 31, 793–801 (2013). [CrossRef]

Modeling of second-harmonic generation in periodic nanostructures by the Fourier modal method with matched coordinates

Abstract

1. Introduction

2. Second-harmonic generation in modal methods

3. Second-harmonic generation in the Fourier modal method

3.1. Fourier factorization rules

3.2. Adaptive spatial resolution and matched coordinates

4. Numerical results

4.1. One-dimensional grating

4.2. Two-dimensional grating

5. Conclusion

A. Basic principles of the Fourier modal method

B. The scattering matrix formalism

Funding

Acknowledgement

References and links

Cited By

Figures (4)

Equations (36)

Optics Express