Modal method for second harmonic generation in
nanostructures

Sébastien Héron, Fabrice Pardo, Patrick Bouchon, Jean-Luc Pelouard, Riad Haïdar

► To cite this version:

Sébastien Héron, Fabrice Pardo, Patrick Bouchon, Jean-Luc Pelouard, Riad Haïdar. Modal method for second harmonic generation in nanostructures. Journal of the Optical Society of America B, 2015, 32 (2), pp.275-280. 10.1364/JOSAB.32.000275 . hal-01221788

HAL Id: hal-01221788

https://hal.science/hal-01221788

Submitted on 5 Jul 2021

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Modal method for second harmonic
generation in nanostructures

Sébastien Héron,<sup>1,2</sup> Fabrice Pardo,<sup>2</sup> Patrick Bouchon,<sup>1,*</sup> Jean-Luc Pelouard,<sup>2</sup> and Riad Haïdar<sup>1,3</sup>

<sup>1</sup>ONERA—The French Aerospace Lab, F-91761 Palaiseau, France

<sup>2</sup>Laboratoire de Photonique et de Nanostructures (LPN-CNRS), Route de Nozay, 91460 Marcoussis, France

<sup>3</sup>École Polytechnique, Département de Physique, 91128 Palaiseau, France

*Corresponding author: patrick.bouchon@onera.fr

A modal method for the study of second harmonic (SH) generation in nanostructures is depicted both for plane waves and focused beam illumination. It is based on the B-spline modal method, and relies on the undepleted pump approximation. The nonlinear polarization induced by an incident plane wave or a focused beam generates a source term at the doubled frequency. The latter is divided into a finite number of sub-sources, and the SH field is subsequently computed by integration of these sub-source contributions.

OCIS codes: (000.4430) Numerical approximation and analysis; (050.1755) Computational electromagnetic methods; (050.1950) Diffraction gratings; (050.6624) Subwavelength structures; (190.2620) Harmonic generation and mixing.

# 1. INTRODUCTION

In the near- and the mid-infrared, nanostructured devices are giving promising results for controlling the localization of the electric field and for strong intensity field enhancement [1–4]. Since nonlinear optics processes benefit from a large enhancement of the local electric field, such nanostructures are very promising [5–8]. Their design and optimization need fast numerical simulation tools.

So far, various numerical methods have been presented to account for nonlinear processes in periodic nanostructures [9]. In particular, methods have been developed for computing the second harmonic (SH) generation based on the Green function method [10], the Dirichlet-to-Neumann map [11], the boundary element method [12], the surface integral method, [13] or the finite difference time-domain method [14]. In addition, modal methods are known to be efficient for the study of a stack of nanostructured layers, and the Fourier modal method has been adapted for the study of SH generation [15–17].

More recently, the B-spline modal method (BMM) was introduced to solve linear electromagnetic problems [18–21]. In this method, the eigenvalue equation describing the electromagnetic problem is sparse, allowing the calculation speed to increase by solving the problem for only a few eigenmodes. Additionally, the B-splines can be nonuniform which allows for a better description of permittivity breaks. In this article, we adapt the BMM to map the SH field throughout periodic nanostructured layers for plane waves (of arbitrary incident angle) under the undepleted pump approximation. In addition, to simulate experimental setups with a focused incident light, we develop a computation of SH field for a cylindrical focused beam at the fundamental frequency. In Section 2, the consecutive steps involved in the computation are depicted, both for plane waves and focused regimes. Then, in Section 3,

numerical examples are described and the convergence of the method is investigated.

# 2. THEORETICAL APPROACH

In this section, a layered nanostructure, invariant along the *y* direction and periodic along the *x* direction (period *d*), is considered (see Fig. 1). The layer B, structured along the *x* direction, is surrounded by two semi-infinite media, A and C. The interface between A and B (resp. B and C) is noted *α* (resp. *β*). The incident plane wave is TE-polarized (i.e., the electric field is along the *y* direction) and impinging on the grating at oblique incidence, making an angle *θ* with the *z* axis. The SH wave is also TE-polarized. The theoretical method to compute the generated SH field in periodic nanostructures is detailed first for a single plane wave of frequency *ω*, and is then generalized to a focused beam.

## A. Method for Plane Wave SH Calculation

### 1. Construction of the Fundamental Field

