Revisiting the boundary conditions for second-harmonic generation at metal-dielectric interfaces

TL;DR

This paper revisits the boundary conditions for second-harmonic generation at metal-dielectric interfaces, clarifying the discontinuity of electric fields, and provides a correction method for numerical simulations.

Contribution

It introduces a spectral decomposition method that simplifies boundary condition considerations and offers a correction recipe for commercial simulation tools.

Findings

Confirmed the discontinuity of the tangential electric field component at the interface.

Reaffirmed the standard convention for surface dipole radiation.

Provided a correction method for COMSOL Multiphysics to accurately model SHG.

Abstract

We study second-harmonic generation (SHG) arising from surface nonlinearity at a metal-dielectric interface using a spectral decomposition method. Since our method avoids the need to consider the generalized boundary condition across the metal-dielectric interface in the presence of a perpendicular surface source, we retrieve the known discontinuity of the tangential component of the electric field ($E_{\parallel}^ {2\omega}$) for a general geometry, based on a purely mathematical argument. Further, we reaffirm the standard convention of the implementation of this condition, namely, that the surface dipole source radiates as if placed outside the metal surface for arbitrary geometries. We also study and explain the spectral dependence of the discontinuity of the tangential component of the electric field at second harmonic. Finally, we note that the default settings of the commercial numerical package COMSOL Multiphysics fail to account for the $E_{\parallel}^ {2\omega}$-discontinuity. We provide a simple recipe that corrects the boundary condition within these existing settings.

Equations

$${\mathbf{P}}^{(2)}(\mathbf{r})={\mathbf{P}}_S^{(2)}({\mathbf{r}}_{\parallel })\delta (r_{\perp }),$$

$$\Delta D_{\perp }$$

$$\Delta {\mathbf{E}}_{\parallel }$$

$$\Delta B_{\perp }$$

$$\Delta {\mathbf{H}}_{\parallel }$$

$${\mathbf{P}}_S^{(2)}={\epsilon }_0{\chi }_{S,\perp ,\perp ,\perp }^{(2)}{E}_{\perp }^{\omega }({\mathbf{r}}_{\parallel })E_{\perp }^{\omega }({\mathbf{r}}_{\parallel }){\hat{r}}_{\perp },$$

$$\nabla \times (\nabla \times \mathbf{E})-k_0^2{\epsilon }_{bg}\mathbf{E}=k_0^2\left[\epsilon (\mathbf{r})-{\epsilon }_{bg}\right]\mathbf{E}+i\omega {\mu }_0{\mathbf{J}}_f,$$

$$\mathbf{E}(\mathbf{r})={\mathbf{E}}_0(\mathbf{r})+k_0^2\left[{\epsilon }_{in}-{\epsilon }_{bg}\right]\int {\overline{\stackrel{ˉ}{\mathcal{G}}}}_0(|\mathbf{r}-{\mathbf{r}}^{\prime }|)\Theta ({\mathbf{r}}^{\prime })\mathbf{E}({\mathbf{r}}^{\prime })\text{leavevmode}d{\mathbf{r}}^{\prime },$$

$${\mathbf{E}}_0(\text{r})=i\omega {\mu }_0\int {\overline{\stackrel{ˉ}{\mathcal{G}}}}_0(|\mathbf{r}-{\mathbf{r}}^{\prime }|)\text{leavevmode}{\mathbf{J}}_f({\text{r}}^{\prime })d{\mathbf{r}}^{\prime }.$$

$${\overline{\stackrel{ˉ}{\mathcal{G}}}}_0(|\mathbf{r}-{\mathbf{r}}^{\prime }|)=\left[\overline{\stackrel{ˉ}{\mathbf{I}}}-\frac{1}{k_0^2{\epsilon }_{bg}}\nabla \nabla \cdot \right]{\mathcal{G}}_0.$$

$$s_l{\mathbf{E}}_l(\mathbf{r})=k_0^2\int {\overline{\stackrel{ˉ}{\mathcal{G}}}}_0(|\mathbf{r}-{\mathbf{r}}^{\prime }|)\Theta ({\mathbf{r}}^{\prime }){\mathbf{E}}_l({\mathbf{r}}^{\prime })d{\mathbf{r}}^{\prime },$$

$$s_l\equiv \frac{1}{{\epsilon }_l-{\epsilon }_{bg}}.$$

$$|\mathbf{E}⟩=|{\mathbf{E}}_0⟩+\sum _{l\in \mathbb{N}}|{\mathbf{E}}_l⟩\left(\frac{{\epsilon }_{in}-{\epsilon }_{bg}}{{\epsilon }_l-{\epsilon }_{in}}\right)⟨{\mathbf{E}}_l|\hat{\theta }|{\mathbf{E}}_0⟩,$$

$${\mathbf{J}}_m^{(2)}=-\frac{i\omega {\epsilon }_0}{\pi }{\chi }_{S,\perp \perp \perp }^{(2)}\left[{\int }_0^{2\pi }{\left[E_{\rho }^{\omega }\left(a,{\theta }^{\prime }\right)\right]}^2{e}^{-im{\theta }^{\prime }}d{\theta }^{\prime }\right]$$

$$\times e^{im\theta }\text{leavevmode}\delta (\rho -a)\text{leavevmode}\hat{\rho }.$$

$$|{\mathbf{E}}_0^{2\omega }(\rho <a,\theta )⟩$$

$$|{\mathbf{E}}_0^{2\omega }(\rho >a,\theta )⟩$$

