## Abstract

We present the frequency-domain boundary element formulation for solving surface second-harmonic generation from nanoparticles of virtually arbitrary shape and material. We use the Rao-Wilton-Glisson basis functions and Galerkin’s testing, which leads to very accurate solutions for both near and far fields. This is verified by a comparison to a solution obtained via multipole expansion for the case of a spherical particle. The frequency-domain formulation allows the use of experimentally measured linear and nonlinear material parameters or the use of parameters obtained using ab-initio principles. As an example, the method is applied to a non-centrosymmetric L-shaped gold nanoparticle to illustrate the formation of surface nonlinear polarization and the second-harmonic radiation properties of the particle. This method provides a theoretically well-founded approach for modelling nonlinear optical phenomena in nanoparticles.

© 2011 OSA

## Corrections

Jouni Mäkitalo, Saku Suuriniemi, and Martti Kauranen, "Boundary element method for surface nonlinear optics of nanoparticles: erratum," Opt. Express**21**, 10205-10206 (2013)

https://www.osapublishing.org/oe/abstract.cfm?uri=oe-21-8-10205

## 1. Introduction

Second-harmonic generation (SHG) is a nonlinear optical phenomenon, in which a field frequency component oscillating at double the frequency of the exciting field is generated. In the electric dipole approximation of material response, SHG vanishes in the bulk of materials with inversion symmetry [1]. At the interface, this symmetry is broken, which gives rise to surface SHG, which is very sensitive tool to probe surfaces.

Bulk SHG is allowed in centrosymmetric media if magnetic dipole and electric quadrupole responses are considered. The surface and bulk contributions to SHG are not fully separable in the sense that part of the bulk response behaves as the surface contribution [2]. This *inseparable contribution* can be included in the surface response, which is then described by an effective surface susceptibility. In general the surface and bulk responses of centrosymmetric media can be of similar magnitude. However, theoretical considerations estimate that for materials with high permittivity, such as metals, the surface contribution can be an order of magnitude greater than bulk contribution [3]. Also, experimental measurements suggest that the effective surface suceptibility is sufficient to describe the total SHG from gold structures [4].

The development of nanofabrication processes has enabled the study of nanoscale structures and their optical properties. From subwavelength structures, it is possible to construct metamaterials with tailored optical properties not found in nature. These properties depend on the size, shape and material of the particles, the properties of the surrounding medium [5] and the possible array structure [6]. Metal nanostructures exhibit plasmon resonances, which result in high local fields near particle surfaces. This can further amplify any nonlinear processes. Because SHG is sensitive to inversion symmetry, the observable second-harmonic (SH) far field also depends on the symmetry properties of the particles and the particle array.

Linear scattering problems of plasmonic nanoparticles have been solved numerically by a variety of methods including the Finite-Difference Time-Domain method, the Finite-Element Method, the Discrete Dipole Approximation, the Fourier Modal Method, the Volume Integral Method and the Surface Integral Method/Boundary Element Method. For many of these methods, problems arise at resonant conditions, as they make the problem sensitive and cause rapid spatial variations of fields, which requires careful discretization. As the linear plasmonic scattering problems can already be problematic, solving nonlinear plasmonic problems has been restricted to special cases.

Solving SHG in multilayer planar geometries is well-established as it can be done essentially in closed form [7]. Several numerical methods for modelling SHG and other nonlinear phenomena in nanoparticles of different geometries have been proposed. SHG and sum-frequency generation have been studied for particles of, in some sense, arbitrary shape by assuming vanishing or low refractive index contrast between the particle and the surrounding medium [8, 9]. In these approaches, the nonlinear scattering amplitude is deduced indirectly by employing the Lorentz reciprocity theorem. This theorem has also been applied to SHG from spheres of arbitrary shape and material [10]. SHG from spherical particles has also been modelled directly by multipole expansion of the nonlinear sources both for small and arbitrary size parameters [11–13]. SHG from 2-D cylindrical structures has been studied for time-harmonic and broad-band pulse excitation by employing the multiple scattering matrix method [14–16]. The Rayleigh equation method has been used for studying SHG from interfaces with translation or rotation symmetrical defects [17, 18]. Further, SHG and higher harmonic generation from layered periodic structures has been modelled by using the Fourier Modal Method (or Rigorous Coupled-Wave Theory) [19–21].

The aforementioned methods all either approximate the linear optical response or are based on a field expansion suitable to a very restricted class of geometries. There have been attempts to model SHG from arbitrarily shaped 3-D particles by employing numerical schemes directly to partial differential equations. These include the Finite-Difference Time-Domain method [22, 23] and the Finite Element Method [24]. In these methods one resorts to artificial absorbing boundaries to simulate an unbounded domain.