In layer B, the field at the fundamental frequency is computed thanks to the BMM [19]. In this method, the electric field is interpolated on *N* samples on the *x* axis with a nonuniform basis of B-splines, which are piecewise continuous polynomials of degree *p* [22]. This field is determined through a modal equation $M_B e_m = k_m^2 e_m$ in TE-polarization, where $M_B$ is a $(N \times N)$ matrix deduced from the Maxwell and constitutive equations, the scalar $k_m$ is the *m*th mode wave vector along *z*, and the vector $e_m$ is the complex amplitude of the electric field of the *m*th mode at the sampling points. The eigenvectors are gathered in a matrix $E_B^m$, whose *m*th column corresponds to the eigenvector $e_m$. $E_B^m$ being computed, $H_B^m$ for the *x* component of the magnetic field can also be calculated. The complex *y* component $\mathcal{E}$ of the electric field and *x* component $\mathcal{H}$ of the magnetic field, at the sampling coordinates *x* and fixed
<figure><img src="image_2.png" /><figcaption>Fig. 1. Typical setup of a grating layer (B) surrounded by two semi-infinite media (A, C) separated by two interfaces α and β. It is illuminated by a TE-polarized plane wave $\mathcal{E}$. The incident wave vector $k$ makes an angle $\theta$ with the optical axis $z$ in the $xOz$ plane. The structure is infinite along the $y$ direction.</figcaption></figure>

height $z$, can be represented as superpositions of two waves, one upgoing and the other downgoing in the structure:

$$ \mathcal{E}(z) = \mathcal{E}^{+}(z) + \mathcal{E}^{-}(z), \quad (1) $$

$$ \mathcal{H}(z) = \mathcal{H}^{+}(z) - \mathcal{H}^{-}(z). \quad (2) $$

These terms are vectors of size $N$ because of the sampling along $x$ over one period. The fields are written at the time frequency $\omega$, thanks to their projections on the eigenmodes basis:

$$ \mathcal{E}^{\pm}(z) = E_{B}^{\omega} \operatorname{diag}(e^{\pm ik_{m}(z-z_{0})}) U^{\pm}(z_{0}), \quad (3) $$

$$ \mathcal{H}^{\pm}(z) = H_{B}^{\omega} \operatorname{diag}(e^{\pm ik_{m}(z-z_{0})}) U^{\pm}(z_{0}). \quad (4) $$

To find the total fields, it suffices to find the two vectors of the mode amplitudes $U^{\pm}$, which are given relatively to a height $z_0$, and to propagate them in the eigenmodes basis ($e^{-i\omega t}$ is omitted). To this end, the scattering matrices at each interface are determined thanks to the continuity of the tangential electric and magnetic fields (see Fig. 2) [23]. For instance, $S_{BA}^a$ stands for the transmission scattering matrix from medium A to layer B, and $S_{BB}^a$ for the reflection scattering matrix on interface $\alpha$ in layer B.

The electric field is thus obtained in layer B at any $z$:

$$ \mathcal{E}_B(z) = E_B^\omega [P_z^\alpha + P_z^\beta S_{BB}^\beta P_\beta^\alpha] \Gamma S_{BA}^a U^-(z_\alpha), \quad (5) $$

where $\Gamma = [I - S_{BB}^a P_z^\alpha S_{BB}^\beta P_\beta^\alpha]^{-1}$, and $P_z^\alpha$ (resp. $P_z^\beta$) is the propagator from interface $\alpha$ (resp. $\beta$) to height $z: P_{z_2}^{z_1} = \operatorname{diag}(e^{ik_z|z_2-z_1|})$. $U^-(z_\alpha)$ is the vector for the incident downgoing wave at interface $\alpha$ in medium A.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Scattering matrices for a three-layered structure.</figcaption></figure>

## 2. Construction of the SH Field

In the following, all the equations and fields are given at the doubled frequency $2\omega$. The grating layer B comprises a non-linear material for which it is assumed a second order nonlinear susceptibility of the type $\chi_{yy}^{(2)} = \chi^{(2)}$, and $\chi_{xyy}^{(2)} = \chi_{zyy}^{(2)} = 0$. Thus, it induces the complex second order nonlinear polarization [24]:

$$ \mathcal{P}^{(2)}(2\omega, z) = \epsilon_0 \chi^{(2)}(2\omega, z) \mathcal{E}_B^2(\omega, z). \quad (6) $$