$$A_m=\frac{ia}{8{\epsilon }_{bg}^{2\omega }}{\chi }_{S,\perp ,\perp ,\perp }^{(2)}\left[{\int }_0^{2\pi }{\left[E_{\rho }^{\omega }\left(a,{\theta }^{\prime }\right)\right]}^2{e}^{-im{\theta }^{\prime }}d{\theta }^{\prime }\right].$$

$$\Delta E_{0,\theta }^{2\omega }=\Delta E_{\theta }^{2\omega }=\sum _{m\in \mathbb{Z}}\frac{-im}{2\pi a{\epsilon }_{bg}^{2\omega }}{\chi }_{S,\perp ,\perp ,\perp }^{(2)}{e}^{im\theta }$$

$$\times \left[{\int }_0^{2\pi }{\left[E_{\rho }^{\omega }\left(a,{\theta }^{\prime }\right)\right]}^2{e}^{-im{\theta }^{\prime }}d{\theta }^{\prime }\right].$$

$${\kappa }_m$$

$$\kappa =\text{max}\left|\frac{E_{\theta }^{2\omega }(\rho =a+\delta a,\theta )-E_{\theta }^{2\omega }(\rho =a-\delta a,\theta )}{\text{max}\left[E_{\theta }^{2\omega }(\rho =a+\delta a,\theta )+E_{\theta }^{2\omega }(\rho =a-\delta a,\theta )\right]}\right|,$$

$${\rho }^{(2)}=-\nabla \cdot {\mathbf{P}}^{(2)},{\sigma }^{(2)}={\mathbf{P}}^{(2)}\cdot r_{\perp },$$

$$\int {\rho }^{(2)}dV$$

$$=$$

$$=\int {\sigma }^{(2)}dS$$

$${\varrho }^{(2)}=-\left[{\nabla }_{\parallel }\cdot {\mathbf{P}}_S^{(2)}\right]\delta (r_{\perp }),{\varsigma }^{(2)}=-{\nabla }_{\parallel }\cdot {\mathbf{P}}_S^{(2)}.$$

$$\nabla \times (\nabla \times {\mathbf{E}}^{2\omega })-k^2(\mathbf{r}){\mathbf{E}}^{2\omega }=2i{\mu }_0\omega J_{S,\perp }^{(2)},k^2(\mathbf{r})=4k_0^2{\epsilon }^{2\omega }(\mathbf{r}).$$

$${\mathbf{E}}^{2\omega }={\mathbf{E}}_{bg}^{2\omega }+{\mathbf{E}}_{sc}^{2\omega },$$

$$\nabla \times (\nabla \times {\mathbf{E}}_{bg}^{2\omega })-(k_{bg}^{2\omega }{)}^2{\mathbf{E}}_{bg}^{2\omega }=2i\omega {\mu }_0{J}_{S,\perp }^{(2)}.$$

$$\nabla \times (\nabla \times {\mathbf{E}}_{sc}^{2\omega })-k^2(\mathbf{r}){\text{E}}_{sc}^{2\omega }=[k^2(\mathbf{r})-(k_{bg}^{2\omega }{)}^2]{\mathbf{E}}_0^{2\omega }.$$

Full text

††journal: josab

1 Introduction

Second-order optical nonlinear processes in centrosymmetric materials have been studied since the early days of nonlinear optics [1, 2, 3, 4, 5, 6, 7]. It is well known that these processes are symmetry forbidden in the local bulk response limit, and thus are governed by higher-order nonlocal bulk and surface symmetry breaking effects [8, 9, 10]. Recently, there has been renewed interest in surface nonlinearities from metal-dielectric interfaces, in which electric fields can be greatly enhanced due to plasmon resonances. Specifically, advances in fabrication of nanostructures have allowed the experimental investigations of second-harmonic generation (SHG) from single metallic nanoparticles [11, 12, 13, 14, 15, 16], split-ring resonators [17, 18], periodic nanostructured metal films  [19], and other nano-geometries [see Ref. 20, and references therein for a more thorough review]. The corresponding theoretical modeling relied either on surface nonlinearity [21, 22, 23] or non-local bulk nonlinearity [24, 25], while some studies have incorporated both effects [26, 27, 28, 29, 30, 31, 32, 33, 34, 35]. Importantly, even the bulk nonlinearity can be mapped to a surface current source, so that the surface nonlinearity provides a simple, general model for second-order nonlinear phenomena in metals [36, 35, 33, 34].

In Ref. [37], Heinz showed that a non-uniform, normal component of the surface current results in a discontinuous tangential electric field ($\mathbf{E}_{\parallel}$) across the interface between the two media. Computing the nonlinear polarization at the second-harmonic (SH) can involve the normal component of the electric field at the fundamental frequency (FF). Since the normal component of the electric field is discontinuous across the interface, it was a-priori not clear whether one has to choose the fields inside or outside the metal for evaluating the SH electromagnetic fields accurately. This was resolved by Sipe et al. as follows: one has to evaluate the nonlinear polarization based on the fundamental fields inside the metal; the dipole layer source generated from the nonlinear polarization should then be placed just outside the metal surface [10].