Although modelling nonlinear electromagnetic phenomena usually requires time-domain formulations, these present one significant drawback. They lead to cumbersome descriptions of the material dispersion and one is required to use simplified dispersion models to perform response convolutions. Thus it is difficult to use measured material parameters directly. Similar drawback applies to modelling the nonlinear response as one needs to use full time-domain models and ab-initio parameters even though only time-harmonic responses are often of interest. In most practical cases harmonic generation can be decomposed into coupled linear problems, so that frequency-domain methods can be used and thus measured material parameters can be directly utilized.

Recently, frequency-domain integral operator methods have gained popularity in the study of the linear response of nanoparticles. Depending on the electromagnetic properties of the media, these can be formulated as volume integral or surface integral methods. The finite-dimensional form of surface integral operator based scattering problems is usually referred to as the Boundary Element Method (BEM) or Surface Integral Equation method. BEM has been used for solving scattering from dielectric and ideally conducting bodies for a long time [25], but recently it has been applied also to scattering by plasmonic nanostructures [26–31].

The Volume Integral Method has been used for modelling SHG in 2-D nanostructures [32,33]. For the case of piece-wise homogeneous media it is possible to utilize BEM, where the unknowns of the problem are defined entirely on a compact boundary surface. This approach allows a natural introduction of surface nonlinear sources. Because in BEM only the fields on a surface need to be discretized, it is more scalable than volume element based methods. In certain formulations of BEM, the solution is not entirely unique, but can contain an unphysical contribution related to an isolated cavity resonator problem defined in the particle. This can happen especially near resonances in the physical solution and this can lead to instability. The Poggio-Miller-Chang-Harrington-Wu (PMCHW) formulation of BEM suppresses the unphysical solutions and thus guarantees stability at resonant conditions, which makes it ideal for plasmonic structures. The plasmon resonances are non-singular due to losses in metals. Integral operator based formulation also allows the use of arbitrary excitation sources, such as focused beams, which are important in the nonlinear microscopy of nanostructures.

In this work, we develop a BEM formulation for solving the surface SHG problem in lossy dielectric particles of virtually arbitrary shape. The formulation builds upon the undepleted-pump approximation, which allows the SH fields to be solved in three steps: 1) Solve a linear scattering problem at the fundamental frequency. 2) Determine the locally-varying source polarization at the second-harmonic frequency from the fundamental fields. 3) Solve a linear scattering problem at the SH frequency by using the polarization as a source. This treatment is valid in the case of low SHG conversion efficiency, which is the usual case for nanoparticles as verified by measurements. Note that although the problems are linear, as a whole they model a nonlinear phenomenon. The formulation allows an arbitrary excitation source, e.g., a plane-wave, focused Gaussian beam or an oscillating dipole. An advantage over reciprocity based approach is that the whole radiation pattern is obtained by solving a single scattering problem and also the SH near fields can be obtained. The integral operator formulation of the problem is developed in Section 2. We then present a finite-dimensional approximate representation in Section 3. For validation purposes, we also develop a multipole solution in Section 4, so that we may test our BEM method for the case of a spherical particle. In Section 5 we bring all this together and show comparison between BEM results and the multipole method and characterize the SH radiation properties of an L-shaped particle. We discuss the properties of the developed method in Section 6 and conclude in Section 7 with reference to future work.

## 2. Problem statement and integral operators

The solution domain of our SH scattering problem is illustrated in Fig. 1(a). The domain of the electric field **E** and magnetic field **H** is ℝ^{3} and it is divided into an unbounded exterior domain *V*_{1} and a compact interior domain *V*_{2}. The compact parts of the domain boundaries are denoted *∂V _{i}* and they are oriented such, that their respective normal vectors

**n**

*point out of*

_{i}*V*. The surface of the scatterer is denoted by

_{i}*S*and it is assumed to be piece-wise smooth and oriented with surface normal

**n**pointing into

*V*

_{1}. The linear electromagnetic properties of the domains are characterized by the permittivity

*ɛ*∈ ℂ and permeability

_{i}*μ*∈ ℝ, which also define the wave impedance

_{i}*η*= (

_{i}*μ*/

_{i}*ɛ*)

_{i}^{1/2}.

The SHG problem involves fields oscillating at two different frequencies: the fundamental fields **e*** _{i}* and

**h**

*at frequency*

_{i}*ω*and the SH fields

**E**

*and*

_{i}**H**

*at frequency Ω = 2*

_{i}*ω*with

*i*= 1,2 denoting the domains. We assume that the surface SH sources may be described in terms of a surface polarization distribution defined as

**P**

*=*

^{S}*ɛ*

_{0}

*χ*

^{(2)}:

**e**

_{2}

**e**

_{2}over the surface

*S*of the particle [3]. Of special interest are locally isotropic surfaces i.e. surfaces with local

*C*symmetry, so that the effective susceptibility tensor

_{∞}_{ν}*χ*

^{(2)}has only seven non-vanishing components: ${\chi}_{nnn}^{\left(2\right)}$, ${\chi}_{nss}^{\left(2\right)}={\chi}_{ntt}^{\left(2\right)}$, ${\chi}_{ssn}^{\left(2\right)}={\chi}_{sns}^{\left(2\right)}={\chi}_{ttn}^{\left(2\right)}={\chi}_{tnt}^{\left(2\right)}$, where

