Surface second-harmonic from metallic nanoparticle configurations - a transformation optics approach

TL;DR

This paper uses transformation optics to analytically study surface second-harmonic generation in metallic nanoparticle structures, revealing new physical insights and conditions affecting SHG efficiency.

Contribution

It introduces an analytical approach linking complex nanoparticle geometries to simpler structures, uncovering new phase-matching and geometric factors influencing SHG.

Findings

Phase-matching is relevant even at sub-wavelength scales.

A previously unidentified geometric factor affects SHG efficiency.

Analytical solutions are supported by numerical simulations.

Abstract

We study surface second-harmonic generation (SHG) from a singular plasmonic structure consisting of touching metallic wires. We use the technique of transformation optics and relate the structure to a rather simpler geometry, a slab waveguide. This allows us to obtain an analytical solution to the problem, revealing rich physical insights. We identify various conditions that govern the SHG efficiency. Importantly, our analysis demonstrates that apart from the mode-matching condition, phase-matching condition is relevant even for this sub-wavelength structure. Furthermore, we identify a geometric factor which was not identified before. We support our analysis with numerical simulations.

Equations

$$P_{S,\perp }^{2\omega }={\epsilon }_0{\chi }_{S,\perp \perp \perp }^{(2)}{E}_{\perp }^{\omega }{E}_{\perp }^{\omega }\delta (x^2+y^2\mp 2ax),$$

$$\Delta {\mathbf{E}}_{\parallel }^{2\omega }=-\mathbf{n}\times {\mathbf{J}}_{mS}^{2\omega },\text{with}{\mathbf{J}}_{mS}^{2\omega }=\frac{1}{{\epsilon }_{bg}^{2\omega }}\mathbf{n}\times \left({\nabla }_{\parallel }{P}_{S,\perp }^{2\omega }\right),$$

$$J_{z,r/l}(x,y)=\frac{{\chi }_{S,\perp \perp \perp }^{(2)}}{{\epsilon }_{bg}^{2\omega }}{\partial }_{\parallel }(E_{\perp }^{\omega }{E}_{\perp }^{\omega })\delta (x^2+y^2\mp 2ax),$$

$${\nabla }^2{H}_{z,TW}^{2\omega }=0.$$

$$\mathbf{n}\times ({\mathbf{H}}_{bg}-{\mathbf{H}}_m)=H_{z,TW}^{2\omega }{|}_{bg}-H_{z,TW}^{2\omega }{|}_m=0,$$

$$\left[\mathbf{n}\cdot \nabla \left(\frac{H_{z,TW}^{2\omega }{|}_{bg}}{{\epsilon }_{bg}^{2\omega }}-\frac{H_{z,TW}^{2\omega }{|}_m}{{\epsilon }_m^{2\omega }}\right)\right]=-2i\omega {\epsilon }_0{J}_{z,r/l}.$$

$$\left(\frac{x\pm a}{a}{\partial }_x+\frac{y}{a}{\partial }_y\right)\left[\frac{H_{z,TW}^{2\omega }{|}_{bg}}{{\epsilon }_{bg}^{2\omega }}-\frac{H_{z,TW}^{2\omega }{|}_m}{{\epsilon }_m^{2\omega }}\right]=-2i\omega {\epsilon }_0{J}_{z,r/l}.$$

$$\tilde{w}=\frac{g^2}{w^{\ast }}.$$

$$\tilde{x}=\frac{g^{2}x}{x^2+y^2},\tilde{y}=\frac{g^{2}y}{x^2+y^2}.$$

$$[{\partial }_{\tilde{x}}^2+{\partial }_{\tilde{y}}^2]H_{z,sl}^{2\omega }=0,$$

$$H_{z,sl}^{2\omega }{|}_{bg}=H_{z,sl}^{2\omega }{|}_m,$$

$$J_{z,r}(x,y)$$

$$J_{z,l}(x,y)$$

$${\alpha }^{\omega }=\frac{1}{2}\log \left(\frac{{\epsilon }_m^{\omega }-{\epsilon }_{bg}^{\omega }}{{\epsilon }_m^{\omega }+{\epsilon }_{bg}^{\omega }}\right)\text{when Re}({\epsilon }_m)<-{\epsilon }_{bg}^{\omega }.$$

$$J_{z,r/l}\to {\Delta }_{r/l}(\tilde{x},\tilde{y})e^{2i{\alpha }^{\omega }|\tilde{y}|/d}\delta (\tilde{x}=\pm d),$$

$$\left(\frac{d^2+{\tilde{y}}^2}{g^2}\right){\partial }_{\tilde{x}}\left[\frac{H_{z,sl}^{2\omega }{|}_{bg}}{{\epsilon }_{bg}^{2\omega }}-\frac{H_{z,sl}^{2\omega }{|}_m}{{\epsilon }_m^{2\omega }}\right]$$