This term generates a SH field $\mathcal{E}^{2\omega}$, which is calculated in the following. The total field at $2\omega$ is obtained by linear superposition of the fields generated in the structure by independent polarization layers. The undepleted pump approximation is used to state that the fundamental field stays unperturbed through the nonlinear process.

At every height $z$, a sub-layer of thickness $dz$ acts as a source creating a fraction of the SH field. Thus, to determine the field generated by a single sub-layer at height $z$, we first consider that it is the only source in the whole system. For the generated fraction of the $x$ component of the magnetic field $d\mathcal{H}$ at $2\omega$, the integration of the Maxwell-Ampere equation over regions of thickness $dz_0$ on a period $d$ gives

$$ d\mathcal{H}(z_+) - d\mathcal{H}(z_-) = -2i\omega \mathcal{P}^{(2)}(z) dz, \quad (7) $$

where $z_{\pm} = z \pm \frac{dz}{2}$. Then, Eq. (7) is projected on the eigen-mode basis thanks to Eqs. (2) and (4):

$$ \begin{aligned} & dU^{+}(z_{+}) - dU^{-}(z_{+}) - dU^{+}(z_{-}) + dU^{-}(z_{-}) \\ & = -2i\omega(\mathcal{H}_{B}^{2\omega})^{-1}\mathcal{P}^{(2)}(z)dz. \end{aligned} \quad (8) $$

In addition, the continuity of the electric field across the sub-layer, and the reflection scatterings at interfaces $\beta$ and $\alpha$ give the following equations (see Fig. 3):

<figure><img src="image_4.png" /><figcaption>Fig. 3. Scheme of the reflections and propagations undergone by the upgoing waves $dU^+$ (resp. the downgoing waves $dU^-$) in the nonlinear medium from height $z^+$ (resp. $z^-$). The sub-layer of thickness $dz$ at height $z$ generates a source term $d\Pi(z) = -i\omega \mathcal{P}^{(2)}(z) dz$.</figcaption></figure>
$$dU^{+}(z_{+}) + dU^{-}(z_{-}) = dU^{+}(z_{-}) + dU^{-}(z_{-}), \quad (9)$$

$$dU^{+}(z_{-}) = M^{\beta}dU^{-}(z_{-}), \quad (10)$$

$$dU^{-}(z_{+}) = M^{\alpha}dU^{+}(z_{+}), \quad (11)$$

where $M^i = P_z^i S_{BB}^i P_i^z$ are functions of $z$ ($i = \alpha, \beta$). Equations (8)–(11) can be reduced to the following system:

$$M^{\beta}dU^{-}(z_{-}) + d\Pi(z) = dU^{+}(z_{+}), \quad (12)$$

$$M^{\alpha}dU^{+}(z_{+}) + d\Pi(z) = dU^{-}(z_{-}), \quad (13)$$

where $d\Pi(z) = -i\omega\mathcal{P}^{(2)}(z)dz$. The contribution of a single sub-layer is then given by

$$dU^{+}(z_{+}) = [I - M^{\beta}M^{\alpha}]^{-1}[M^{\beta} + I]d\Pi(z), \quad (14)$$

$$dU^{-}(z_{-}) = [I - M^{\alpha}M^{\beta}]^{-1}[M^{\alpha} + I]d\Pi(z), \quad (15)$$

Working in the limit $dz \to 0$, one gets $z_+ = z_- = z$. Each sub-layer is considered linearly independent from the others. Then, the contributions of all the sub-layers are combined to obtain the value of the total SH field at $2\omega$:

$$\mathcal{E}^{+}(z) = E_{B}^{2\omega} \int_{z_{\beta}}^{z} P_{z}^{z'} dU^{+}(z'), \quad (16)$$

$$\mathcal{E}^{-}(z) = E_{B}^{2\omega} \int_{z}^{z_{\alpha}} P_{z}^{z'} dU^{-}(z'). \quad (17)$$

Nevertheless, the problem must be discretized for numerical processing, preventing the different sub-layers from being infinitely small. The infinitesimal quantity $dz$ must then be replaced by a finite one, $\Delta z$.

## B. Method for SH in Focused Regime

To design nanostructures, it is more significant to use focused beams as inputs rather than plane waves to take their experimental characterization into account. Here, a method is depicted to numerically model beams focused along the $x$ direction as they can be decomposed in a plane wave basis.

### 1. Construction of the Fundamental Field