*n*refers to local normal and

*s,t*refer to the two mutually orthogonal tangent vectors.

Our formulation builds upon the *undepleted-pump* approximation, where the SH fields do not couple back to the fundamental field. This is justified by the fact that the measured SH signals are always orders of magnitude weaker than the source at the fundamental frequency. The problems at both frequencies are then linear and we can first solve the fundamental fields and then calculate the polarization **P*** ^{S}*, which acts as a source for the SH fields.

The time-dependence of harmonic fields is taken to be exp(−*iωt*). The fundamental fields in the domains *V _{i}* satisfy the Helmholtz equation ∇

*×*∇ ×

**f**–

*k*(

_{i}*ω*)

^{2}

**f**=

**0**, where

*k*(

_{i}*ω*)

^{2}=

*ω*

^{2}

*ɛ*(

_{i}*ω*)

*μ*(

_{i}*ω*) and

**f**∈ {

**e**,

**h**}. The expressions for the SH fields are the same with substitutions

**e**→

**E**,

**h**→

**H**and

*ω*→ Ω. In the SH problem, the surface polarization implies the following interface conditions over

*S*for tangential components of the SH fields [34]:

*ɛ*′ is called the

*selvedge region permittivity*[35]. This model assumes that the medium in the infinitely thin layer between

*∂V*

_{1}and

*∂V*

_{2}can be modelled macroscopically. In practice it may not be possible to determine this permittivity and one needs to resort to

*ad hoc*values. The problems of

*ɛ*′ are not in the scope of this work.

By assuming that the media in domains *V _{i}* are homogeneous, we may express the fields at any point by specifying the fields over the compact domain boundaries

*∂V*. This is governed by the Stratton-Chu equations [36]. To express the equations in our particular case, we define the Green’s function

_{i}*G*(

_{l}**r**,

**r**′) = exp(

*ik*)/(4

_{l}R*πR*) with

*R*= ||

**r**

*–*

**r**′||

_{2}and the following integro-differential linear operators:

*V*can now be expressed as

_{i}*V*

_{1}so that the fields represent purely outgoing waves [37]. The nonlinear source is not present in the integrals, because the source is introduced only in the interface conditions in Eqs. (1) and (2) and does not affect the representation of the fields in

*V*

_{1}and

*V*

_{2}.

It is customary to introduce the *equivalent* surface current densities, which are defined as
${\mathbf{J}}_{i}^{S}=\left(-{\mathbf{n}}_{i}\right)\times {\mathbf{H}}_{i}$ and
${\mathbf{M}}_{i}^{S}={\mathbf{E}}_{i}\times \left(-{\mathbf{n}}_{i}\right)$. Unlike in the linear problem, we now have four surface current densities to solve for due to the discontinuities in the tangential fields on *S*.

We take the limit where the boundaries *∂V _{i}* approach

*S*and enforce the interface conditions. By using the above representations of the fields, we obtain the infinite-dimensional problem: Given

**P**

*, seek ${\mathbf{J}}_{i}^{S}$, ${\mathbf{M}}_{i}^{S}:S\to {\u2102}^{3}$, such that Eqs.*

^{S}*S*. Note that we could substitute e.g. ${\mathbf{J}}_{2}^{S}$ and ${\mathbf{M}}_{2}^{S}$ from the Eqs. (9) and (10) into Eqs. (7) and (8), but this would result in numerically cumbersome integrals.

A similar formulation for the fundamental fields can be obtained by setting **P*** ^{S}* =

**0**and adding the tangential components of the incident fields to the right-hand-side of Eqs. (5) and (6) for

*V*

_{1}. This leads to the traditional formulation with only two unknown surface current densities as shown in e.g. [29]. Denoting the current densities in this case with lower case letters, it follows that

**j**

_{1}= −

**j**

_{2}and

**m**

_{1}= −

**m**

_{2}. The polarization

**P**

*is then evaluated from the field*

^{S}**e**

_{2}, whose components are given directly by the relations to corresponding equivalent surface current densities:

**m**

_{2}= −

**e**

_{2}×

**n**

_{2}and ∇

*·*

_{S}**j**

_{2}= −

*iωɛ*

_{2}

**n**

_{2}·

**e**

_{2}.

## 3. Finite dimensional formulation: the Method of Moments

To obtain approximate solutions to the SH problem, we employ the Method of Moments [38]. We seek such solutions from a finite-dimensional space that ensures proper continuity. A good choice for such a space is the one spanned by the Rao-Wilton-Glisson (RWG) functions, which are divergence-conforming affine functions that implicitly enforce the continuity of surface current and the conservation of charge [39]. These functions are geometrically attributed to the edges of a triangular mesh and their support, denoted by *S _{n}*, is a pair of adjacent triangles, as is illustrated in Fig. 1(b) (by the support of a function

*f*we mean the set

*S*= {

_{f}*x*∈

*S*|

*f*(

*x*) ≠ 0}). The unknown surface current densities are expanded using the RWG basis functions