Here, we revisit the derivation of the generalized boundary condition (BC) in the presence of surface sources and provide an alternative description of surface SHG using the spectral decomposition method (SDM) in a full electrodynamic setting [38]. We show that the generalized BC is naturally reproduced in the SDM formulation without invoking Heinz’s generalized BC [37]. Based on SDM, we show from a purely mathematical point of view and for a general geometry that the dipole layer radiates indeed as if placed outside the metal surface, in agreement with Sipe et al. [10]. As a specific example, we apply SDM to an infinite metallic cylinder under normal incidence and provide physical insight into the spectral dependence of the discontinuity of the parallel component of the SH electric field, $\textbf{E}_{\parallel}^{2\omega}$. In particular, we show that the $\textbf{E}_{\parallel}^{2\omega}$-discontinuity is maximal away from the SH resonance, yielding 50%–100% errors if not implemented correctly; yet, the discontinuity is smaller at a SH resonance, its value diminishing as the quality factor of the resonance improves. Finally, we note that the default settings of COMSOL Multiphysics, a commonly used commercial solver of Maxwell’s Equations, neglect the normal component of the surface current $J_{S,\perp}^{(2)}$ and fail to account for the $\mathbf{E}_{\parallel}^{2\omega}$-discontinuity. As a remedy, we provide a simple recipe that accounts for the error within the existing settings using the insight provided by SDM.

The paper is organized as follows: In Sec. 2 we review the generalized BCs for electromagnetic waves and their implementation. Sec. 3 briefly introduces the spectral decomposition method (SDM) and its application to SHG. We calculate the SHG from a metallic infinite cylinder with a surface nonlinearity using SDM in Sec. 4 and study the spectral dependence of the $\textbf{E}_{\parallel}^{2\omega}$-discontinuity field in Sec. 5. An alternative approach for the numerical implementation (in COMSOL Multiphysics) of the generalized BC is discussed in Sec. 6. We conclude with a discussion and outlook.

2 GENERALIZED BOUNDARY CONDITIONS AND THEIR IMPLEMENTATION

Consider a surface SH polarization vector at a metal-dielectric interface given by

[TABLE]

where $\delta({{r}}_{\perp})$ is the Dirac-delta function [37]. Note that we denote the vector and scalar quantities in bold and italic font, respectively. The subscript $\parallel$ ($\perp$) corresponds to the direction tangential (normal) to the interface. For example, for a flat interface at $z=0$, one would have $\mathbf{r}_{\parallel}=(x,y)$ and $r_{\perp}=z=0$. The SH polarization (1) can also be formulated as a surface current density, $\mathbf{J}_{S}^{(2)}=-2i\omega\mathbf{P}_{S}^{(2)}$, or as surface charge $\varsigma^{(2)}$ or volume charge $\varrho^{(2)}$ densities, see Appendix A.

In the standard formulation of the Maxwell´s Equations for SH fields, $\mathbf{P}^{(2)}_{S}$ appears only in the BCs at the metal-dielectric interface as an external source. The tangential component of $\mathbf{P}_{S}^{(2)}$, denoted by $\mathbf{P}_{S,\parallel}^{(2)}$, gives rise to a surface current $\mathbf{J}_{S,\parallel}^{(2)}$ tangential to the interface. Surface currents $\mathbf{J}_{S,\parallel}^{(2)}$ are very common in electrodynamic problems and the way to incorporate them in the BCs across the interface is described in any standard textbook, see e.g., [39]. On the other hand, the normal component ${P}_{S,\perp}^{(2)}$ corresponds to a perpendicular surface zero-thickness current, ${J}_{S,\perp}^{(2)}$, or equivalently, to a dipole layer pointing normal to the interface. Such a component is unusual and requires a more careful treatment. In the electrostatic limit, Jackson [39] derived the scalar electric potential for a uniform surface dipole layer oriented along the normal direction in free space [40]. Heinz [37] generalized this BC to an electrodynamic setting (i.e., for Maxwell´s Equations rather than Laplace´s Equation). Specifically, he modeled the surface sources in the zero thickness limit of a thin “pill box” in which the dipole layer is placed. Heinz’s theory generalizes Jackson’s BCs [39] in the following aspects:

The permittivity of the pill box $\varepsilon^{\prime}$ is allowed to be different from vacuum. 2. 2.

The surface dipole layer is allowed to be inhomogeneous (along the transverse coordinates).

The BCs across the interface then read

[TABLE]

where $\Delta$ denotes the difference of the variable across the interface. $\nabla_{\parallel}$ and $\hat{r}_{\perp}$ correspond to the tangential component of the vector derivative operator and the unit vector normal to the interface, respectively. We see that Eq. (2), and Eqs. (4)-(5) are just the usual BCs found in standard textbooks (see e.g., [39]). However, Eq. (3) is more general. In particular, it shows that the spatial variation of the nonlinear polarization along the transverse coordinates results in discontinuity of $\mathbf{E}_{\parallel}$. Importantly, Eq. (3) reduces to the standard BC when the normal component of $\mathbf{P}^{(2)}_{S}$ has zero gradient along the tangential direction, i.e., for $\nabla_{\parallel}{P}^{(2)}_{S,\perp}=0$. Moreover, note that the normal component of the electric field $E_{\perp}$ in Eq. (2) becomes singular at the due to the presence of the dipole layer. Finally, Equation (3) was also derived based on Green’s function analysis in the context of planar structures by Sipe [41].

We now discuss the implementation of the BC given by Eq. (3). To capture the unusual aspects of this general BC for SH radiation in a simple configuration, we choose to work with a $\mathbf{P}^{(2)}_{S}$ given by

[TABLE]

where $E_{\perp}^{\omega}({\mathbf{r}}_{\parallel})$ is the normal component of the electric field at the FF. Such a polarization corresponds to a purely surface source, which was shown to be the dominant contribution to the $\bar{\bar{\chi}}^{(2)}_{S}$ tensor [42, 16].