$$H_{z,sl}^{2\omega }(\tilde{x},\tilde{y})=\{\begin{array}{l}{\mathcal{A}}_{sl}{e}^{2i{\alpha }^{\omega }|\tilde{y}|/d}\sinh \left(2{\alpha }^{\omega }\tilde{x}/d\right),|\tilde{x}|\le d,\\ {\mathcal{B}}_{sl}{e}^{2i{\alpha }^{\omega }|\tilde{y}|/d}\text{sgn}[\tilde{x}]e^{-2{\alpha }^{\omega }|\tilde{x}|/d},|\tilde{x}|>d.\end{array}$$

$${\mathcal{B}}_{sl}={\mathcal{A}}_{sl}{e}^{2{\alpha }^{\omega }}\sinh \left(2{\alpha }^{\omega }\right).$$

$${\mathcal{A}}_{sl}=\frac{-i\omega {\epsilon }_0{\epsilon }_{bg}^{2\omega }{g}^{2}d}{{\alpha }^{\omega }\mathcal{P}}\left(\frac{{\Delta }_r(\tilde{x}=d,\tilde{y})}{d^2+{\tilde{y}}^2}\right),$$

$$\mathcal{P}=\cosh \left(2{\alpha }^{\omega }\right)+\frac{{\epsilon }_{bg}^{2\omega }}{{\epsilon }_m^{2\omega }}\sinh \left(2{\alpha }^{\omega }\right).$$

$$\cosh {\alpha }^{2\omega }+\frac{{\epsilon }_{bg}^{2\omega }}{{\epsilon }_m^{2\omega }}\sinh {\alpha }^{2\omega }=0,$$

$$-{\epsilon }_{bg}^{2\omega }<\text{Re}({\epsilon }_m^{2\omega })<{\epsilon }_{bg}^{2\omega },$$

$$H_{z,TW}^{2\omega }=\{\begin{array}{l}{\mathcal{A}}_{TW}(x,y)\exp \left(\frac{4ia{\alpha }^{\omega }|y|}{x^2+y^2}\right)\sinh \left(\frac{4a{\alpha }^{\omega }x}{x^2+y^2}\right),\text{when}{x}^2+y^2+2a|x|\ge 0,\\ {\mathcal{B}}_{TW}(x,y)\exp \left(\frac{4ia{\alpha }^{\omega }|y|}{x^2+y^2}\right)\text{sgn}[x]\exp \left(\frac{-4a{\alpha }^{\omega }|x|}{x^2+y^2}\right),\text{when}{x}^2+y^2+2a|x|<0,\end{array}$$

$${\mathcal{A}}_{TW}(x,y)=\frac{-i\omega {\epsilon }_0{\epsilon }_{bg}^{2\omega }}{2{\alpha }^{\omega }\mathcal{P}}\mathcal{C}(x,y),$$

$${\mathcal{B}}_{TW}(x,y)={\mathcal{A}}_{TW}(x,y)e^{2{\alpha }^{\omega }}\sinh 2{\alpha }^{\omega },$$

$$\mathcal{C}(x,y)=\frac{4a(x^2+y^2{)}^2}{4a^2{y}^2+(x^2+y^2{)}^2}\mathcal{R}\left({\tau }_x,{\tau }_y\right),$$

$${\tau }_x=\frac{1/(2a)}{1/(4a^2)+y^2/(x^2+y^2{)}^2},$$

$${\tau }_y=\frac{y/(x^2+y^2)}{1/(4a^2)+y^2/(x^2+y^2{)}^2}.$$

$$H_{z,TW}^{2\omega }(x,y)=$$

$$H_{z,TW}^{2\omega }(x,y)=-i\omega {\epsilon }_0{\epsilon }_{bg}^{2\omega }\left(\frac{\sinh 2{\alpha }^{\omega }}{{\alpha }^{\omega }\mathcal{P}}\right)x{J}_{z,r}(x,y).$$

$${\partial }_x{H}_{z,TW}^{\omega }(x,y)$$

$${\partial }_y{H}_{z,TW}^{\omega }(x,y)$$

$$H_{z,TW}^{\omega }(x,y)=\{\begin{array}{l}\mathcal{N}\text{sgn}[y]\cosh \left(\frac{2a{\alpha }^{\omega }x}{x^2+y^2}\right)\exp \left(\frac{2ia{\alpha }^{\omega }|y|}{x^{\prime 2}+y^{\prime 2}}\right),\text{when}{x}^2+y^2+2|x|a>0,\\ \left(\frac{\mathcal{N}{\epsilon }_m^{\omega }}{{\epsilon }_m^{\omega }+{\epsilon }_{bg}^{\omega }}\right)\text{sgn}[y]\exp \left(\frac{-2a{\alpha }^{\omega }|x|}{x^2+y^2}\right)\exp \left(\frac{2ia{\alpha }^{\omega }|y|}{x^2+y^2}\right),\text{when}{x}^2+y^2+2|x|a<0,\end{array}$$

Full text

Suggested keywords

I Introduction

Transformation optics (TO) is a theoretical tool recently developed that allows an unprecedented control over light propagation and confinement [1, 2, 3]. TO relies on the form-invariance of Maxwell’s equations under coordinate transformations to provide the connection between a given electromagnetic effect, coded into a geometric transformation, and the material parameters ($\varepsilon$, $\mu$) required for its realization [4, 3].

Recently, it was shown that most of the complexity associated with TO schemes can be avoided by using a special class of transformations called conformal transformations (CT). These are two-dimensional (2D) mappings which have a very convenient characteristic - they preserve the material and spectral characteristics of the original system [5, 2]. This enables the description of the optical response of a wide variety of geometries by cascading conformal transformations [5, 6, 2, 7, 8]. The analytic predictions of TO have been found to be in excellent agreement with numerical results (for subwavelength particles), and verified experimentally by several groups [9, 10, 11, 12]. In addition, CT allowed addressing various other physical effects such as the moulding of surface plasmon polariton propagation [13], the emergence of non-local effects in metallic nanostructures [14, 15], van der Waals interactions at the nano-scale [16], graphene [17] and singular [18] metasurfaces, and plasmon-exciton interactions [19] and strong coupling [20].

The most fundamental contribution of CT is, perhaps, to establish the equivalence between the (linear) scattering from isolated particles (supporting localized plasmon resonances) and light propagation along straight waveguides (supporting surface plasmon resonances) [5, 7, 8, 21], at least within the quasi-static approximation. This enabled calculation of the field distribution in complex structures and especially the challenging singular geometries by first transforming to simple, regular solvable geometries, and then transforming back (see Figure 1). Such geometries give rise to extremely high local-field enhancement near the singularities. The CT also showed that the spectrum of the singular geometry is “inherited” from the broadband plasmonic waveguide structure, thus, resolving the long-standing arguments about the origin of the broadband spectral response from rough (disordered) and/or multiscale structures [5, 22]. This class of structures thus allow for breaking the conventional resonance-based limitations of bandwidth in plasmonic devices.

To date, all the studies employing TO were limited to media that have a linear response to the incoming electromagnetic field, namely, for low intensities 111The only exception is a paper that focused on the general formulation and provided a single example of third-order nonlinearity [69]. . However, extending TO to nonlinear wave interactions is appealing, since both the high local-field enhancement as well as the unusually wide bandwidth make the singular structures potentially useful for such interactions.

Here, we employ TO for a singular nanoparticle configuration that includes nonlinear media (specifically, for Second Harmonic Generation, SHG from metal nanoparticles). We demonstrate the strength of our approach via an analytic solution of the near-field distribution and conversion efficiency, which is found to be in excellent agreement with exact electrodynamic numerical simulations. Most importantly, we provide deep insights of the physics underlying frequency conversion processes.

In particular, we show that the theoretical description of nonlinear wave mixing from nanoparticles is more complex than one would have expected compared to the metal waveguides. Indeed it is well known that in the latter case, efficient SHG has three requirements - a strong source, phase matching (PM) and mode matching (MM) [24, 25]. Yet, for subwavelength structures, e.g., single nanoparticles or nanoparticle clusters, it is customary to dismiss the need to achieve phase matching [26], justified by the claim that the phase accumulation across the structure is small, so that the source and generated wave are nearly unaffected by any phase mismatch. This implies that there is a fundamental difference between extended structures (such as waveguides) and particle geometries. However, the intrinsic analogy between the infinitely extended flat geometry and particle(s) geometry (see Figure 1) [2] implies that this is not justified in all cases. Indeed, the CT map shows that the modes propagating an infinite distance in the flat geometry are the same modes that repeatedly circulate around the touching wires (TWs). Thus, since it is trivial that PM is essential in the flat geometry, then, it is necessarily as important in the particle geometry. Indeed, particle-induced wavelength compression makes the accumulated phase substantial despite the short propagation distance along the circumference and regardless of any detuning. In that sense, a critical observation is that for the purpose of frequency conversion, the particle geometry is very different from a bulk material or homogenized composite of the same dimensions for which the phase mismatch is not important.

Below, we show explicitly that PM does manifest itself in small NPs, and connect it to the localized plasmon resonance. Then, we show, somewhat unexpectedly, that these so-called three requirements are not sufficient, and that in practice there is an additional consideration. We discuss how this additional term modifies the SH response.

The article is organized as follows: Section II.1 describes the identical TW system, SH surface sources and their symmetry aspects and formulates the SHG problem in the TW frame. To obtain the SHG solution, we follow the route of conformal transformation to transform TW geometry to a rather simpler geometry, i.e., slab geometry and solve for SHG in Section II.2 . We then transform the obtained solution back to the TW geometry in Section II.3 and interpret it. Finally, we conclude with a discussion and outlook.

II Touching wire system and surface nonlinear sources

II.1 Formulation and methodology

Consider two identical metallic TWs of radii $a$ with permittivity $\varepsilon_{m}$, touching each other at the coordinate origin of the $x-y$ plane; they are embedded in a homogeneous dielectric background of permitivitty $\varepsilon_{bg}$ (see Figure 1(a)). The TW structure is assumed to be illuminated by a TM polarized plane wave at a frequency $\omega$ by a spatially uniform $x$-polarized electric field with amplitude $E^{\omega}_{0x}$. The linear electric field response of the structure at the fundamental frequency (FF), $\mathbf{E}^{\omega}$, can be evaluated analytically using the technique of conformal transformation as in [5, 27]. Firstly, the TW geometry and the plane wave source are transformed to a metal-dielectric-metal slab structure and a line dipole source, respectively. Then, the dipole radiation is coupled to plasmonic modes of the slab system in the momentum space. The dominant contribution to the linear response arises from the surface plasmon pole and the contribution from the lossy surface (radiative) waves was neglected. Having obtained the closed form solution in the slab frame, the solution was then transformed back to the TW frame. The analytic solution revealed several interesting physical phenomena occurring close to the singular (touching) point. Most notable of them are the unusually wide spectral response, strong spatial field confinement and large field enhancement close to the touching point. The latter effect can give rise to efficient nonlinear optical phenomena.

In this article, we consider a second-order process, specifically, second-harmonic generation (SHG). It is well known that the second-order nonlinear processes are symmetry forbidden in centro-symmetric materials. However, due to the broken symmetry at an interface, the metal-dielectric interface has a non-zero second-order surface tensor $\bar{\bar{\chi}}^{(2)}_{S}$ [28, 29]. Apart from purely surface effects, non-locality can give rise to nonlocal bulk nonlinearity which 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 [30, 31, 32, 33]. We assume that the surface SH polarization $\mathbf{P}^{2\omega}_{S}$ at the metal-dielectric interface is given by

[TABLE]

where $E_{\perp}^{\omega}$ and $P^{2\omega}_{S,\perp}$ correspond to the normal component (to the metal-dielectric interface) of linear electric field and surface SH polarization, respectively. The metal-dielectric interface of the right (left) wire corresponds to $-$ ($+$) sign in the argument of the Dirac-delta function $\delta$.

In writing $P^{2\omega}_{S,\perp}$, we assume that the dominant contribution to the surface nonlinear polarization arises from the $\chi^{(2)}_{S,\perp\perp\perp}$ element, see justification in Ref. [34]. To obtain $P^{2\omega}_{S,\perp}$, we use the analytical expressions of the FF electric field $E^{\omega}_{\perp}$ derived in Ref. [27]. We note that $E^{\omega}_{\perp}$ is discontinuous across the interface, thus, leading to an ambiguity in choosing $E^{\omega}_{\perp}$ either on the metal or dielectric side of the interface. Following Sipe et al. [35], we choose $E^{\omega}_{\perp}$ on the metal side of the interface.

In the solution of SH fields, Eq. (1) appears as a source term in Maxwell’s equations, and can be solved in this form [36, 37]. But in surface SHG, the sources coincide with the boundary, so alternatively $P^{2\omega}_{S,\perp}$ can be incorporated into the boundary conditions at the metal-dielectric interface [29]. We shall adopt this convenient approach.

First, note that the SH polarization $P^{2\omega}_{S,\perp}$ is normal to the interface, corresponding to an electric dipole layer pointing normal to the interface. In such a case, the tangential component of the SH field $\mathbf{E}_{\parallel}^{2\omega}$ is discontinuous. The discontinuity across the interface, denoted as $\Delta\mathbf{E}_{\parallel}^{2\omega}$, is given by the generalized boundary condition as $\Delta\mathbf{E}_{\parallel}^{2\omega}=-\frac{1}{\varepsilon_{bg}^{2\omega}}\nabla_{\parallel}P^{2\omega}_{S,\perp}$, where $\nabla_{\parallel}$ is the tangential derivative along the interface  [29, 36]. This generalized boundary condition can be conveniently reformulated using the divergence-free magnetic surface current density $\mathbf{J}_{mS}^{2\omega}$ [38, 36]. In such a case, the generalized boundary condition takes the form

[TABLE]

where $\mathbf{n}$ is the unit normal to the interface. The magnetic surface current density generated by the SH polarization $P^{2\omega}_{S,\perp}$ considered in Eq. (1) for the TW geometry is given by

[TABLE]

where $J_{z}$ is the out-of-plane (i.e., $z$) component and furthermore, the only component of the magnetic current $\mathbf{J}_{mS}^{2\omega}$ and we have dropped the superscript $2\omega$ and the subscript $mS$ of $J_{mS,z}^{2\omega}$ for the sake of brevity. The subscript $r$($l$) of $J_{z,r/l}$ denotes the right(left) wire and $\partial_{\parallel}$ denotes the tangential derivative along the interface.

Since our TW structure and the incidence is symmetric, we now investigate the symmetry relations of the $J_{z,r/l}$ on the left and right wires. Figure 2 shows various fields (on the metal side) on the circumference of the TW close to the touching point. As shown in Figure 2(a), the linear response $E_{\perp}^{\omega}$ has a definite symmetry which is anti-symmetric in $x$ and symmetric in $y$ below the surface plasma resonance frequency. Since $E_{\perp}^{\omega}E_{\perp}^{\omega}$ is quadratic in $E_{\perp}^{\omega}$, then necessarily $E_{\perp}^{\omega}E_{\perp}^{\omega}$ is symmetric in $x$ and $y$ (see Figure 2(b)), as expected. The tangential derivative $\partial_{\parallel}$ of $E_{\perp}^{\omega}E_{\perp}^{\omega}$ along the circumference of the left and right wires yields $J_{z,r/l}$ (3) which is an anti-symmetric source in $x$ and $y$. As a consequence, $J_{z,r}(x,y)=-J_{z,l}(-x,y)$ (see Figure 2(c)). Note that the SH source symmetry described above holds well for the choice of $\chi^{(2)}_{S,\perp\perp\perp}$ element at any frequency of operation and all material parameters.

Having understood the SH source and its symmetries, we now formulate the SH problem for the TW. Since the SH response is usually weak when compared to the FF response, we employ the undepleted pump approximation [39]. Furthermore, we also assume the physical dimensions of the wire $a$ to be much smaller than both the FF and SH wavelengths such that one can safely employ the quasi-static approximation (namely, we take the limit $(\omega/c)a\ll 1$). Since we consider in-plane $TM$ polarized illumination, it is convenient to formulate a scalar SHG problem in terms of the SH magnetic field $H_{z,TW}^{2\omega}$. We emphasize that for a quasi-static structure, the magnetic response is rather weak yet non-zero (see Appendix A for further discussion). Upon obtaining the SH magnetic field, the corresponding SH electric fields can then be evaluated using Ampere’s law.

The SH magnetic field $H_{z,TW}^{2\omega}$ can be obtained by solving the homogeneous Helmholtz equation in each domain accompanied with two boundary conditions. However, in the quasi-static limit Helmholtz equation in each domain simplifies to the homogeneous Laplace equation given by

[TABLE]

Futher, the two boundary conditions are as follows. Continuity of $H_{z,TW}^{2\omega}$ across the interface is given by

[TABLE]

where the subscripts $bg$ and $m$ correspond to the dielectric and metal sides of the interface, respectively. The generalized boundary condition (see Eq. (2)) in terms of $H_{z,TW}^{2\omega}$ and the magnetic surface density $J_{z,r/l}$ for the TW is given by

[TABLE]

Invoking the definition of $\mathbf{n}$ for the right and left wires, Eq. (6) can be rewritten as

[TABLE]

We solve the SH problem analytically using conformal transformation by relating the TW to a much simpler geometry, a slab waveguide, as for the linear case [5].

II.2 Transformation to the slab frame

Consider the conformal coordinate transformation

[TABLE]

Here, we use the complex number notation $w=x+iy$ for the TW coordinates and $\tilde{w}={\tilde{x}}+i{\tilde{y}}$ for the new coordinate system ($*$ denotes the complex conjugate). The parameter $g$ in Eq. (8) is a scaling constant with dimensions of length. The coordinates of the old and new systems are related by

[TABLE]

The touching wire system (in $x-y$ plane) transforms to a metal-dielectric-metal slab waveguide (in ${\tilde{x}}-{\tilde{y}}$ plane) of half-width $d=g^{2}/2a$ (see Figure 1). The permitivitty of the materials remains unchanged under this 2D conformal transformation [27].

Having transformed the TW system to slab frame, we now formulate the SH problem in the slab geometry by transforming the governing equation (4) and the two boundary conditions Eqs. (5) and (7). Firstly, we note that the magnetic field is invariant under 2D conformal transformation (see Appendix B) as it is an out-of-plane component [1, 40]. Therefore, $H_{z,TW}^{2\omega}(x,y)=H_{z,sl}^{2\omega}({\tilde{x}}(x,y),{\tilde{y}}(x,y))$). Now transforming Eqs. (4)-(5) and expressing them in terms of the slab variables, we arrive at

[TABLE]

[TABLE]

respectively.

Before transforming the boundary condition in Eq. (7), let us examine $J_{z,r/l}$ in the TW frame. We decompose $J_{z,r/l}$ in terms of amplitude and phase as

[TABLE]

where $\mathcal{R}$ and $\mathcal{L}$ correspond to the amplitude of the SH source on the right and left wires, respectively (see Appendix C). We now describe in detail the origin of the functional dependence of the phase part of $J_{z,r/l}$ in Eqs. (12).

Recall that the linear response was analytically obtained by transforming the TW geometry to the slab [5]. This slab structure in the quasi-static limit exhibits two surface plasmon modes. The symmetric mode extends from zero frequency to the surface plasmon resonance, while the anti-symmetric mode exists above the surface plasmon resonance. The linear solution involved only the symmetric mode whose dimensionless (i.e., the product of the propagation constant and slab half-width) propagation constant at FF is given by [41] 222Lei et al.[27] consider a slab of width $d$. Therefore, their dimensionless propagation constant differs by a factor 2 from ours. We adopt this change of notation in order to conform to standard PM notations, see below.

[TABLE]

Despite transforming the FF solution back to the TW frame, $\alpha^{\omega}$ still continues to act as the propagation constant for the TW system [43]. As a consequence, the phase of the linear field $E_{\perp}^{\omega}$ on the circumference of the identical TW is given by $\exp{\left(2ia\alpha^{\omega}|y|/(x^{2}+y^{2})\right)}$ [27, Eqs. (32)-(37)]. Since the SH source is obtained upon squaring the linear fields (see Eq. (3)), the phase of SH source is exactly twice the FF phase. Furthermore, as $J_{z,r}(x,y)=-J_{z,l}(-x,y)$, we have $\mathcal{R}(x,y)=-\mathcal{L}(-x,y)$.

The CT unfolds the considerably complicated SH source $J_{z,r/l}$ to a much simpler one in the slab frame. Specifically, transforming $J_{z,r/l}$ (12) to the slab frame yields (see Appendix D)

[TABLE]

where $\Delta_{r}$ and $\Delta_{l}$ are the amplitudes of the magnetic surface currents in the slab geometry obtained by transforming $\mathcal{R}$ and $\mathcal{L}$, respectively. Thus, $\Delta_{r}({\tilde{x}}=d,{\tilde{y}})=-\Delta_{l}({\tilde{x}}=-d,{\tilde{y}})$ as $\mathcal{R}(x,y)=-\mathcal{L}(-x,y)$, the phase part of $J_{z,r/l}$ (see Eqs. (12)) simplifies to $e^{2i\alpha^{\omega}|{\tilde{y}}|/d}$ with $2\alpha^{\omega}$ as dimensionless propagation constant of the SH source in the slab frame, and the SH surface source $J_{z,r}$ ($J_{z,l}$) on the circumference of the right (left) wire transforms to line source placed at the right (left) interface of the slab.

The generalized boundary condition (7) now transforms to the slab frame as (see Appendix D)

[TABLE]

We note that the spatial factor $(d^{2}+{\tilde{y}}^{2})$ on the left-hand-side of Eqs. (15) arises from transforming $\mathbf{n}\cdot\nabla$; this term thus encodes information of the original TW geometry (see Eqs. (6)-(7)).

We can now calculate the SHG response in the slab geometry by solving Eq. (10) with boundary conditions in Eq. (11) and Eqs. (15) using an ansatz. The mathematical form of our ansatz for $H_{z,sl}^{2\omega}$ is dictated by the governing physical phenomena, namely,

The ansatz should satisfy the Laplace equation (10) in all the domains. 2. 2.

Since we are interested in the near-field of the nanostructure which is expected to be governed by surface plasmon waves rather than photonic waveguide modes, the ansatz should exhibit exponential decay along the transverse coordinate ${\tilde{x}}$ away from the interfaces at ${\tilde{x}}=\pm d$. 3. 3.

Since the slab structure and SH source have definite symmetry in the ${\tilde{x}}-{\tilde{y}}$ plane, the ansatz should also have definite a symmetry. In fact, the anti-symmetric magnetic line sources at the interfaces ($\Delta_{r}=-\Delta_{l}$) dictates that the ansatz should be anti-symmetric in ${\tilde{x}}$. 4. 4.

The longitudinal dependence of the ansatz should be the same as the SH source, i.e., $e^{2i\alpha^{\omega}|{\tilde{y}}|/d}$ as it is specified for all ${\tilde{y}}$. 5. 5.

Having set the propagation constant of the ansatz along ${\tilde{y}}$ to be $\alpha^{\omega}$, the dispersion relation gives us the propagation constant along the transverse direction ${\tilde{x}}$ as $k_{{\tilde{x}},bg/m}=2\sqrt{(\alpha^{\omega})^{2}-\left(\frac{\omega}{c}\right)^{2}\varepsilon^{2\omega}_{bg/m}}$. However, we operate in the quasi-static regime, i.e., in the limit $\omega/c\rightarrow 0$, therfore $k_{{\tilde{x}},bg/m}$ simplifies to $2\alpha^{\omega}$, making the transverse propagation constant the same as that of the longitudinal one in all the regions.

Following the above considerations and assuming that $H_{z,sl}^{2\omega}$ is variable separable, we arrive at the ansatz

[TABLE]

The amplitudes $\mathcal{A}_{sl}$ and $\mathcal{B}_{sl}$ in Eq. (16) are the slowly varying amplitudes  333The ansatz (16) satisfies Eq. (10) only when the amplitudes $\mathcal{A}_{sl}$ and $\mathcal{B}_{sl}$ are slowly varying such that their spatial derivatives can be neglected. in ${\tilde{y}}$ such that the ansatz satisfies Eq. (10) and Eqs. (15).

Demanding the continuity of $H_{z,sl}^{2\omega}$ across the interfaces (see Eq. (11)) gives

[TABLE]

Implementing the transformed generalized boundary condition at the right and left interfaces (see Eq. (15)) and using Eq. (17) yields

[TABLE]

Let us now interpret the SH solution (16)-(19). First, we note that the spatial dependence of the solution (16) has a striking resemblance to an anti-symmetric mode of the slab [41, 25, 45, 24]. However, potentially unexpectedly (e.g., based on analysis such as coupled mode theory, see e.g., [25, 24, 46]), it does not necessarily correspond to the excitation of a single mode of the structure because the phase accumulation rate (i.e., the ${\tilde{y}}$ dependence) is determined by the source rather than by a propagation constant of a given (anti-symmetric) mode; in the quasi-static limit, these dimensionless propagation constants $\alpha^{2\omega}$ are given by the roots of the dispersion relations [41]444For the electrodynamic case, a root search is required [70, 71].

[TABLE]

which has solutions only for

[TABLE]

which is above the surface plasmon resonance frequency.

A remarkable aspect of conformal transformation is that the source exhibits harmonic phase variation along ${\tilde{y}}$ (see Eq. (14)), which is preserved in the solution via momentum conservation. This also enables elegant phase matching interpretations. Firstly, one can appreciate that the absence of propagation constants of the SH modes from the solution (16) is in fact expected, since the configuration we study is analogous to the scattering problem from a Fabry-Pérot etalon. Indeed, the spatial profile of the electric fields is known in each domain, and the corresponding unknown amplitudes are calculated by ensuring the satisfaction of the boundary conditions. In that procedure, the scattered fields have the same longitudinal momentum as that of the source (incident field) and the excitation is resonant only when the source momentum matches that of a mode of the etalon.

We now employ the powerful phase matching interpretation of the result. A SH mode is resonantly excited in the slab only if the longitudinal momentum component of the SH source matches the momentum of the (anti-symmetric) mode. This can be inferred by comparing Eq. (20) to $\mathcal{P}$ in Eq. (19) where we see that PM is obtained when $2\alpha^{\omega}=\alpha^{2\omega}$. Thus, the above reveals the analogy between a resonant excitation and PM for SHG. Further, we note that PM results in the enhancement of the SH fields as $H_{z,sl}^{2\omega}$ diverges for $\mathcal{P}\rightarrow 0$. Unfortunately, a comparison of the expressions of $\alpha^{\omega}$ (Eq. (13)) and $\alpha^{2\omega}$ (Eq. (20)) shows that due to material dispersion, the zeros of $\mathcal{P}$ cannot be reached without special means (see e.g., [25, 24]); in practice the values taken by $\mathcal{P}$ are even greater in the presence of absorption, as these push the poles to the complex plane of the parameter domain.

An additional condition for efficient SHG usually found in the literature is the so-called mode matching (MM) condition. In contrast to PM discussed above, which essentially involves the spatial overlap of the source and mode in the longitudinal coordinate, MM refers to the spatial overlap of the mode and source in the direction perpendicular to the propagation; maximal overlap is desired to achieve efficient conversion. For symmetric structures, the overlap integral has a “binary” interpretation - a symmetric (anti-symmetric) source can only excite the symmetric (anti-symmetric) modes. Furthermore, for the line sources arising from the surface nonlinearity studied here, the overlap simplifies drastically to the product of the amplitudes on each of the interfaces. As explained, in our case, the anti-symmetric combination of the SH line sources at the interfaces ($\Delta_{r}$ and $\Delta_{l}$) dictates that the SH ansatz (16) has to be anti-symmetric along the transverse coordinate ${\tilde{x}}$. As a consequence, the condition $-\varepsilon_{bg}^{2\omega}<\textrm{Re}(\varepsilon_{m}^{2\omega})<\varepsilon_{bg}^{2\omega}$ in Eq. (21) reveals that PM can be achieved only for frequencies above the surface plasmon resonance frequency. Since the modes occurring above the surface plasmon resonance frequency are very lossy, the SH efficiency is expected to be low.

All the above clarifies the role of phase-matching and mode-matching conditions in the slab geometry for SHG with surface sources. Apart from the above, we also note an additional feature in our solution when compared to the usual slab solution, namely, the term $(d^{2}+{\tilde{y}}^{2})$ in the denominator of the SH solution (18). As mentioned earlier, it originates from the transformation of the boundary condition (6) from the TW frame to the slab frame. Such a term is absent when one calculates the SHG of the waveguide geometry in isolation (as e.g., in [24, 25][48, Ch. 4]), meaning that it is not intrinsic to the waveguide, but rather, is a result of transforming the TW structure. The $1/(d^{2}+{\tilde{y}}^{2})$ dependence determines the spatial variation of the SH amplitudes $\mathcal{A}_{sl},~{}\mathcal{B}_{sl}$ (see Eqs. (16)-(18)) along ${\tilde{y}}$ together with the SH source $\Delta_{r}$; it suppresses SH source strength and thereby, the strength of the SH fields for large ${\tilde{y}}$.

II.3 Transformation back to the TW frame

Our next step is to transform the SH solution $H_{z,sl}^{2\omega}$ back to the TW frame. Since $H_{z,sl}^{2\omega}$ is preserved under the transformation, see Appendix B, $H_{z,TW}^{2\omega}$ can be obtained by just rewriting the slab coordinates in terms of the TW coordinates. Thus, $H_{z,TW}^{2\omega}$ is given by

[TABLE]

where

[TABLE]

Comparing the spatial dependence of $\mathcal{R}$ in Eq. (12a) to Eq. (25) reveals that $\mathcal{R}$ now gains a different spatial dependence denoted by $\tau_{x,y}$; $\tau_{x,y}$ in Eqs. (26)-(27) map every point in the TW plane ($x-y$ plane) onto the circumference of the right wire while, the circumference of the right wire (i.e., when $x^{2}+y^{2}=2ax$) is mapped onto itself, i.e., $\tau_{x}=x,~{}\tau_{y}=y$. This map is a consequence of Eq. (16) and Eq. (18) as the amplitude $\mathcal{A}_{sl}$ in the slab frame is determined by the spatial variation of $\Delta_{r}$ at ${\tilde{x}}=d$.

We now focus on the spatial variation of $H_{z,TW}^{2\omega}$ on the circumference of the right wire. Since this curve maps onto itself, the functional dependence of $\mathcal{R}(\tau_{x},\tau_{y})$ in Eq. (25) simplifies to $\mathcal{R}(x,y)$. Also the additional pre-factor in Eq. (25) simplifies to $2x$. Thus, $H_{z,TW}^{2\omega}$ is given by

[TABLE]

and upon using Eq. (12a), the above can be rewritten as

[TABLE]

This analytic solution is the main result of the paper.

In order to validate our analytical solution, we compare it with a numerical solution, obtained using a commercially-available finite element method software package COMSOL Multiphysics (see Appendix E for further details). Figure 3(a) shows an excellent agreement between the analytic and numeric solutions of $H_{z,TW}^{2\omega}$ on the circumference of the right wire close to the touching point ($\theta=180^{\circ}$). In Figure 3(b), we see that this agreement extends even to angles as large as several tens of degrees away from the touching point. This validates our analytical solution.

Similar to the linear response, the SH response exhibits interesting physical phenomena like wavelength compression, slow light, field enhancement and energy accumulation close to the touching point [27], the latter being inhibited by the absorption. In what follows, we interpret the SH solution of the TW, and highlight how the terms that appeared in the solution for the slab manifest themselves in the solution for the TW geometry and how they affect the above mentioned phenomena.

First, the momentum mismatch between the mode and the SH source is preserved under the transformation, showing that (unlike previous claims in the literature (e.g., [26])), the phase-matching condition is relevant even for this sub-wavelength structure 555For the set of parameters considered in Figure 3, $|\mathcal{P}|\approx 1$. This implies that the SHG indeed occurs under phase-mismatched condition.. Indeed, the link between the TW geometry and the slab geometry revealed by the coordinate transformation shows that if PM is important for the slab geometry, then, it must be important also for the particle (i.e., TW) geometry. Such an effect is possible due to the (non-uniform) wavelength compression induced by the (touching) wires which makes the optical length along the wires effectively infinite [5, 27]. We note that such a result is nontrivial to reproduce by a multipolar expansion, despite the deep subwavelength nature of the TW geometry [50]. Second, MM manifests through the symmetry of the SH sources for the identical TW system. The $J_{z,r/l}$ is anti-symmetric in nature (see Figure 2(c)), thus generating an anti-symmetric $H_{z,TW}^{2\omega}$ field.

Apart from the phase-matching and mode-matching conditions, the SH solution also exhibits an additional spatial dependence which occurs as the prefactor in $\mathcal{C}$ (see Eq. (25)). We coin this term as the geometric factor. The geometric factor stems from transforming $1/(d^{2}+{\tilde{y}}^{2})$ in Eq. (18) to the TW frame. Recall that $1/(d^{2}+{\tilde{y}}^{2})$ was obtained from transforming generalized boundary condition from TW to slab frame, however, its signature still remains in the final solution as we transform only the SH magnetic field back from the slab to the TW frame leaving out the transformation of generalized boundary condition (15).

As mentioned, the geometric factor on the circumference of the TW simplifies to $2x$, i.e., it attains small values close to the touching point ($x\approx 0$). Since the magnetic field solution (see Eq. (28)) is given by the product of the SH source and the geometric factor, the SH source (and thereby, the SH response) is further suppressed by the geometric factor close to the touching point. Therefore, the SH solution decays faster to the touching point compared with the linear solution.

The origin of the geometric factor can be understood as follows. In the usual Green’s function approach, a distributed source should be integrated over its spatial coordinates along with the Green’s tensor corresponding to the structure, namely, $\int\bar{\bar{G}}_{TW}(\mathbf{r},\mathbf{r}^{\prime})\mathbf{P}^{2\omega}(\mathbf{r}^{\prime})d\mathbf{r}^{\prime}$. This integral contains all the information about the structure (in particular, the modal structure) and the spatial distribution of the SH source, and the convolution naturally yields PM and MM. However, prior to the evaluation of this integral, the final spatial distribution of the SH solution is not known. The conformal transformation reveals that the SH solution (22)-(27) has a remarkably simple form, namely, it is linearly proportional to the SH source, phase-matching and mode-matching factors. Additionally, we encounter the geometric factor which can be thought of as a simplification of the Green’s function integral with the distributed source corresponding to the TW geometry. We note that such a simplification/factorization of SH solution is realizable due the power of conformal transformation by transforming the seemingly difficult SH problem to the much simpler system, i.e., the slab frame and thereby, enabling us to relate the local SH solution in terms of the the local source. Yet, this simplicity pertains only exactly at the metal-dielectric interface of the TW, and is far more complicated elsewhere.

III Discussion and Outlook

Using the technique of conformal transformations, we have related the SHG in waveguides to the SHG in a far more complicated geometry - the touching wire dimer. Our ansatz approach is simpler compared to the treatment of the linear TW problem with CT, as we avoid the need to transform to momentum space and perform contour integrations. The transformation allows us to unfold the hidden symmetries, separating variables such that the roles of PM and MM become obvious also in the TW frame, making the interpretation of the final solution quite straightforward.

In particular, this showed the equivalence of phase-matching in the slab geometry to tuning to the localized surface plasmon resonance at the SH for the TW structure. However, our analysis shows rather surprisingly, that PM and MM and source strength are accompanied by an additional factor. This additional factor was not identified before and we refer to it as a geometric factor. Indeed, its presence implies that different original geometries might give rise to different sources, hence different geometric factors. This shows that the solution for the particle geometry exhibits richer physics compared to the slab. In particular, the revealed complexity demonstrates the limitations of approaches based on just evaluating the strength of the source $P^{2\omega}$, or just on PM 666In [72], the SHG efficiency was optimized via an integral over the longitudinal coordinate involving the surface charges. While in [72] it was referred to as mode matching, adopting the link between the particle and waveguide geometries revealed in this manuscript, we argue that it in fact represents phase matching instead..

This work introduces new tools and insights into nonlinear optical wave mixing for nanoscale structures (and on the treatment of distributed sources on the nanoscale [52, 53, 54]). It differs from many studies by providing a unique analytic near-field solution rather than on a qualitative description of the (experimentally-accessible) far-field pattern, and by going beyond qualitative symmetry breaking arguments [55, 56, 57, 58].

The solution approach adopted in this manuscript can serve as a means to study and optimize the near-field and the spectral response of the touching wire dimer. It can also be used to calculate the SH field distributions analytically for various other singular nanoparticle structures like asymmetric touching wires, crescent structures, wedge structures, circular protrusion from a planar interface, etc. [8]. Additionally, this approach can also be extended to calculate the SH response from non-singular structures such as non-touching wires [59], blunt crescent structures [60].

The solution approach adopted in this manuscript is suitable also for other elements of the surface polarization tensor as well as for (non-local) bulk polarizations (by mapping them to a surface polarization [30, 31, 32, 33]).

Other problems that can be solved with the same approach are 3D structures [7, 61], SHG at other spectral regimes, e.g., in the THz regime [62], optical rectification, and more complex nonlinear wave interactions, such as 3 wave mixing [63, 64, 65], phase-conjugation [66] or even additional effects such as the role of non-locality [14, 15] on SH and many more.

Acknowledgements.

The authors would like to thank A. Niv and A. Isha’aya for many useful discussions. KNR and YS were partially supported by Israel Science Foundation (ISF) grant (899/16). KNR and AIFD would like to acknowledge STSM Grant from the COST Action MP1403. 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.

Appendix A Extracting the magnetic fields in the quasi-static regime

The linear electric fields were obtained by setting $\nabla\times\mathbf{E}^{\omega}=0$ [5], thus, the linear magnetic field was set to zero. For the quasi-static structures, the electric response dominates. However, the magnetic response even though negligible, is not truly zero. Here, we demonstrate on how to extract the magnetic field from the quasi-static structures.

In general, this magnetic field can be evaluated in two different ways. One can express the excitation sources in terms of the magnetic sources and solve for the magnetic response. For example, a plane wave illumination was modelled using magnetic line currents and the magnetic response for various quasi-static structures was obtained [67]. Similarly, in this manuscript, we formulate the nonlinear polarization in terms of a surface magnetic current to obtain the SH magnetic response of the TW system. Alternatively, despite setting $\nabla\times\mathbf{E}=0$, one can employ Ampere’s law to extract the magnetic field from the evaluated electric fields. In what follows, we extract the linear magnetic field $H_{z,TW}^{\omega}$ of the TW system under plane wave illumination using the latter.

Using Ampere’s law explicitly for our case yields

[TABLE]

with $\varepsilon_{bg/m}^{\omega}$ considered in the appropriate domains. $E_{x,TW}^{\omega}$ and $E_{y,TW}^{\omega}$ in Eqs. (29a)-(29b) correspond to the linear electric field solutions of the identical TW system in different domains as obtained by Lei et al. [27]777Note that Lei et al. obtained the linear response for $\varepsilon_{bg}^{\omega}=1$. However, one can extend this linear response to arbitrary $\varepsilon_{bg}^{\omega}$ by substituting $\varepsilon_{m}^{\omega}\rightarrow\varepsilon_{m}^{\omega}/\varepsilon_{bg}^{\omega}$. Upon integrating Eqs. (29a)-(29b) and determining the constants of integration we arrive at

[TABLE]

with $\mathcal{N}=-2\pi i\omega\varepsilon_{0}aE^{\omega}_{0x}$.

Alternatively, one can also arrive at magnetic field $H_{z,TW}^{\omega}$ by obtaining the magnetic field in the slab geometry can then transforming to the TW frame. The out-of-plane component of the magnetic field is preserved under the 2D conformal transformation (see Appendix B) [67]. We have verified that both the methods are consistent.

Figure 4 shows the comparison between the analytical and numerical solutions of $H_{z,TW}^{\omega}$ on the circumference of the right wire of TW system. The numerical solution was obtained using COMSOL Multiphysics by solving for electrodynamic response from TW under plane wave illumination. We find that the analytical solution is in excellent agreement with the numerical solution.

Appendix B Transformation of the magnetic field

Here, we derive the transformation rules of the magnetic field under 2D conformal transformation. It is important to note that the transformation considered here is completely in-plane, i.e., 2D, thus, the out-of-plane coordinate $z$ remains unaffected. Following the TO formulas for a vector quantity [1], the magnetic fields in the TW frame $\mathbf{H}_{TW}$ and the slab frame $\mathbf{H}_{sl}$ are related as

[TABLE]

where $\Lambda_{ij}$ is the Jacobian of the transformation. Equating the $z$-components in Eq. (31) we arrive at

[TABLE]

Thus, we have obtained that the $z$-component of the magnetic field is invariant under the transformation [67]. This conclusion holds true for both FF and SH magnetic fields.

It also follows from Eqs. (31)-(32) that any vector pointing normal to the $x-y$ plane or ${\tilde{x}}-{\tilde{y}}$ plane, for example, SH source $J_{z,r/l}$, is preserved under 2D conformal transformation.

Appendix C Amplitudes of SH source - $\mathcal{R}$ and $\mathcal{L}$

The amplitudes of the SH sources $\mathcal{R}$ and $\mathcal{L}$ for the TW as defined in Eqs. (12) can be evaluated as follows. The FF response $E_{\perp}^{\omega}$ at the metal-dielectric interface on the right and left cylinders(on the metal side) is given by [27]

[TABLE]

respectively. Other variables in Eqs. (33) are defined as

[TABLE]

The variable $\theta_{r}$ ($\theta_{l}$) in Eq. (34) (Eq. (35)) corresponds to the angular co-ordinate on the circumference of the right (left) cylinder of radius $a$ centered at $x=a$ ($x=-a$). The SH source (3) on the right (left) cylinder in terms of $E_{\perp,r}^{\omega}$ ($E_{\perp,l}^{\omega}$) is given by

[TABLE]

respectively. In writing Eqs. (36) we have used the fact that the tangential derivative $\partial_{\parallel}$ on the right (left) cylinder is given by $\frac{1}{a}\partial_{\theta_{r}}$ ($\frac{1}{a}\partial_{\theta_{l}}$). Rewriting $J_{z,r}$ and $J_{z,l}$ (in Eqs. (36)) in ($x,y$) and invoking the decomposition of $J_{z,r/l}$ from Eq. (12) gives us $\mathcal{R}$ and $\mathcal{L}$ as

[TABLE]

respectively.

Appendix D Transformation of the generalized boundary condition to the slab frame

Here, we compute the transformation of the generalized boundary condition (7) from TW to slab frame. In what follows, we first compute the transformation of the left-hand-side of Eq. (7) and then followed by its right-hand-side. The magnetic field $H_{z,TW}^{2\omega}$ can be replaced by $H_{z,sl}^{2\omega}$ (see Appendix B). Replacing the spatial coordinates and the derivatives of TW frame by that of the slab ones shows that the left-hand-side of Eq. (7) transformations as

[TABLE]

where

[TABLE]

We now compute the transformation of the right-hand-side of Eq. (7). The magnetic current $J_{z,r/l}$ is a vector normal to the $x-y$ plane, thus, it is preserved under the transformation. Invoking the decomposition of $J_{z,r/l}$ from Eqs. (12) and rewriting it in terms of the slab variables gives

[TABLE]

In writing the above expressions, we have relabeled $\mathcal{R}$ and $\mathcal{L}$ in slab geometry by

[TABLE]

respectively. Since the conformal transformation (8) preserves the symmetry relation between $\mathcal{R}$ and $\mathcal{L}$, we have $\Delta_{r}({\tilde{x}}=d,{\tilde{y}})=-\Delta_{l}({\tilde{x}}=-d,{\tilde{y}})$. It can seen from Eq. (39a) (Eq. (39b)) that the surface source $J_{z,r}$ ($J_{z,l}$) on the circumference of the right (left) wire transforms as the source placed at right (left) interface of the slab geometry.

Thus, by combining Eqs. (38)-(39), we arrive at the transformed generalized boundary condition given by

[TABLE]

Appendix E Numerical Simulations

The numerical solutions are obtained using the commercially available finite element method software package COMSOL Multiphysics 3.5a. We have used harmonic propagation analysis in the radio frequency module to obtain the SH electrodynamic response of the TW system. The FF analytical expressions were used to compute the SH source $J_{z,r/l}$ (3) and were given as an input to the COMSOL 3.5a solver. Specifically, SH source $J_{z,r/l}$ was inserted as a source in the ‘Boundary setting - Equation system’ node. The region close to the touching point was resolved with mesh sides below $10^{-4}$ nm to sample the rapidly oscillating $J_{z,r/l}$ accurately and to tackle the geometric singularity. A perfectly matched layer was used to suppress the reflections from the boundaries of the simulation domain. The convergence of the numerical solution as a function of mesh size and size of the simulation domain was verified. Despite the availability of the latest versions (COMSOL 5.3a), we were forced to resort to the old version COMSOL 3.5a due to its lower sensitivity to the ill-conditioning associated with the extremely fine mesh near the touching point. A detailed account concerning the sensitivity of COMSOL 5.3a to the small mesh and evaluation of the linear and SH response is presented in Ref. [48, Appendix D].