**f**

*as*

_{n}*N*unknowns

*α*,

_{n}*β*,

_{n}*γ*,

_{n}*δ*∈ ℂ, that can be arranged into a vector

_{n}**x**= (

*α*

_{1},...,

*α*,

_{N}*β*

_{1},...,

*β*,

_{N}*γ*

_{1},...,

*γ*,

_{N}*δ*

_{1},...,

*δ*)

_{N}*. To obtain these coefficients, we use Galerkin’s testing [38] with the inner-product 〈*

^{T}*f,g*〉 = ∫

*·*

_{S}f*g*d

*S*.

The testing procedure leads to a linear system of equations, which can be expressed as Z**x** = **b**, where we have the system matrix Z:

*G*are evaluated at Ω. By enforcing both the electric and magnetic field interface conditions (the PMCHW testing), we avoid the internal resonance problem of BEM and thus also ensure robustness against plasmonic resonances [29].

_{l}The nonlinear surface polarization appears only in the source vector **b**:

**b**

^{n2}reduces to a contour integral by the identity

**n**is discontinuous, which implies that this identity should be applied piece-wise over each triangle. The curl of RWG-basis functions vanishes over their support and we obtain

*b*

^{n1}and

*b*

^{n2}clearly vanish if ${P}_{n}^{S}$ is constant, which is expected as the source in its original form depends on the surface gradient of ${P}_{n}^{S}$.

In case only the component
${P}_{n}^{S}$ is considered significant, we have **b**^{t1} = **b**^{t2} = **0** and Eq. (9) implies that
${\mathbf{J}}_{2}^{S}=-{\mathbf{J}}_{1}^{S}$. Then the problem is reduced to:

In the presented formulation, all the integrals can be evaluated with high precision by utilizing the singularity subtraction technique [40] and by using a high-order Gaussian quadrature for the resulting integrals with smooth kernels. In our calculations, a two-term singularity subtraction and 13-point Gauss-Legendre quadrature were used.

## 4. Solution in multipoles

To validate our BEM implementation, we develop an analytic solution for the same problem for the special case of a spherical particle by using the multipole expansion. This has been done previously in the small particle limit with the same interface conditions as here [13] and for an arbitrarily large sphere with different interface conditions [12]. An indirect far field solution has also been developed by using the Lorentz reciprocity [10].

A solution to the source-free Maxwell’s equations can be expanded in multipoles as [41]

**X**

*are related to the scalar spherical harmonics*

_{lm}*Y*via the angular momentum operator

_{lm}**L**by

**X**

*= [*

_{lm}*l*(

*l*+ 1)]

^{−1/2}

**L**

*Y*. Functions

_{lm}*f*and

_{l}*g*are linear combinations of spherical Bessel functions as discussed in [41]. The complex coefficients

_{l}*A*and

_{lm}*B*that fix the fields can be determined conveniently by using the orthogonality of

_{lm}**X**

*with respect to an inner-product defined in [41].*

_{lm}Next we consider a sphere of radius *a* and size parameter *x* = *k*_{1}*a*. We restrict ourselves here to the case of polarization
${\mathbf{P}}^{S}=2{\varepsilon}_{0}{\chi}_{nnn}^{\left(2\right)}{e}_{n}^{2}\mathbf{n}$, where

*B*-expansion coefficient of Eq. (21) of the fundamental problem in the interior domain. We use the identity ∇

*= −(*

_{S}*i*/

*r*

^{2})

**r**×

**L**to conveniently express the source function

**S**:

*ψ*(

_{l}*x*) =

*xj*(

_{l}*x*) and

*ξ*(

_{l}*x*)=

*xh*

^{(1)}(

*x*) and the relative refractive index

*N*=

*n*

_{2}(Ω)/

*n*

_{1}(Ω).

## 5. Numerical results

We next concentrate on modelling plasmonic nanostructures. First we validate BEM by comparing results with the ones obtained from the multipole expansion. We then apply the method for modelling SHG from an L-shaped particle, which has been studied experimentally before. In both cases the material is gold, whose refractive index is that from Johnson and Christy [42] and the surrounding medium is taken as vacuum. This yields also *ɛ*′ = *ɛ*_{0}.

#### 5.1. The spherical particle

We take a moderately sized spherical gold nanoparticle of radius *a* = 50 nm in vacuum, whence a plasmonic resonance takes place at the wavelength of *λ* = 520 nm. The excitation source is an *x*-polarized plane wave propagating in *z*-direction. The multipole solution of the fundamental field is sufficiently accurate with *l*_{max} = 4 (*m* = ±1) and the SH solution with *l*_{max} = 3 (all *m*-values included). We choose as the only nonzero tensor component
${\chi}_{nnn}^{\left(2\right)}=1$.