Note that $\chi^{(2)}_{S,\perp,\perp,\perp}$ is in general frequency dependent [9, 16]. However, the SH electric fields in all our results are normalized to their maximum values, thus, the actual value of $\chi^{(2)}_{S,\perp,\perp,\perp}$ is not important. Recall that the RHS of Eq. (3) is non-zero only when $\mathbf{P}^{(2)}_{S}$ has non-zero gradient along the tangential direction, thus yielding a discontinuous parallel SH field, $\mathbf{E}^{2\omega}_{\parallel}$. This is the case in Eq. (6) when applied to any curved interface and a flat interface under oblique FF incidence.

The generalized BC in Eq. (3) leads to an ambiguity since $\mathbf{P}^{(2)}_{S}$ in Eq. (6) depends on ${E}_{\perp}^{\omega}$ at the boundaries, which is inherently discontinuous due to the induced charges at the FF. Therefore, it is a-priori not clear if ${E}_{\perp}^{\omega}$ should be evaluated on the metal or on the dielectric side of the interface. In order to resolve this seeming arbitrariness, Sipe et al. argued that the FF fields in the nonlinear polarization are to be chosen inside the metal [10, 43]. It was also shown that the resulting dipole originating from $\mathbf{P}^{(2)}_{S}$ radiates as if placed outside the metal surface. The placement of this dipole layer outside the metal corresponds to $\varepsilon^{\prime}=\varepsilon_{bg}^{2\omega}$ in the generalized BC in Eq. (3). Mizrahi et al. established a proper correspondence between Heinz’s and Sipe’s derivation for planar slab geometries [44].

One of the main goals of this paper is to provide a simple method to calculate SHG without invoking Eq. (3) and also to confirm its correspondence with the emergence of a SH dipole layer source outside the metal from a purely mathematical point of view for arbitrary geometries. We use a spectral decomposition method to achieve this.

3 Spectral Decomposition method (SDM)

We now outline our approach to calculate the SHG based on the spectral decomposition method [38]. In short, SDM is an eigenmode expansion solution for the electrodynamic fields produced by an arbitrary current source even in the presence of a lossy resonator in an open system. In what follows, we briefly discuss the full electrodynamic formulation of SDM for non-magnetic and isotropic media. We then outline the full electrodynamic SDM solution for SHG arising from $\mathbf{P}^{(2)}_{S}$ [Eq. (6)].

The SDM [38] is based on the electrodynamic analogue of the Lippmann-Schwinger equation. It begins with the standard inhomogeneous Maxwell’s equations for a monochromatic wave, i.e., the vector Helmholtz equation, but rearranged to treat the response of the structure as another source term, namely,

[TABLE]

where $\varepsilon(\mathbf{r})$ is the spatial distribution of the permittivity and $k_{0}=\omega/c$. Here, we focus our attention on structures that have only two constituent materials, loosely defined as the inclusion ($\varepsilon_{in}$) and background ($\varepsilon_{bg}$). Thus, the free current source ($\mathbf{J}_{f}$) and the bound (displacement) current response (proportional to [$\varepsilon_{in}-\varepsilon_{bg}$]) terms are treated on an equal footing. The Green’s function solution to Eq. (7) is the Lippmann-Schwinger equation, and has the form

[TABLE]

where $\Theta(\bf{r}^{\prime})$ is a Heaviside step function which is $1$ within the structure and [math] elsewhere. The first term on the RHS of Eq. (8) represents the electric field due to the free current source in a homogeneous background,

[TABLE]

The second term in Eq. (8) captures the response of the system, and describes the field scattered by the structure. Due to the simple LHS of Eq. (7), both terms in Eq. (8) share the same simple dyadic Green’s function for a homogeneous medium,

[TABLE]

$\mathcal{G}_{0}$ is available analytically, being the scalar Green’s function for a homogeneous medium of premittivity $\varepsilon_{bg}$, which depends only on the dimension of the geometry investigated. Thus, $\mathbf{E}_{0}(\mathbf{r})$ in Eq. (8) is relatively straightforward to find.

The efficient solution of the second term of the Lippmann-Schwinger Eq. (8) is the key purpose of SDM. In quantum mechanics, the Lippmann-Schwinger equation is typically solved via either the Born series or simply the Born approximation. In electrodynamics, Eq. (8) is the basis of the volume integral method of moments [45], and the discrete dipole approximation [46, 47]. In contrast, SDM uses Eq. (8) as the basis of an analytic formulation, expanding the response of the structure by its eigenmodes. The eigenmodes are defined by the source-free version of the Lippmann-Schwinger equation Eq. (8),

[TABLE]

where the eigenvalue is

[TABLE]

Here, we choose to fix the frequency ($k_{0}$) and vary the inclusion permittivity until Eq. (11) is satisfied. Such eigenmodes are not the usual eigenmodes found in the literature, where the permittivity is fixed and frequency is taken as the eigenvalue, and instead correspond to a fixed frequency with the permittivity of the structure being the eigenvalue. The eigenvalue equation is thus defined by replacing $\epsilon_{in}$ of Eq. (8) by $\varepsilon_{l}$ in Eq. (12), the latter corresponding to the eigen-permittivity of the mode. Note that $\varepsilon_{l}$ differs from, and is independent of the actual permittivity of the structure/inclusion ($\varepsilon_{in}$), and the set of all modes each with a different $\varepsilon_{l}$ provides a basis applicable to all possible $\epsilon_{in}$. For complete details related to the evaluation of the eigenvalues and normalized eigenmodes, see Ref. [48].