The structure is supposed $d$-periodic along $x$, and the wavelength of the beam is fixed at $\lambda$. The incident beam is described by its electric field distribution $\mathcal{E}^i(x, z)$ in the $xOz$ plane as it is considered infinite along the $y$ direction (see Fig. 4). To gain

<figure><img src="image_5.png" /><figcaption>Fig. 4. Scheme of a cylindrical collimated beam. It is infinite along the $y$ axis.</figcaption></figure>

access to its plane wave spectrum, it can be written as an inverse Fourier transform in the reference plane $z = z_0$:

$$\mathcal{E}^i(x, z = z_0) = \frac{1}{2\pi} \int_{\mathbb{R}} \tilde{\mathcal{E}}^i(k_x) e^{ik_x x} d\mathbf{k}_x, \quad (18)$$

where the wave vector is $\mathbf{k} = (k_x, k_z)$ such that $|\mathbf{k}| = 2\pi/\lambda$. The Fourier coefficients are given by

$$\tilde{\mathcal{E}}^i(k_x) = \int_{\mathbb{R}} \mathcal{E}^i(x, z = z_0) e^{-ik_x x} dx. \quad (19)$$

An incident beam optically focused on a nanostructure does necessarily contain only propagative orders in its Fourier decomposition. It then reduces the space of integration to $k_x$ respecting the condition $|k_x| \le 2\pi/\lambda$, ensuring the existence of an associated real $k_z$:

$$\mathcal{E}^i(x, z = z_0) = \int_{|k_x| \le 2\pi/\lambda} \tilde{\mathcal{E}}^i(k_x) e^{ik_x x} d\mathbf{k}_x. \quad (20)$$

For the sake of numerical simulations, a discrete spectrum of $k_x$ is needed. To this end, a super-period of size $D = Md$ is considered. It contains a finite number $M$ of periods, such that $D$ is much larger than the spatial extent of the beam at $z = 0$. This super-period $D$ stands for the repetition period of the focused field: $\mathcal{E}(x+D, z) = \mathcal{E}(x, z)$. This way, the fields can be written as a Fourier series over a discrete spectrum. Equation (18) for the incident $D$-periodic field $\mathcal{E}_D^i$ takes the form

$$\mathcal{E}_D^i(x, z = z_0) = \sum_{|m|<D/\lambda} \tilde{\mathcal{E}}_m^i e^{ik_m x}, \quad (21)$$

where $k_m^x = 2\pi m/D$ and the wave vector is $\mathbf{k} = (k_x^m, k_z)$ such that $|\mathbf{k}| = 2\pi/\lambda$. Note that the condition for propagating modes now turns to $|m| < D/\lambda$. The Fourier coefficients are given by

$$\tilde{\mathcal{E}}_m^i = \frac{1}{D} \int_{-D/2}^{D/2} \mathcal{E}_D^i(x, z = z_0) e^{-ik_m^x x} dx. \quad (22)$$

The fundamental $m$ components of the field $\mathcal{E}_m^f(x, z)$, resulting from the interaction between the incoming plane waves $e^{ik_m^x x}$, and the nanostructures are then constructed through the modal method previously described. Each mode of the field respects the condition of pseudo-periodicity: $\mathcal{E}_m^f(x+d, z) = \mathcal{E}_m^f(x, z)e^{ik_m^x d}$. This gives the relative phase between the various mode wave vectors $k_m^x$ from one period of the super-period to another. In a given period $n$, the total field can then be written thanks to the linearity of Maxwell equations:

$$\mathcal{E}_n^\omega(x, z) = \sum_{|m|<D/\lambda} \tilde{\mathcal{E}}_m^i \mathcal{E}_m^f(x, z) e^{ik_m^x nd}. \quad (23)$$

This description working for cylindrical propagating beams, whose field distribution respects the $1D$ wave equation, can be generalized to $2D$ beams [25].

### 2. Construction of the SH Field

The following computation is similar to the previous construction of the fundamental field with some substantial
differences: the first being the wavelength which is now equal to $\lambda/2$. The first step consists of writing the complex second order polarization as a Fourier series:

$$P_D^{(2)}(x, z) = \sum_{m'=-\infty}^{m'=+\infty} \tilde{P}_{m'}^{(2)}(z) e^{ik_{x}^{m'} x}, \quad (24)$$