For the BEM we used two regularly triangulated icosahedral meshes, with 1280 and 5120 triangles and a more irregular mesh with 1454 triangles generated by GMSH’s MeshAdapt algorithm. The first mesh is shown in Fig. 1(c).

The computed radiated power per unit solid angle *σ* and the relative errors between the multipole and BEM solutions are shown in Fig. 2. The radiated power, defined by the far fields, is very accurate, the relative error being practically sub one percent. Surprisingly, the finest regular mesh yields the largest error, while the irregular mesh results in almost lowest errors in general. Thus the irregularity of the mesh does not deteriorate far field accuracy, but the maximum obtainable accuracy might be limited.

The SH electric field’s *x*-component amplitude near the sphere and its relative errors are shown in Fig. 3. Again the relative error is mostly below one percent, peaking at the point of highest field enhancement. Now the finest regular mesh yields markedly most accurate results and the irregular mesh is practically on par with the regular mesh of slightly lower triangle density. Overall the mesh refinement does not significantly remedy the high error at the point of highest enhancement.

We note that while the error in the far field quantity is largest for the densest mesh, it is already very low. Because the total size of the system is in the deep sub-wavelength regime, the far field can be very insensitive to small variations in the near field and the error cannot be expected to display a regular behaviour when the mesh is changed. The behaviour of the error is specific to each problem, this being just one particular case.

The normal component of the fundamental field *e _{n}* is represented in the method with piecewise constant functions, and the SH response depends on the surface gradient of
${e}_{n}^{2}$, which is only differentiable in the weak sense. Considering this, the results are surprisingly accurate, although it is quite well known that the BEM can produce highly accurate far fields even if using very coarse meshes.

#### 5.2. The L-shaped particle

The second-order nonlinear properties of L-shaped gold nanoparticles have been extensively studied experimentally (e.g. [43]). Although SHG from these particles has been measured, the simulations have been limited to modelling the linear response [44] and making indirect estimations of SHG. Here we apply our method to characterize the surface second-harmonic response of a single gold L-shaped nanoparticle in vacuum at plasmon resonances.

The L-shaped particle is illustrated in the inset of Fig. 4(a), where the chosen coordinate system is also shown. The particle’s arm length is 150 nm, arm width is 50 nm and the particle height is 20 nm. These dimensions are typical for the samples that have been used in measurements. We consider a plane wave propagating in *z*-direction, in which case the response of the particle is most clearly seen by considering *x*- and *y*-polarized incident waves. These lead to electrically antisymmetrical and symmetrical solutions to the scattering problem with respect to plane *x* = 0. The extinction spectra are shown in Fig. 4(a), with plasmon resonances at different wavelengths for the different incident polarizations. These resonances are related to charge oscillations along the arms and are sensitive to the arm length. Small peaks can also be observed at 530 nm for both polarizations. These are related to plasmon oscillations along the arm width. The extinction tail at even smaller wavelengths is mostly due to interband transitions of gold.

The effective second-order susceptibility components for gold have been measured from thin films, and their relative values are
${\chi}_{nnn}^{\left(2\right)}=250$,
${\chi}_{ntt}^{\left(2\right)}=1$ and
${\chi}_{ttn}^{\left(2\right)}=3.6$ when the polarization is evaluated using the fields inside the particle [4]. By using these values, the nonlinear surface polarization was computed by using the two polarizations for the incident wave at the corresponding resonance wavelengths. The results are plotted in Fig. 4(b). It is evident, that the nonlinear source polarization is driven into the corners of the particle. This suggests that only a small fraction of the particle surface gives rise to significant SHG, although the whole particle can affect the formation of the polarization in the first place. Small defects in these corners could drastically alter the SHG. The localization also implies that one must take care that the discretization of the problem is sufficiently smooth around sharp edges and corners. For this reason we used a mesh where the edges are carefully rounded as depicted in Fig. 1(d). As a last remark, the polarizations are practically symmetrical with respect to the *z* = 0 plane. It will be worthwhile to investigate if e.g. altering the height of the particle will induce phase retardation to the fundamental fields and give rise to less symmetrical source polarization.

The surface polarization does not directly tell how the system actually radiates second-harmonic waves. The full radiation patterns were computed with the developed method for the two incident polarizations at their resonance wavelengths. We attempt to gain insight into the importance of the different susceptibility components by calculating the radiated power by using the full susceptibility tensor and each tensor component separately. The results are shown in Fig. 5, where it can be seen, that at least in this case,
${\chi}_{nnn}^{\left(2\right)}$ clearly dominates the response. It is also clear that approximately the highest power per unit solid angle goes to forward and backward directions in the case of *x*-polarized input, but this is not the case for *y*-polarized input. The plots are symmetrical with respect to the plane *x* = 0, which is required for a valid solution to the problem. Symmetry considerations also dictate that SHG intensity in the forward and backward directions must be the same and this was verified to hold within 1 % relative error margin, except for the case of Fig. 5(g), for which the error was 2.8 %. The actual intensity in this case is, however, very weak compared to the other cases.