Once the eigenmode basis has been obtained, Eq. (8) can be solved by projecting onto the set of eigenmodes of Eq. (11). After some manipulation  [48], this procedure yields

[TABLE]

where for convenience, we have introduced the Dirac notation of quantum mechanics. Equation (13) expresses the response of the system in terms of the eigenmodes, governed by two weights. First is the detuning of the inclusion’s actual permittivity, $\varepsilon_{in}$, from the modal eigen-permittivity, $\varepsilon_{l}$. For $\varepsilon_{l}=\varepsilon_{in}$, one encounters a resonance of the structure [49]. Close to this resonance, the response of the system (for well-separated resonances) is governed by the mode $|\mathbf{E}_{l}\rangle$, and the contribution from all the other modes and $|\mathbf{E}_{0}\rangle$ can be neglected. Second, $\langle\mathbf{E}_{l}|\hat{\theta}|\mathbf{E}_{0}\rangle$ determines the efficiency of excitation of each mode $|\mathbf{E}_{l}\rangle$ as a function of overlap of the specific shape of the illumination field with the mode $l$. Finally, Eq. (13) demonstrates an important advantage of SDM, since the $|\mathbf{E}\rangle$ produced by the current source can be obtained without further simulation once the eigenmodes are known.

The quasi-static limit of SDM was used in the past to demonstrate the nano-focusing properties of surface plasmons [50, 51, 52, 53, 54] and to calculate SHG from metal nano-particles [30, 31]. Here, we proceed a step further by applying SDM to describe surface plasmon assisted SHG in the electrodynamics limit within the so-called the undepleted pump approximation [55]. In our SHG calculations, the following variables $\{\mathbf{E}(\mathbf{r})$, $|\mathbf{E}_{0}\rangle$, $|\mathbf{E}_{l}\rangle$, $\varepsilon_{(in,bg,l)}$, $\omega,\leavevmode\nobreak\ k_{0}\}$ are replaced by $\{\mathbf{E}^{2\omega}(\mathbf{r})$, $|\mathbf{E}_{0}^{2\omega}\rangle$, $|\mathbf{E}_{l}^{2\omega}\rangle$, $\varepsilon_{(in,bg,l)}^{2\omega}$, $2\omega,\leavevmode\nobreak\ 2k_{0}\}$. The free current ($\mathbf{J}_{f}$) dependence of $|\mathbf{E}_{0}^{2\omega}\rangle$ in Eq. (9) is now given by $\mathbf{J}^{(2)}=-2i\omega\mathbf{P}^{(2)}$. It is then straightforward to evaluate $|\mathbf{E}_{0}^{2\omega}\rangle$ from Eq. (9) using the free space Green’s function. Our main observation is that, since the eigenmodes are obtained for a source-free configuration, the tangential components of the eigenmodes, i.e., $\langle\mathbf{E}_{l,\parallel}^{2\omega}|$, are continuous across the interface. Thus, it follows from Eq. (13) that the discontinuity of $\mathbf{E}_{\parallel}^{2\omega}$ across the interface for perpendicular surface current sources can arise only from $|\mathbf{E}_{0}^{2\omega}\rangle$. In that sense, in SDM, the parameter $\varepsilon^{\prime}$ in Eq. (3) is naturally and unambiguously set to the background dielectric function ($\varepsilon^{2\omega}_{bg}$) by construction (without invoking the derivation of Sipe et al. [10]).

In the following Section, we demonstrate the power of SDM by considering SHG from the simple geometry of a single metal cylinder.

4 SHG from Single cylinder with surface nonlinearity: SDM implementation

We study SHG arising due to surface nonlinearity from an infinitely long metallic cylindrical inclusion (of permittivity $\varepsilon_{in}=\varepsilon_{cyl}$) of radius $a$ placed in a dielectric background. In this geometry, $(\mathbf{r}_{\parallel},{r}_{\perp})$ take the form $(\theta,\rho)$ corresponding to cylindrical coordinates in two dimensions. We assume a normally incident ($\partial_{z}=0$) in-plane TE polarized plane wave at FF. The total FF fields generated by the cylinder serve as the source for the SHG. The surface nonlinear polarization $\mathbf{P}^{(2)}_{S}$ is now placed at the interface $\rho=a$. We use vector cylindrical harmonics to project $\mathbf{J}^{(2)}$ on onto various angular orders $m$, giving

[TABLE]

In writing Eq. (14), we chose the FF fields from the metal, in accordance with Ref. [10]. Once $\mathbf{J}_{m}^{(2)}$ is evaluated, it is then straightforward to calculate $|\mathbf{E}_{0}^{2\omega}\rangle$ from Eq. (9). The free space scalar Green’s function in two-dimensions is given by the Hankel function of first kind $\mathcal{H}^{(1)}_{0}(k_{bg}^{2\omega}|\rho-\rho^{\prime}|)$ with $k_{bg}^{2\omega}=2\sqrt{\varepsilon_{bg}^{2\omega}}k_{0}$. We use Graf’s addition theorem to shift $\mathcal{H}^{(1)}_{0}$ to the coordinate origin [56]. The analytical solution of $|\mathbf{E}_{0}^{2\omega}\rangle$ is thus

[TABLE]

where $\mathcal{J}_{m}$ ($\mathcal{H}^{(1)}_{m}$) and $\mathcal{J}_{m}^{\prime}$ ($\mathcal{H}^{\prime(1)}_{m}$) correspond to the Bessel (Hankel) function of order $m$ and its derivative, respectively. The constant $A_{m}$ is given by