where $k_x^{m'} = 2\pi m'/D$ and the wave vector is $\mathbf{k} = (k_x^{m'}, k_z)$ such that $|\mathbf{k}| = 4\pi/\lambda$. All the modes, including the nonpropagative ones, are taken into account. The Fourier coefficients are given by

$$\tilde{P}_{m'}^{(2)}(z) = \frac{1}{D} \int_{-D/2}^{D/2} P_D^{(2)}(x, z) e^{-ik_x^{m'} x} dx. \quad (25)$$

In Section 2.A.2, the nonlinear medium was decomposed into sub-layers generating fractions of the SH field from the source term $P^{(2)}(x, z)$. Here each sub-layer generates a fraction of the SH field which is the sum of SH fields created from source terms for mode $m'$: $P_{m'}^{(2)}(z)e^{-ik_{x}^{m'} x}$. The term $d\Pi(z)$ now depends on $m'$:

$$d\Pi_{m'}(z) = -i\omega \tilde{P}_{m'}^{(2)}(z) e^{-ik_{x}^{m'} x} dz. \quad (26)$$

The SH field is then computed in one period for every mode $m'$, thanks to the method described in Section 2.A.2. The latter is found in the other periods through the pseudo-periodicity condition written for this SH field: $\mathcal{E}_{m'}^{\text{shg}}(x+d, z) = \mathcal{E}_{m'}^{\text{shg}}(x, z)e^{ik_{x}^{m'} d}$. The total field in the nth period eventually reads

$$\mathcal{E}_n^{2\omega}(x, z) = \sum_{m'=-\infty}^{m'=+\infty} \mathcal{E}_{m'}^{\text{shg}}(x, z) e^{ik_x^{m'} nd}. \quad (27)$$

This result does not depend on the form of the nonlinear polarization, but is nonetheless limited to the undepleted pump regime. Propagation throughout the structure is also supposed to remain unperturbed by the presence of a nonlinear polarization term.

# 3. NUMERICAL RESULTS

## A. Maker Fringes

Numerical simulations are performed in the case of Maker fringes in a 110 crystalline GaAs membrane (see Fig. 5), where a square sine function at a spatial frequency $\Delta k = k_{2\omega} - 2k_{\omega} = 2(n_{2\omega} - n_{\omega})\omega/c$ is expected in the direction $z$ of the incident plane wave [24,26]. The incident electric field is polarized along a 111 axis, and the following parameters are used: $\lambda = 2\pi c/\omega = 5$ µm, $n_{\omega}^2 = \epsilon_{\text{GaAs}}(\omega) = 11.01$, $n_{2\omega}^2 = \epsilon_{\text{GaAs}}(2\omega) = 11.19$ [27] and $\chi^{(2)} = 2/\sqrt{3} \times 370$ pm/V (for polarization along the 111 axis) [24]. The thickness of the crystal is equal to three times the coherence length $L_{\text{coh}} = \pi/\Delta k$, and both faces have an anti-reflective coating at $\omega$ and $2\omega$, giving

$$|\mathcal{E}^{2\omega}(z)|^2 \propto \left| e^{-ik_{2\omega}z} \frac{\sin\left(\frac{\Delta k z}{2}\right)}{\Delta k} - \frac{\sin\left(\frac{\Delta k z}{2}\right)}{\Sigma k} \right|^2, \quad (28)$$

where $\Sigma k = k_{2\omega} + 2k_{\omega}$ appears to be the frequency of the negligible oscillations which can be seen on the inset. They

<figure><img src="image_6.png" /><figcaption>Fig. 5. Absolute square value of the SH electric field, normalized by the maximum value of the analytical curve for Maker fringes in 110 crystalline GaAs. The inset zooms in the oscillations at frequency $\Sigma k$. Involved parameters are $2\pi c/\omega = 5$ µm, $\epsilon_{\text{GaAs}}(\omega) = 11.01$, and $\epsilon_{\text{GaAs}}(2\omega) = 11.19$, $\theta = 0^{\circ}$.</figcaption></figure>

correspond to the sum of the contributions of upgoing waves ($dU^+$).

A Gaussian focused beam at $\omega$ is now incident onto the same structure, with a waist $W = 5$ µm and carrying a power of $P = 10$ mW. The expected solution for an incoming 2D collimated Gaussian beam is a Bessel function multiplied by a Gaussian envelope [24]. Here the result must be adapted to take 1D collimation into account, giving an approximate normalized SH field distribution:

$$\mathcal{E}_{\text{foc}}^{2\omega}(x, z) \propto \frac{A(x, z)}{\sqrt{q_0 + z}} e^{-ik_{\omega}x^2/(q_0 + z) - i2k_{\omega}z}, \quad (29)$$

where $q_0 = i\pi W^2/\lambda$, and

$$A(x, z) = \frac{i\omega}{\sqrt{\epsilon_2 c}} \chi^{(2)} \int_0^z \frac{e^{i\Delta k z'}}{\sqrt{q_0 + z'}} dz'. \quad (30)$$

Figure 6(c) shows the evolution of the square of the SH field absolute value at the center of the beam with respect to $z$. The waist of the beam, at the entrance of the medium, is equal to $W = 5$ µm. The computation is performed over a depth equal to three times the coherence length of the SHG process. The analytical and numerical curves are found to be in good agreement. The difference is mainly because of the fact that the analytical solution relies on the hypothesis that the SH field has the same confocal parameter $b = k_{\omega}W^2$ as the fundamental beam. The absolute values of the fundamental and SH fields are also displayed in Figs. 6(a) and 6(b) as functions of $x$ and $z$ to show the shape of the beam. These amplitudes are normalized to the maximum value of the field at $\omega$.

## B. SHG Exaltation in a Grating of Nanorods

As this method aims at designing nanostructures able to greatly enhance the nonlinear processes, a single grating of cylindrical nanorods is studied. The behavior of this architecture has already been investigated both in linear [28,29] and nonlinear cases [10]. These works highlighted a need for a strong linear exaltation of the electric field inside the cylinders. For a given wavelength $\lambda$ and diameter $R$, such a nanostructure exhibits reflectivity resonances for peculiar values
<figure><img src="image_7.png" /><figcaption>Fig. 6. Fields maps (a) of the focused fundamental field and (b) of the SH field in 110 crystalline GaAs, normalized to the maximum amplitude of the fundamental field. (c) Absolute square value of the SH electric field at $x = 150$ µm as a function of $z$, normalized to the maximum value of the analytical curve. The waist $W$ of the Gaussian beam is located at the entrance of the crystal and is equal to 5 µm. Involved parameters are $2\pi c/\omega = 5$ µm, $\epsilon_{GaAs}(\omega) = 11.01$, $\epsilon_{GaAs}(2\omega) = 11.19$, $\chi^{(2)} = 2/\sqrt{3} \times 370$ pm/V, and $\theta = 0^\circ$. The beam carries a power $P = 10$ mW.</figcaption></figure>

of the period $d$. At $2\pi c/\omega = 5$ µm and $R = 0.2d$, a resonance is found for a sub-wavelength period, $d = 4.36$ µm. Other involved parameters are $\epsilon_{rods} = 2$, $\theta = 8.25^\circ$, $\chi^{(2)} = 2$ pm/V incident field, and $|\mathcal{E}^\omega| = 1 \times 10^6$ V/m. For the parameters $\Delta z = 0.01R$ and $N = 200$ (corresponding to the situation $\Delta z \sim 0.002\lambda$), the resulting field maps to be compared with Ref. [11] are displayed in Figs. 7(a) and 7(b). Circles are fitted with crosses as drawn on the field maps: the cross-shaped section of the rods has total sides of length $R$ and branches of length $R/4$. The difference between these results and those of the Ref. [11] can be explained by the fact that circular sections are identified with cross-shaped sections.

<figure><img src="image_8.png" /><figcaption>Fig. 7. Fundamental and SH absolute values of the electric field in the plane perpendicular to the rod axis of a grating of GaAs nanorods. The electric field is polarized along the rods. Involved parameters are period $d = 4.36$ µm, $R = 0.2d$, $2\pi c/\omega = 5$ µm, $\chi^{(2)} = 2$ pm/V, $\epsilon_{rods} = 2$, $\theta = 8.25^\circ$, and incident field $|\mathcal{E}^\omega| = 1 \times 10^6$ V/m. The cross-shaped section of the rods has total sides of length $R$ and branches of length $R/4$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Relative error on the value of the output fundamental intensity as a function of the number of interpolation points $N$, with the reference value $I^\omega(N = 1000) = 17.33$ GW/m$^2$. The structure is a grating of GaAs nanorods; involved parameters are period $d = 4.90$ µm, $R = 0.2d$, $2\pi c/\omega = 5$ µm, $\epsilon_{GaAs}(\omega) = 11.01$, $\epsilon_{GaAs}(2\omega) = 11.19$, $\theta = 1^\circ$, and incident field $|\mathcal{E}^\omega| = 1 \times 10^6$ V/m.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Relative error on the value of the output zeroth order SH intensity in the zero diffractive order (a) as a function of the number $N$ of samples along $x$ in a period for $\Delta z = 10^{-4}R$, and (b) as a function of $\Delta z$, the discretization pace along $z$ direction for $N = 1000$. The reference value is $I^2(\Delta z = 10^{-4}R, N = 1000) = 1.817757$ W/m$^2$. The structure is a grating of GaAs nanorods; involved parameters are the same as those in Fig. 8 with $\chi^{(2)} = 150$ pm/V.</figcaption></figure>
### C. Convergence for a Grating of GaAs Nanorods

The convergence of the method is now investigated in the case of a grating of nanorods as a function of the number *N* of samples in a period along *x* and as a function of the discretization pace along *z* direction $\Delta z$. A grating (period *d* = 4.9 µm) of GaAs nanorods with square section (*R* = 0.2d) is studied. These nanorods are fabricated out of a 110 GaAs substrate and they are oriented along a 111 direction. The incoming light ($2\pi c/\omega$ = 5 µm) is TE-polarized and quasi-normally incident ($\theta$ = 1°), and the amplitude of the incident electric field is $|\mathcal{E}^\omega| = 1 \times 10^6$ V/m. First, to assess the convergence of the linear computed field, the relative error on the fundamental transmitted intensity is plotted in Fig. 8 as a function of *N*. Second, the relative error made on the value of the output zeroth order SH intensity as a function of *N* is displayed in Fig. 9 for $\Delta z = 10^{-4}R$, where $I^{2\omega}(\Delta z = 10^{-4}R, N = 1000)$ is taken as the reference value. This error decreases as *N* increases and goes under $10^{-4}$ from $N \approx 200$. This relative error is also displayed as a function of $\Delta z$ for a fixed $N = 1000$, with the same reference value. The same behavior is observed, highlighting the fact that the choices of $\Delta z$ and *N* have to be balanced to obtain an accurate result.

## 4. CONCLUSION

A method for computing nonlinear effects in multilayered nanostructured media has been presented in this paper. On one hand, the method for plane waves works under the undepleted pump approximation and is efficient for a properly chosen *z*-sampling. It has been shown to converge with an increasing number of interpolations points and decreasing discretization pace along *z*. On the other hand, numerical simulations under a focused beam have been achieved by decomposing the incident light onto a plane wave basis. These results can be extended to 3D focused beams. Recent developments in the B-spline modal method have also introduced a Galerkin description which was shown to give even better convergence rates for linear problems [21].

## REFERENCES

1. P. Bouchon, F. Pardo, B. Portier, L. Ferlazzo, P. Ghenuche, G. Dagher, C. Dupuis, N. Bardou, R. Haïdar, and J. Pelouard, "Total funneling of light in high aspect ratio plasmonic nanoresonators," Appl. Phys. Lett. 98, 191109 (2011).

2. T. Grosjean, M. Mivelle, F. Baida, G. Burr, and U. Fischer, "Diabolo nanoantenna for enhancing and confining the magnetic optical field," Nano Lett. 11, 1009-1013 (2011).

3. L. Novotny and N. van Hulst, "Antennas for light," Nat. Photonics 5, 83-90 (2011).

4. P. Chevalier, P. Bouchon, J.-J. Greffet, J.-L. Pelouard, R. Haïdar, and F. Pardo, "Giant field enhancement in electromagnetic Helmholtz nanoantenna," Phys. Rev. B 90, 195412 (2014).

5. S. Kim, J. Jin, Y.-J. Kim, I.-Y. Park, Y. Kim, and S.-W. Kim, "High-harmonic generation by resonant plasmon field enhancement," Nature 453, 757-760 (2008).

6. J. A. Schuller, E. S. Barnard, W. Cai, Y. C. Jun, J. S. White, and M. L. Brongersma, "Plasmonics for extreme light concentration and manipulation," Nat. Mater. 9, 193-204 (2010).

7. P. Genevet, J.-P. Tetienne, E. Gatzogiannis, R. Blanchard, M. A. Kats, M. O. Scully, and F. Capasso, "Large enhancement of nonlinear optical phenomena by plasmonic nanocavity gratings," Nano Lett. 10, 4880-4883 (2010).

8. M. Kauranen and A. V. Zayats, "Nonlinear plasmonics," Nat. Photonics 6, 737-748 (2012).

9. J. Bravo-Abad, S. Fan, S. G. Johnson, J. D. Joannopoulos, and M. Soljačić, "Modeling nonlinear optical phenomena in nanophotonics," IEEE J. Lightwave Technol. 25, 2539-2546 (2007).

10. A. Borisov, D. Marinica, and S. Shabanov, "Second harmonic generation from arrays of sub-wavelength cylinders," Phys. Rev. B 76, 085311 (2007).

11. L. Yuan and Y. Y. Lu, "Analyzing second harmonic generation from arrays of cylinders using dirichlet-to-neumann maps," J. Opt. Soc. Am. B 26, 587-594 (2009).

12. J. Mäkitalo, S. Suuriniemi, and M. Kauranen, "Boundary element method for surface nonlinear optics of nanoparticles," Opt. Express 19, 23386-23399 (2011).

13. J. Butet, B. Gallinet, K. Thyagarajan, and O. J. Martin, "Second-harmonic generation from periodic arrays of arbitrary shape plasmonic nanostructures: a surface integral approach," J. Opt. Soc. Am. B 30, 2970-2979 (2013).

14. J. Francés, J. Tervo, and A. Márquez, "Tensorial split-field finite-difference time-domain approach for second-and third-order nonlinear materials," J. Opt. Soc. Am. B 30, 1711-1719 (2013).

15. W. Nakagawa, R.-C. Tyan, and Y. Fainman, "Analysis of enhanced second-harmonic generation in periodic nanostructures using modified rigorous coupled-wave analysis in the undepleted-pump approximation," J. Opt. Soc. Am. A 19, 1919-1928 (2002).

16. 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).

17. 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).