In the measurements, one is usually interested in the *x*- and *y*-components of the second-harmonic signal in the forward direction. Symmetry considerations dictate that only *y*-component can be nonzero in the case of an ideal particle. The fulfilment of this condition has been of considerable interest in the measurements, and it has been observed that small defects can easily brake this rule. For validation purposes, we made sure that our method gives rise to SHG that obeys this symmetry rule

## 6. Discussion

The developed method is applicable to a very general class of particle geometries: it is required that the particle is topologically homeomorphic with a sphere. The method is also applicable to a broad range of frequencies. This range is bounded by the low-frequency breakdown of BEM at very low frequencies and by the computational cost at very high frequencies. The frequency domain formulation allows the direct use of measured permittivity and second-order susceptibility values, which is a great advantage in scattering problems. The low refractive index contrast limit of [9] is obtained by making the Born approximation in Eqs. (5) and (6) i.e. by operating on incident fields to obtain an explicit solution.

It is also important to consider the computational burden of new numerical techniques. The presented BEM yields a dense, non-Hermitian system matrix. Thus, if we have *N* basis functions, the matrix build time and memory requirements scale proportionally to *N*^{2}. However, due to the block structure of the matrix, the matrix build time is not essentially different from the case of the linear scattering problem. If direct methods are used to solve the linear system of equations, the time complexity is of the order *O*(*N*^{3}). The general ”rule of thumb” is that at least ten basis functions are needed per wavelength, but sharp geometrical features may require locally denser mesh. However, BEM tends to yield high far field accuracy with very few basis functions [45]. Real samples usually consist of arrays of particles on a substrate. BEM can be conveniently generalized for modelling also this type of systems by using a periodic Green’s function, which can be efficiently evaluated by the Ewald’s method. It is also possible to significantly reduce the memory requirements and solution time by utilizing the Adaptive Cross Approximation or the Fast Multipole Method. The latter can reduce matrix-vector product time complexity to *N* log*N* in iterative solution methods.

As has been pointed out before, it is possible to use the Lorentz reciprocity for obtaining the SH scattered far field. Assume that we want to find out **u**^{*} · **E** (**u**^{*} denotes complex conjugate of **u**) far away for some unit vector **u**. We first solve the linear scattering problem at frequency *ω* for given excitation to obtain the nonlinear surface polarization **P*** ^{S}*. Then we solve another linear problem at frequency Ω, where the excitation is a plane wave incident from the observation direction with polarization

**u**yielding solution

**E**′. The Lorentz reciprocity [46] then states that

*χ*

^{(S)}and is at least four times less time consuming

*if only a single scattering direction is of interest*. If one desires to solve the whole radiation pattern, then BEM will be superior. This can be useful e.g. when we wish to simulate the SHG signal collected by an objective with large numerical aperture in nonlinear microscopy of nanostructures. Also, attempting to obtain the SH fields near the particle by using the reciprocity is not very convenient and even for the far fields, only relative scattering amplitude is obtained. As a final note, because in BEM one solves directly the fields on the boundary of the particle, the integral (28) can be evaluated conveniently in closed form in the case of RWG-basis.

## 7. Conclusion

We presented a BEM formulation for solving surface second-harmonic generation from particles of arbitrary shape and material. The comparison of the results to accurate multipole solutions for a sphere revealed that the developed method has potential in modelling SHG in complicated structures. Both near and far fields exhibit relative errors in the range 0.1 %–10 % when using a practically feasible number of basis functions.

We also characterized the SHG response of an L-shaped gold nanoparticle, whose second-order nonlinear properties have been studied experimentally. The calculations suggest that the nonlinear surface polarization is driven into the sharp edges of the particle and that the second-order nonlinear susceptibility component ${\chi}_{nnn}^{\left(2\right)}$ dominates the second-harmonic response as suggested by its large relative magnitude.

Although the treatment here was focused on surface SHG, it is in principle possible to extend the method for modelling bulk SHG that originates from higher microscopic multipoles. This would be done by considering the general Stratton-Chu equations with source volume current densities, which depend on gradients of the fundamental electric field. The accurate evaluation of these gradients near the particle surface is nontrivial since it requires proper treatment of hypersingular integral kernels and thus requires further analysis.

The treatment presented here could also be easily extended to modelling e.g. sum-frequency generation and higher harmonic generation. It is also possible to give up the undepleted-pump approximation and seek a solution to a fully coupled problem. This will give rise to a large nonlinear system of equations, so that it will become necessary to exploit geometrical symmetry, advanced matrix compression or possibly the Fast Multipole Method.

The BEM can also be extended to modelling spatially periodic structures by employing periodic Green’s functions, which can be rapidly evaluated by e.g. the Ewalds method. Also scattering problems consisting of multiple bodies of different media are straight-forward to implement in BEM. These additional features should be implemented before a realistic comparison to measurements can be carried out.

To conclude, the presented method enables accurate simulation of nonlinear phenomena in plasmonic nanostructures. This can be used in the design of new kinds of nanostructures and metamaterials with special nonlinear optical properties.