[TABLE]

Having evaluated $|\mathbf{E}_{0}^{2\omega}\rangle$, the eigenmodes and eigen-permittivities need to be found in order to complete the solution in Eq. (13). These are defined by Eq. (11), which corresponds to a cylinder at resonance. For cylindrical geometries however, it is not necessary to solve Eq. (11) directly, since solutions of the standard step-index fiber dispersion relation also produce eigen-permittivities that satisfy Eq. (11). This transcendental equation requires a root search in the complex plane, but we employed an efficient method based on contour integration, guaranteed to capture all the roots [57]. The corresponding modal fields are given in Ref. [48], which also describes the necessary normalization and projection steps to use Eq. (13), all performed analytically.

In what follows, we compare the SH electric field solutions as evaluated by SDM to the nonlinear Mie-type solution (invoking Eq. (3)) [21]. We consider a typical cylindrical silver inclusion of radius $a=30$ nm, whose permittivity is given by the Drude model and the background medium is considered to be air. Other parameters are provided in the caption of Fig. 1. We consider the dominant contribution arising from the angular orders $-4\leq m\leq 4$ for the incident wavelength $\lambda_{FF}=700$ nm. We plot the normalized $D_{\rho}^{2\omega}$ and $E_{\theta}^{2\omega}$ as evaluated by nonlinear Mie-type solutions in Fig. 1(a) and Fig. 1(b), respectively. We see that $D_{\rho}^{2\omega}$ is continuous, whereas $E_{\theta}^{2\omega}$ is discontinuous across the interface, as expected. The angular variation of $D_{\rho}^{2\omega}$ and $E_{\theta}^{2\omega}$ is dominated by the quadrapolar mode, i.e., $m=\pm 2$. The relative errors between nonlinear Mie-type and SDM solution for $D_{\rho}^{2\omega}$ and $E_{\theta}^{2\omega}$ are shown in Fig. 1(c) and Fig. 1(d), respectively, showing excellent agreement between both solutions. It is important to note that this excellent agreement is achieved by incorporating only the first radial order (i.e., $l=1$ [58]) in computing $|\mathbf{E}^{2\omega}\rangle$ in Eq. (13) for the angular modes $-4\leq m\leq 4$. Given the frequency of operation, the response of a generic plasmonic sub wavelength structure is dominated by the combination $|\mathbf{E}_{0}^{2\omega}\rangle$ and the first radial mode, which is usually plamonic in nature [48]. Therefore, the first radial order would suffice to achieve such excellent agreement.

It follows from Eqs. (15)-(16) that the normal component of $|\mathbf{E}_{0}^{2\omega}\rangle$, ${E}_{0,\rho}^{2\omega}$, is continuous across the interface, whereas the tangential component, $E_{0,\theta}^{2\omega}$ is discontinuous. Specifically, the $\mathbf{E}_{\parallel}^{2\omega}$-discontinuity is given by [59]

[TABLE]

Since the tangential components of the eigenmodes are continuous (as the eigenmodes are obtained from a source-free configuration), the discontinuity $\Delta{E}_{\theta}^{2\omega}$ arises only from the free space contribution $\Delta{E}_{0,\theta}^{2\omega}$. Importantly, even though we did not employ surface-source terms in the boundary conditions, $\Delta{E}_{0,\theta}^{2\omega}$ [Eq. (17)] has exactly the same form as in Eq. (3) with $\varepsilon^{\prime}=\varepsilon^{2\omega}_{bg}$.

5 The spectral dependence of discontinuity in ${E}_{\parallel}^{2\omega}$

We now study the spectral dependence of the discontinuity in the tangential component of electric field $\Delta{E}_{\theta}^{2\omega}$ and exploit SDM to provide physical insight into it. Specifically, for a given angular order $m$, it can be seen from Eq. (13) that close to the resonance of the first radial order ($\varepsilon^{2\omega}_{l=1,m}\rightarrow\varepsilon^{2\omega}_{cyl}$), the first weight factor attains a large value. Hence, the mode $|\mathbf{E}_{l=1,m}^{2\omega}\rangle$ contributes dominantly to $|\mathbf{E}^{2\omega}\rangle$, while the contributions from $|\mathbf{E}_{0}^{2\omega}\rangle$ and other radial modes are negligible. Thus, close to this resonance, the discontinuity $\Delta{E}_{\theta}^{2\omega}$ will be minimal. On the other hand, away from this resonance, the $|\mathbf{E}_{0}^{2\omega}\rangle$ contribution is significant, hence $\Delta{E}_{\theta}^{2\omega}$ will be significant. Such a simple physical insight is not provided by the nonlinear Mie-type solution [21] and this is crucial in understanding the importance of the errors associated with the incorrect implementation of the generalized BC (3). In order to quantify this trend, we define $\kappa_{m}$ as the ratio of the difference and the sum of $E_{m,\theta}^{2\omega}$ across the interface, namely,

[TABLE]

where $\rho=a\pm\delta a$ correspond to a point outside/inside the cylinder with $\delta a/a\ll 1$. Note that unlike $\Delta{E}_{\theta}^{2\omega}$, which depends on $\theta$, $\kappa_{m}$ does not depend on $\theta$ as it is defined for a specific angular order.