18. K. Edee, P. Schiavone, and G. Granet, "Analysis of defect in extreme uv lithography mask using a modal method based on nodal b-spline expansion," Jpn. J. Appl. Phys. 44, 6458-6462 (2005).

19. P. Bouchon, F. Pardo, R. Haïdar, and J.-L. Pelouard, "Fast modal method for subwavelength gratings based on b-spline formulation," J. Opt. Soc. Am. A 27, 696-702 (2010).

20. M. Walz, T. Zebrowski, J. Küchenmeister, and K. Busch, "B-spline modal method: a polynomial approach compared to the Fourier modal method," Opt. Express 21, 14683-14697 (2013).

21. G. Granet, "Efficient implementation of b-spline modal method for lamellar gratings," J. Opt. Soc. Am. A 31, 332-337 (2014).

22. C. De Boor, *A Practical Guide to Splines* (Springer-Verlag, 1978), Vol. 27.

23. L. Li, "Multilayer modal method for diffraction gratings of arbitrary profile, depth, and permittivity," J. Opt. Soc. Am. A 10, 2581-2591 (1993).

24. R. W. Boyd, *Nonlinear Optics* (Academic, 2002).

25. P. Chevalier, P. Bouchon, F. Pardo, and R. Haïdar, "Electromagnetic modeling of spherical focusing on a one-dimensional grating thanks to a conical b-spline modal method," J. Opt. Soc. Am. A 31, 1692-1697 (2014).

26. P. Maker, R. Terhune, M. Nisenoff, and C. Savage, "Effects of dispersion and focusing on the production of optical harmonics," Phys. Rev. Lett. 8, 21-22 (1962).

27. A. Kachare, W. Spitzer, J. Fredrickson, and F. Euler, "Measurements of layer thicknesses and refractive indices in high-energy ion-implanted GaAs and GaP," J. Appl. Phys. 47, 5374-5381 (1976).

28. R. Gómez-Medina, M. Laroche, and J. J. Sáenz, "Extraordinary optical reflection from sub-wavelength cylinder arrays," Opt. Express 14, 3730-3737 (2006).

29. P. Ghenuche, G. Vincent, M. Laroche, N. Bardou, R. Haïdar, J.-L. Pelouard, and S. Collin, "Optical extinction in a single layer of nanorods," Phys. Rev. Lett. 109, 143903 (2012).