## Acknowledgments

JM acknowledges support from the Graduate School of Tampere University of Technology. We thank Stefan Kurz from the Department of Electronics of Tampere University of Technology for helpful discussions about theoretical details.

## References and links

**1. **Y. Shen, “Surface properties probed by second-harmonic and sum-frequency generation,” Nature **337**, 519–525 (1989). [CrossRef]

**2. **J. E. Sipe, V. Mizrahi, and G. I. Stegeman, “Fundamental difficulty in the use of second-harmonic generation as a strictly surface probe,” Phys. Rev. B **35**, 9091–9094 (1987). [CrossRef]

**3. **P. Guyot-Sionnest, W. Chen, and Y. Shen, “General considerations on optical second-harmonic generation from surfaces and interfaces,” Phys. Rev. B **33**, 8254 (1986). [CrossRef]

**4. **F. Wang, F. Rodríguez, W. Albers, R. Ahorinta, J. Sipe, and M. Kauranen, “Surface and bulk contributions to the second-order nonlinear optical response of a gold film,” Phys. Rev. B **80**, 233402 (2009). [CrossRef]

**5. **K. Kelly, E. Coronado, L. Zhao, and G. Schatz, “The optical properties of metal nanoparticles: the influence of size, shape, and dielectric environment,” J. Phys. Chem. B **107**, 668–677 (2003). [CrossRef]

**6. **H. Husu, J. Mäkitalo, R. Siikanen, G. Genty, H. Pietarinen, J. Lehtolahti, J. Laukkanen, M. Kuittinen, and M. Kauranen, “Spectral control in anisotropic resonance-domain metamaterials,” Opt. Lett. **36**, 2375–2377 (2011). [CrossRef] [PubMed]

**7. **J. Sipe, “New Green-function formalism for surface optics,” J. Opt. Soc. Am. B **4**, 481–489 (1987). [CrossRef]

**8. **S. Roke, M. Bonn, and A. Petukhov, “Nonlinear optical scattering: The concept of effective susceptibility,” Phys. Rev. B **70**, 115106 (2004). [CrossRef]

**9. **A. de Beer, S. Roke, and J. Dadap, “Theory of optical second-harmonic and sum-frequency scattering from arbitrarily shaped particles,” J. Opt. Soc. Am. B **28**, 1374–1384 (2011). [CrossRef]

**10. **A. de Beer and S. Roke, “Nonlinear Mie theory for second-harmonic and sum-frequency scattering,” Phys. Rev. B **79**, 155420 (2009). [CrossRef]

**11. **J. Dewitz, W. Hübner, and K. Bennemann, “Theory for nonlinear Mie-scattering from spherical metal clusters,” Zeitschrift für Physik D Atoms, Molecules and Clusters **37**, 75–84 (1996). [CrossRef] [PubMed]

**12. **Y. Pavlyukh and W. Hübner, “Nonlinear Mie scattering from spherical particles,” Phys. Rev. B **70**, 245434 (2004). [CrossRef]

**13. **J. Dadap, J. Shan, and T. Heinz, “Theory of optical second-harmonic generation from a sphere of centrosymmetric material: small-particle limit,” J. Opt. Soc. Am. B **21**, 1328–1347 (2004). [CrossRef]

**14. **C. Biris and N. Panoiu, “Second harmonic generation in metamaterials based on homogeneous centrosymmetric nanowires,” Phys. Rev. B **81**, 195102 (2010). [CrossRef]

**15. **C. Biris and N. Panoiu, “Nonlinear pulsed excitation of high-Q optical modes of plasmonic nanocavities,” Opt. Express **18**, 17165–17179 (2010). [CrossRef] [PubMed]

**16. **C. Biris and N. Panoiu, “Excitation of linear and nonlinear cavity modes upon interaction of femtosecond pulses with arrays of metallic nanowires,” Appl. Phys. A pp. 1–5 (2011).

**17. **L. Cao, N. Panoiu, and R. Osgood Jr, “Surface second-harmonic generation from surface plasmon waves scattered by metallic nanostructures,” Phys. Rev. B **75**, 205401 (2007). [CrossRef]

**18. **L. Cao, N. Panoiu, R. Bhat, and R. Osgood Jr, “Surface second-harmonic generation from scattering of surface plasmon polaritons from radially symmetric nanostructures,” Phys. Rev. B **79**, 235416 (2009). [CrossRef]

**19. **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). [CrossRef]

**20. **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]

**21. **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]

**22. **W. Schaich, “Second harmonic generation by periodically-structured metal surfaces,” Phys. Rev. B **78**, 195416 (2008).

**23. **Y. Zeng, W. Hoyer, J. Liu, S. Koch, and J. Moloney, “Classical theory for second-harmonic generation from metallic nanoparticles,” Phys. Rev. B **79**, 235109 (2009). [CrossRef]