Figure 2(a) plots $\kappa_{m}$ as function of incident wavelength $\lambda_{FF}$ for a typical silver cylinder with a=30\leavevmode\nobreak\nm for $-4\leq m\leq 4$. In order to gain additional physical insight, we plot the first weight of Eq. (13), i.e., $|\varepsilon_{l=1,m}^{2\omega}-\varepsilon_{cyl}^{2\omega}|$ in Fig. 2(b). The location of the minimum value attained by $|\varepsilon_{l=1,m}^{2\omega}-\varepsilon_{cyl}^{2\omega}|$ defines the resonance of that particular angular order $m$. For example, the quadrapolar mode ($m=\pm 2$) is resonant around $\lambda_{FF}=667$ nm and exactly at the same location $\kappa_{m=\pm 2}$ attain a minimum value of about 2.5%. $|\varepsilon_{l=1,m}^{2\omega}-\varepsilon_{cyl}^{2\omega}|$ also determines the quality factor of the resonance, i.e., the smaller the minimum value of $|\varepsilon_{l=1,m}^{2\omega}-\varepsilon_{cyl}^{2\omega}|$, the sharper the resonance. Importantly, away from the resonance, $\kappa_{m}$ takes large values ($>$50%), meaning that the discontinuity becomes comparable to the field itself.

Since different angular orders $m$ have different weights in computing $E_{\theta}^{2\omega}$, it is not possible to determine the relative magnitude of $\Delta{E}_{\theta}^{2\omega}$ from $\kappa_{m}$ alone. Thus, we also define $\kappa$ as

[TABLE]

and plot its spectral dependence in Fig. 2(c). We see that close to resonance, i.e., for \lambda_{FF}=665\leavevmode\nobreak\nm, $\kappa$ attains a minimum value of about 2.5%, but away from resonance it attains values higher than 50%. For very large values of $\lambda_{FF}$, both $\kappa$ and $\kappa_{m}$ attain asymptotic values (not shown in the plot).

Finally, we compare the magnitude of the discontinuity for different materials. We now consider a gold cylinder (see Fig. 3) and the permittivity of the gold is obtained from the experimental data [60]. As before, we plot three quantities, namely, $\kappa_{m}$, $|\varepsilon_{l=1,m}^{2\omega}-\varepsilon_{cyl}^{2\omega}|$ and $\kappa$ in Fig. 3 and these quantities exhibit similar trends. However, the locations of the resonances and the minimum values attained by various quantities are substantially different, as expected. Indeed it is well known that gold is more lossy than silver. As a consequence, the minimum values attained by all the three quantities at the resonances for a gold inclusion are much higher than that of the silver inclusion. The minimum value attained by $\kappa$ for the gold inclusion is about 67%, implying that the $\mathbf{E}^{2\omega}_{\parallel}$-discontinuity is significant throughout the spectrum.

6 Numerical implementation

Many studies have implemented the generalized BC in SH calculations using the nonlinear Mie-type solution [21], the weak form of the differential equations [16], surface integral methods [22, 35, 23] and various other methods [27, 28, 61, 62].

We now discuss the implementation of the generalized BC in the commercial numerical package COMSOL Multiphysics [63]. Having solved for FF fields using the scattering wave formulation, $J_{S,\perp}^{(2)}$ can be constructed following Eq. (6). However, the default surface current settings in the radio frequency (RF) module do not permit the surface current to have a normal component, so the effect of $J_{S,\perp}^{(2)}$(in general, any given $J_{S,\perp}$) is not accounted for [64]. Accordingly, modeling any surface nonlinearity (such as the models proposed in Refs. [33, 34, 32]) with a perpendicular component (i.e., including models more general than the one used here Eq. (6)) based on COMSOL’s existing settings is necessarily only partially correct.

Instead of the default settings, one can choose to work with the weak form of the differential equations or redefine the boundary conditions. An alternative method to implement the generalized BC in Eq. (3) is to approximate the surface with a very thin layer and use the external bulk current sources option. Note that the thin layer approximation is a numerical implementation of Heinz’s derivation [37]. Following Sipe et al., the FF fields (needed for the calculation of $\mathbf{P}^{(2)}_{S}$) have to be chosen from the metal side of the interface and the thin layer (of thickness $t$) should have the permittivity of the background, i.e., $\varepsilon^{\prime}=\varepsilon_{bg}^{2\omega}$ [41]. Hence, this implementation results in the cylinder radius at SH being $a-t$. Thus, while this approach can be applied to most structures, it would fail to mimic the correct physics of complicated structures, such as kissing cylinders [65], as the structures at SH will always be non-touching, thus, modifying the response significantly [66]. Moreover, this implementation can be computationally expensive as the thin layer requires extremely fine numerical discretization.

In order to circumvent this limitation, and as a simpler alternative to adding the thin layer, we exploit the physical interpretation of $|\mathbf{E}_{0}^{2\omega}\rangle$ from the SDM. We use $|\mathbf{E}_{0}^{2\omega}\rangle$ [see Eqs. (9)] generated by $J_{S,\perp}^{(2)}$ as our background volume incident field in the scattering field formulation instead of the surface source $J_{S,\perp}^{(2)}$ itself, see Appendix B. For the single cylinder case, we have used the analytical expressions of $|\textbf{E}_{0}^{2\omega}\rangle$ [see Eqs. (15)-(16)] in the scattered wave formulation. We have tested the numerical convergence and find the field distributions to be in excellent agreement with the nonlinear Mie-type solution (not shown). Furthermore, we note that when it is not possible to obtain an analytical expression for $|\mathbf{E}_{0}^{2\omega}\rangle$ for an arbitrary structure, one may rely on surface integral methods to evaluate $|\mathbf{E}_{0}^{2\omega}\rangle$ [see Eq. (9)] arising from the surface current sources in free space/homogeneous background [67, 22].