**24. **G. Bachelier, I. Russier-Antoine, E. Benichou, C. Jonin, and P. Brevet, “Multipolar second-harmonic generation in noble metal nanoparticles,” J. Opt. Soc. Am. B **25**, 955–960 (2008). [CrossRef]

**25. **K. Umashankar, A. Taflove, and S. Rao, “Electromagnetic scattering by arbitrary shaped three-dimensional homogeneous lossy dielectric objects,” IEEE Trans. Antennas Propag. **34**, 758–766 (1986). [CrossRef]

**26. **J. Aizpurua, P. Hanarp, D. S. Sutherland, M. Käll, G. W. Bryant, and F. J. Garcia de Abajo, “Optical properties of gold nanorings,” Phys. Rev. Lett. **90**, 57401 (2003). [CrossRef]

**27. **I. Romero, J. Aizpurua, G. W. Bryant, and F. J. G. D. Abajo, “Plasmons in nearly touching metallic nanoparticles: singular response in the limit of touching dimers,” Opt. Express **14**, 9988–9999 (2006). [CrossRef] [PubMed]

**28. **G. Bryant, F. de Abajo, and J. Aizpurua, “Mapping the plasmon resonances of metallic nanoantennas,” Nano Lett. **8**, 631–636 (2008). [CrossRef] [PubMed]

**29. **A. Kern and O. Martin, “Surface integral formulation for 3D simulations of plasmonic and high permittivity nanostructures,” J. Opt. Soc. Am. A **26**, 732–740 (2009). [CrossRef]

**30. **B. Gallinet and O. Martin, “Scattering on plasmonic nanostructures arrays modeled with a surface integral formulation,” Phot. Nano. Fund. Appl. **8**, 278–284 (2010). [CrossRef]

**31. **B. Gallinet, A. Kern, and O. Martin, “Accurate and versatile modeling of electromagnetic scattering on periodic nanostructures with a surface integral approach,” J. Opt. Soc. Am. A **27**, 2261–2271 (2010). [CrossRef]

**32. **A. Benedetti, M. Centini, C. Sibilia, and M. Bertolotti, “Engineering the second harmonic generation pattern from coupled gold nanowires,” J. Opt. Soc. Am. B **27**, 408–416 (2010). [CrossRef]

**33. **M. Centini, A. Benedetti, C. Sibilia, and M. Bertolotti, “Coupled 2D Ag nano-resonator chains for enhanced and spatially tailored second harmonic generation,” Opt. Express **19**, 8218–8232 (2011). [CrossRef] [PubMed]

**34. **T. F. Heinz, “Second-order nonlinear optical effects at surfaces and interfaces,” in *Nonlinear Surface Electromagnetic Phenomena*, H.-E. Ponath and G. I. Stegeman (Elsevier, Amsterdam, 1991) p. 353.

**35. **J. E. Sipe, V. C. Y. So, M. Fukui, and G. I. Stegeman, “Analysis of second-harmonic generation at metal surfaces,” Phys. Rev. B **21**, 4389 (1980). [CrossRef]

**36. **J. Stratton, *Electromagnetic theory* (New York and London: McGraw-Hill, 1941).

**37. **D. Colton and R. Kress, *Inverse acoustic and electromagnetic scattering theory*, vol. 93 (Springer Verlag, 1998).

**38. **R. Harrington, *Field computation by moment methods* (Wiley-IEEE Press, 1993). [CrossRef]

**39. **S. Rao, D. Wilton, and A. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” IEEE Trans. Antennas Propag. **30**, 409–418 (1982). [CrossRef]

**40. **I. Hänninen, M. Taskinen, and J. Sarvas, “Singularity subtraction integral formulae for surface integral equations with RWG, rooftop and hybrid basis functions,” Prog. Elec. Res. **63**, 243–278 (2006). [CrossRef]

**41. **J. Jackson, *Classical electrodynamics* (John Wiley & Sons inc., 1999).

**42. **P. Johnson and R. Christy, “Optical constants of the noble metals,” Phys. Rev. B **6**, 4370 (1972). [CrossRef]

**43. **S. Kujala, B. Canfield, M. Kauranen, Y. Svirko, and J. Turunen, “Multipole interference in the second-harmonic optical radiation from gold nanoparticles,” Phys. Rev. Lett. **98**, 167403 (2007). [CrossRef] [PubMed]

**44. **H. Husu, J. Mäkitalo, J. Laukkanen, M. Kuittinen, and M. Kauranen, “Particle plasmon resonances in L-shaped gold nanoparticles,” Opt. Express **18**, 16601–16606 (2010). [CrossRef] [PubMed]

**45. **A. F. Peterson, D. R. Wilton, and R. E. Jorgenson, “Variational nature of Galerkin and non-Galerkin moment method solutions,” IEEE Trans. Antennas Propag. **44**, 500–503 (1996). [CrossRef]

**46. **L. D. Landau and E. M. Lifshits, *The electrodynamics of continuous media* (Pergamon, Oxford, 1960).