7 Discussion and outlook

In order to optimize the SHG efficiency from nano plasmonic structures, one needs to make sure that both FF and SH frequencies are tuned to different resonances supported by the structure. Our results show that if only the excitation (FF) is resonant, then, an incorrect implementation of the generalized BC (3) can lead to severe errors in the computed SH near-fields. On the other hand, if the SH emission is resonant, then, the $\mathbf{E}_{\parallel}^{2\omega}$-discontinuity could be small (for instance, see Fig. 2) and thus, even if the $\mathbf{E}_{\parallel}^{2\omega}$-discontinuity is not implemented via Eq. (3), the parallel component of electric field can still be reasonably accurate. However, the correct implementation of the generalized BC is crucial throughout the spectrum for gold inclusion. The effect of an incorrect implementation of the BC (e.g., a la COMSOL) on far-field quantities, such as the absorption and the scattering cross-sections, will depend sensitively on the type of nonlinearity tensor and geometry; such a study lies beyond the scope of the current paper.

Finally, it is clear that the effects demonstrated here for SHG from a single metal wire with a simplistic second-order polarization tensor apply also for other structures (e.g., 3D nanostructures, waveguides etc. [20]), as well as for more sophisticated second-order optical effects, such as three wave mixing [68, 69]. It is very easy to incorporate the models for SH sources in SDM which are more general than the one used in Eq. (6), i.e., including the other components of $\chi^{(2)}_{S}$, nonlocal bulk effects, etc. We hope that future studies will indeed adopt the formulation described here.

Appendix A Sources

We now obtain the relations between the surface and volume charge accumulations generated by the nonlinear surface polarization $\mathbf{P}^{(2)}_{S}$. The surface ($\sigma^{(2)}$) and volume ($\rho^{(2)}$) charge are related via $\int\sigma^{(2)}dS=\int\rho^{(2)}dV$. In order to relate them with a differential relation, we use the usual continuity relation $\partial_{t}\rho^{(2)}=-\nabla\cdot\mathbf{J}^{(2)}$ (with current source $\textbf{J}^{(2)}$ as the external source that arises from SH polarization, $\mathbf{J}^{(2)}=-2i\omega\mathbf{P}^{(2)}$). $\sigma^{(2)}$ and $\rho^{(2)}$ in terms of $\mathbf{P}^{(2)}$ are given by

[TABLE]

respectively. Note that Eq. (20) is very general, however, it does not hold for the surface cases such as Eq. (1). We now derive the surface $\sigma^{(2)}$ and the volume $\rho^{(2)}$ charge accumulations arising from such surface sources. Integrating $\rho^{(2)}$ in Eq. (20) over the volume and invoking Eq. (1) gives us

[TABLE]

Note that in arriving at Eq. (22) from Eq. (21) we have employed the relation $\int\partial_{\perp}\delta({r}_{\perp})d{r}_{\perp}=0$. Thus, necessarily, the surface and volume charges for the surface sources are defined by $\varrho^{(2)}$ and $\varsigma^{(2)}$, respectively, are given by

[TABLE]

It is important to note that Eq. (20) becomes Eq. (24) for the surface sources.

Appendix B Scattered wave formulation for a current source

Since the existing COMSOL settings do not account for $J_{S,\perp}$, we would like to employ the insights provided by SDM to overcome this problem. Specifically, we now show how to derive $|\mathbf{E}_{0}^{2\omega}\rangle$ from $J_{S,\perp}$, and then use it in the scattered wave formulation.

Helmholtz’s equation with a free current source $J_{S,\perp}$ reads as

[TABLE]

Now $\mathbf{E}^{2\omega}$ can be decomposed as

[TABLE]

where $\mathbf{E}_{bg}^{2\omega}$ and $\mathbf{E}_{sc}^{2\omega}$ correspond to the field originating from the current source ${J}_{S,\perp}^{(2)}$ in an homogeneous background and the field scattered from the structure, respectively. Thus, $\mathbf{E}_{bg}^{2\omega}$ obeys

[TABLE]

Note that the above equation in its integral form is the same as Eq. (9) when written at SH, thus, $\mathbf{E}_{bg}^{2\omega}$ is identical to $|\mathbf{E}_{0}^{2\omega}\rangle$. Now, comparing Eq. (13) to Eq. (26) reveals that the scattered field $\mathbf{E}_{sc}^{2\omega}$ corresponds to the second term of RHS in Eq. (13). Rewriting Eq. (25) in terms of $\mathbf{E}_{sc}^{2\omega}$ by invoking Eqs. (26)-(27) and replacing $\mathbf{E}_{bg}^{2\omega}$ by $|\mathbf{E}_{0}^{2\omega}\rangle$ yields

[TABLE]

Eq. (28) is nothing but what is known as the scattered wave formulation, so that now, $\mathbf{E}_{sc}^{2\omega}$ can be evaluated from $|\mathbf{E}_{0}^{2\omega}\rangle$ instead of $J_{S,\perp}^{(2)}$.

Funding Information

KNR, PYC and YS were partially supported by Israel Science Foundation (ISF) (899/16). YS acknowledges the financial support from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007-2013) under REA grant (333790) and the Israeli National Nanotechnology Initiative. AIFD acknowledges funding from EU Seventh Framework Programme under Grant Agreement FP7-PEOPLE-2013-CIG-630996, and the Spanish MINECO under contract FIS2015-64951-R.