A non-singular, field-only surface integral method for interactions between electric and magnetic dipoles and nano-structures

In the frequency domain, the electric and magnetic field due to an electric dipole are, respectively,

where $\boldsymbol{p}$ is the electric dipole momentum vector (with dimension Coulomb meter), $\boldsymbol{x}$ is the observation location and $\boldsymbol{x}_{d}$ the location of the dipole (presumed to be fixed in space, to avoid Doppler effects), $k=\omega/c$ is the wavenumber with $\omega$ the angular frequency and $c$ the speed of light, $\epsilon_{0}$ and $\epsilon_{r}$ are the permittivity of free space and relative permittivity of the medium respectively.

In the domain that excludes the dipole, thus the domain outside the scatterer $S$ and a tiny spherical surface $S_{d}$ that encloses the dipole as shown in Fig. 1, the Maxwell equations in the frequency domain reduce to the following two equations for the electric field $\boldsymbol{E}$

Thus the electric field satisfies a vector Helmholtz equation and is simultaneously divergence free. To obtain the interactions between one or more dipoles and scatterers in a homogeneous medium, it is necessary to efficiently and accurately solve Eq. (2).

It has been demonstrated that Eq. (2a) can be solved effectively and accurately by the recently developed non-singular field only surface integral method [25, 26]. To obtain the electric field on the boundary, each Cartesian component of the $\boldsymbol{E}$ field in Eq. (2a), which obeys the scalar Helmholtz equation (an elliptic partial differential equation), is represented by the non-singular boundary integral equation [29] as (see Appendix A for derivation)

where $\boldsymbol{x}_{0}$ is an observation point and $\boldsymbol{x}$ is an integration point, as shown in Fig. 1. $G_{k}(\boldsymbol{x},\boldsymbol{x}_{0})=\exp(\mathrm{i}k|\boldsymbol{x}-\boldsymbol{x}_{0}|)/|\boldsymbol{x}-\boldsymbol{x}_{0}|$ is the Green’s function for the Helmholtz equation and $G_{0}=1/|\boldsymbol{x}-\boldsymbol{x}_{0}|$ the Laplace equivalent. $H_{k}=\partial G_{k}/\partial n$ is the normal derivative of the Green’s function as in the standard boundary element framework with $\partial/\partial n=\boldsymbol{n}\cdot\nabla$ where $\boldsymbol{n}$ is the unit normal vector pointing into the object. The classical boundary integral equation is made non-singular analytically by subtracting a function $\boldsymbol{F}(\boldsymbol{x})$:

which satisfies the Laplace equation as $\nabla^{2}\boldsymbol{F}(\boldsymbol{x})=\boldsymbol{0}$. Note that both terms in each bracket of Eq. (3) cancel each other out when $\boldsymbol{x}$ approaches $\boldsymbol{x}_{0}$ and effectively eliminate the singularities analytically which arise due to the Green’s function.

The term $\boldsymbol{I}^{\mathrm{Ed}}$ represents the integral on $S_{d}$ which is a tiny spherical surface with vanishing radius $r_{d}$ that encloses the dipole:

In Appendix A, it is shown that the term $\boldsymbol{I}^{\mathrm{Ed}}$ can be written as:

in which $\boldsymbol{x}_{0d}=\boldsymbol{x}_{0}-\boldsymbol{x}_{d}$, $r_{0d}=|\boldsymbol{x}_{0d}|$ and $G_{k}(\boldsymbol{x}_{d},\boldsymbol{x}_{0})=G_{k}(\boldsymbol{x}_{0},\boldsymbol{x}_{d})$ are used.

The above boundary integral equation (3) gives a relationship between the Cartesian components of the electric field ($E_{x}$, $E_{y}$, $E_{z}$) and their normal derivatives ($\partial E_{x}/\partial n$, $\partial E_{y}/\partial n$, $\partial E_{z}/\partial n$) on the surface $S$. However, the divergence free condition in the domain, $\nabla\cdot\boldsymbol{E}=0$, also needs to be satisfied simultaneously. Sun et al. [25] demonstrated that when the divergence free condition is satisfied everywhere on the surface of the object, it is automatically satisfied everywhere in the domain. It has also been shown in Sun et al. [25] that the divergence free condition $\nabla\cdot\boldsymbol{E}=0$ is equal to

Here $\kappa$ is the curvature on the surface, $E_{n}\equiv\boldsymbol{E}\cdot\boldsymbol{n}$ is the normal component of the electric field on the surface, and $E_{t1}\equiv\boldsymbol{E}\cdot\boldsymbol{t}_{1}$, $E_{t2}\equiv\boldsymbol{E}\cdot\boldsymbol{t}_{2}$ are the tangential components with $\boldsymbol{t}_{1}$ and $\boldsymbol{t}_{2}$ the two orthogonal surface unit tangential vectors and $\partial/\partial t_{1}=\boldsymbol{t}_{1}\cdot\nabla$ (and similar for $\partial/\partial t_{2}$).

Since the boundary conditions are given in terms of tangential components of the electric field, it turns out to be easier to use the normal and tangential components as unknowns in our framework. As such, we rewrite the electric field and its normal derivative as

Thus for example $E_{x}$ can be expressed as $E_{x}=E_{n}(\boldsymbol{n}\cdot\boldsymbol{e}_{x})+E_{t1}(\boldsymbol{t}_{1}\cdot\boldsymbol{e}_{x})+E_{t2}(\boldsymbol{t}_{2}\cdot\boldsymbol{e}_{x})=E_{n}n_{x}+E_{t1}t_{1x}+E_{t2}t_{2x}$.

To numerically solve Eq. (3) for the interaction between a dynamic electric dipole and a PEC object, the object surface is discretised with $N$ nodes. Then, the left hand side of Eq. (3) results in a matrix equation relating all $x$-components of the electric field as: ${\cal{H}}E_{x}$ and its normal derivatives on the right hand side of Eq. (3) as: ${\cal{G}}\partial E_{x}/\partial n$. Here, ${\cal{G}}$ and ${\cal{H}}$ are $N\times N$ matrices and $E_{x}$ and $\partial E_{x}/\partial n$ are now column vectors of size $N$. A similar procedure applies for the $y$ and $z$-components. To introduce boundary condition in Eq. (7) into the linear system, the Cartesian components of the electric field and its normal derivatives can be reconstructed from the normal and tangential components with Eq. 8. Also, for a PEC object, the surface tangential components of the electric field are zeros. As such, we get the following $3N\times 3N$ matrix system as

from which the electric field on the PEC surface can be obtained. The interested reader can find more details on the implementation in Sun et al. [25].

For a dielectric object, we need to calculate the field inside the object $\boldsymbol{E}^{\mathrm{in}}$ as well as the external electric field on both sides of the surface $\boldsymbol{E}^{\mathrm{out}}$. The internal domain is characterised by $k_{in}$, $\epsilon_{\mathrm{in}}$, $\mu_{\mathrm{in}}$ while the outer domain by $k_{\mathrm{out}}$, $\epsilon_{\mathrm{out}}$ and $\mu_{\mathrm{out}}$. The field only surface integral equation for the external domain s:

For the internal domain we find:

In Eqs. (10) and (11), $\boldsymbol{F}^{\mathrm{out}}(\boldsymbol{x})=\boldsymbol{E}^{\mathrm{out}}(\boldsymbol{x}_{0})+\boldsymbol{n}(\boldsymbol{x}_{0})\cdot(\boldsymbol{x}-\boldsymbol{x}_{0})\frac{\partial\boldsymbol{E}^{\mathrm{out}}(\boldsymbol{x}_{0})}{\partial n}$ and similar for $\boldsymbol{F}^{in}$ as $\boldsymbol{F}^{\mathrm{in}}(\boldsymbol{x})=\boldsymbol{E}^{\mathrm{in}}(\boldsymbol{x}_{0})+\boldsymbol{n}(\boldsymbol{x}_{0})\cdot(\boldsymbol{x}-\boldsymbol{x}_{0})\frac{\partial\boldsymbol{E}^{\mathrm{in}}(\boldsymbol{x}_{0})}{\partial n}$. Note the difference of a term with $4\pi\boldsymbol{E}^{\mathrm{out}}$ between Eqs. (10) and (11) due to the absence of integrals at infinity for the internal problem in Eq. (11) (see Appendix A).

The boundary conditions for the electric field, expressed in terms of the normal and tangential components, are:

with $\epsilon_{\mathrm{oi}}=\epsilon_{\mathrm{out}}/\epsilon_{\mathrm{in}}$ the ratio of permittivities between the two media. Meanwhile, the divergence condition of Eq. (7) is applied to $\boldsymbol{E}^{\mathrm{out}}$ as $\boldsymbol{n}\cdot\frac{\partial\boldsymbol{E}^{\mathrm{out}}}{\partial n}-\kappa E_{n}^{\mathrm{out}}+\frac{\partial E_{t_{1}}^{\mathrm{out}}}{\partial t_{1}}+\frac{\partial E_{t2}^{\mathrm{out}}}{\partial t_{2}}=0$ and to $\boldsymbol{E}^{\mathrm{in}}$ as $\boldsymbol{n}\cdot\frac{\partial\boldsymbol{E}^{\mathrm{in}}}{\partial n}-\kappa E_{n}^{\mathrm{in}}+\frac{\partial E_{t_{1}}^{\mathrm{in}}}{\partial t_{1}}+\frac{\partial E_{t2}^{\min}}{\partial t_{2}}=0$. In order to get rid of the tangential derivative terms these two equations are subtracted while the boundary conditions of Eq. (12b) and (12c) are used. As such, we can express $E_{n}^{\mathrm{in}}$ in terms of $E_{n}^{\mathrm{out}}$ with Eq. (12a):

This equation relates the normal component of the normal derivatives, $\boldsymbol{n}\cdot(\partial\boldsymbol{E}/\partial n)$, at both sides of the boundary of the object. We still need another equation to relate both tangential components of the normal derivatives: $\boldsymbol{t}\cdot(\partial\boldsymbol{E}/\partial n)$. This can be done by using the tangential continuity of the magnetic field on the dielectric surface: $H_{t1}^{\mathrm{in}}=H_{t1}^{\mathrm{out}}$ and $H_{t2}^{\mathrm{in}}=H_{t2}^{\mathrm{out}}$. Using the convention $\boldsymbol{n}=\boldsymbol{t}_{1}\times\boldsymbol{t}_{2}$, $\boldsymbol{t}_{1}\cdot\boldsymbol{t}_{2}=0$, we express $H_{t1}$ and $H_{t2}$ in terms of the electric field on the surface as

It is worth noting the inversion of $t_{1}$ and $t_{2}$ and the appearance of a minus sign in the $H_{t2}$ term in the above equation.

For simplicity, we assume that $\mu_{\mathrm{in}}=\mu_{\mathrm{out}}$ (if this is not the case, see Sun et al. [26]), and then the tangential continuity of $\boldsymbol{H}$ across the dielectric surface is identical to

The terms with $\boldsymbol{n}\cdot\partial\boldsymbol{E}/\partial t_{1}$ and $\boldsymbol{n}\cdot\partial\boldsymbol{E}/\partial t_{2}$ can be rewritten as¹¹1This can easiest be seen by writing the vector $\boldsymbol{E}$ into normal and tangential components as: $\partial\boldsymbol{E}/\partial t_{1}=\frac{\partial}{\partial t_{1}}[E_{n}\boldsymbol{n}+E_{t1}\boldsymbol{t}_{1}+E_{t2}\boldsymbol{t}_{2}]=\boldsymbol{n}\frac{\partial E_{n}}{\partial t_{1}}+E_{n}\frac{\partial\boldsymbol{n}}{\partial t_{1}}+\boldsymbol{t}_{1}\frac{\partial E_{t1}}{\partial t_{1}}+E_{t1}\frac{\partial\boldsymbol{t}_{1}}{\partial t_{1}}+\boldsymbol{t}_{2}\frac{\partial E_{t2}}{\partial t_{1}}+E_{t2}\frac{\partial\boldsymbol{t}_{2}}{\partial t_{1}}$. Using $\partial\boldsymbol{n}/\partial t_{1}=-\kappa_{1}\boldsymbol{t}_{1}$, $\partial\boldsymbol{t}_{1}/\partial t_{1}=\kappa_{1}\boldsymbol{n}$ and $\partial\boldsymbol{t}_{2}/\partial t_{1}=\boldsymbol{0}$, we get $\boldsymbol{n}\cdot\frac{\partial\boldsymbol{E}}{\partial t_{1}}=\frac{\partial E_{n}}{\partial t_{1}}+\kappa_{1}E_{t1}$ and similar for $\boldsymbol{n}\cdot\frac{\partial\boldsymbol{E}}{\partial t_{2}}=\frac{\partial E_{n}}{\partial t_{2}}+\kappa_{2}E_{t2}$. $\boldsymbol{n}\cdot\partial\boldsymbol{E}/\partial t_{1}=\partial E_{n}/\partial t_{1}+\kappa_{1}E_{t1}$ and $\boldsymbol{n}\cdot\partial\boldsymbol{E}/\partial t_{2}=\partial E_{n}/\partial t_{2}+\kappa_{2}E_{t2}$, with $\kappa_{1}$ the curvature in the $\boldsymbol{t}_{1}$ direction and $\kappa_{2}$ the curvature in the $\boldsymbol{t}_{2}$ direction. Putting this into Eq. (15) where the terms with $\kappa_{1}E_{t1}$ and $\kappa_{2}E_{t2}$ cancel out due to Eqs. (12b) and (12c), and further using Eq. (12a) to express $E_{n}^{\mathrm{in}}$ in terms of $E_{n}^{\mathrm{out}}$, we have $\frac{\partial E_{n}^{\mathrm{out}}}{\partial t_{1}}+\kappa_{1}E_{t1}^{\mathrm{out}}-\frac{\partial E_{n}^{\mathrm{in}}}{\partial t_{1}}-\kappa_{1}E_{t1}^{\mathrm{in}}=\frac{\partial E_{n}^{\mathrm{out}}}{\partial t_{1}}(1-\epsilon_{\mathrm{oi}})$. As such,

With Eq. (12), we now can express all the components of $\boldsymbol{E}^{\mathrm{in}}$ in terms of $\boldsymbol{E}^{\mathrm{out}}$. All components of the normal derivatives of the internal field can also be expressed in terms of outer field variables with Eq. (13) and Eq. (16). The Cartesian $x$, $y$ and $z$ components can then be obtained with Eq. (8).

To numerically solve Eqs. (10) and (11) for the interaction between a dynamic electric dipole and a dielectric object, the object surface is discretised with $N$ nodes. Meanwhile, the $x$, $y$, $z$ components of $\boldsymbol{E}^{\mathrm{in}}$ and $\boldsymbol{E}^{\mathrm{out}}$ and their normal derivatives are expressed in terms of $E_{n}^{\mathrm{out}}$, $E_{t1}^{\mathrm{out}}$, $E_{t2}^{\mathrm{out}}$, $\boldsymbol{n}\cdot(\partial\boldsymbol{E}^{\mathrm{out}}/\partial n)$, $\boldsymbol{t}_{1}\cdot(\partial\boldsymbol{E}^{\mathrm{out}}/\partial n)$ and $\boldsymbol{t}_{2}\cdot(\partial\boldsymbol{E}^{\mathrm{out}}/\partial n)$, which results in a $6N\times 6N$ matrix system as

In the above matrix system, the first line corresponds to ${\cal{H}}_{\mathrm{out}}E_{x}^{\mathrm{out}}={\cal{G}}_{\mathrm{out}}\partial E_{x}^{\mathrm{out}}/\partial n$, the matrix equivalent of Eq. (10) for the $x$-component. Line 2 and 3 correspond to the $y$ and $z$-component, respectively. Lines 4, 5 and 6 correspond to the matrix equivalent of Eq. (11) for the $x$, $y$ and $z$-component of $\boldsymbol{E}^{\mathrm{in}}$, respectively. The implementation is straightforward, expect perhaps for the tangential derivative matrices $\partial/\partial t_{1}$ and $\partial/\partial t_{2}$ (for more details see Sun et al. [26]). Once the matrix system is solved we can get back the Cartesian components with Eq. (8).

As we now have the electric field and its normal derivative on the surface for both the internal and external domain, we can get the electric field anywhere in the domain by postprocessing (see Sun et al. [30]).

The numerical framework was tested with an electric dipole placed at the centre of a dielectric sphere for which an analytical solution exists. The radius is indicated with $a$, and the physical parameters are chosen to be $k_{\mathrm{out}}a=1$, $k_{\mathrm{in}}a=2$, $\epsilon_{\mathrm{in}}/\epsilon_{\mathrm{out}}=4$ and $\mu_{\mathrm{in}}=\mu_{\mathrm{out}}=1$. The number of nodes is 642 and 320 quadratic triangular elements were used to represent the sphere surface in the simulation with the dipole moment being $\boldsymbol{p}=(0,0,1)$. The analytical solution is described in Appendix B. We checked the numerical solutions of $\boldsymbol{E}$ on both sides of the surface and they show excellent agreement with the theory (not shown here). We also chose a circle with radius $1.2a$ at the plane $y=0$ to confirm that the electric field is correct in the domain close to the sphere. 72 sample locations were seeded along that circle with equal angular interval of polar angle $\theta$ where sample 1 corresponds to $\theta=0$. The results are shown in Fig. 2b, where the numerical calculations of the real parts of $E_{x}^{\mathrm{out}}$ (black symbols) and $E_{z}^{\mathrm{out}}$ (red symbols) are plotted as a function of the polar angle $\theta$ (as indicated in Fig. 2a) and compared to their theoretical counterparts (lines). The theoretical and numerical solutions are virtually identical.In Fig. 2c, we investigated the convergence performance of our method with respect to the example demonstrated in Fig. 2b. As we can see, as the number of surface nodes increases, the average relative error between the analytical solution decreases and the computational time²²2All the cases in Fig. 2c were performed on a MacBook Pro with 2.6 GHz 6-Core Intel Core i7 processor and 32 GB 2400 MHz DDR4 Memory without parallel computation. increases. The relative error in that figure is defined as

where the subscript ‘ana’ stands for the analytical solution and ‘num’ for the numerical solution.

Some other tests to examine the correctness and accuracy of the constructed non-singular field only boundary element framework were also performed: if the surface $S$ in Eq. (3) is very small or if this surface is situated far away from both the dipole and the observation point, then the electric field goes back to the field of the undisturbed dipole as it should, because both integrals in Eq. (3) will disappear for these cases.

These tests give us the confidence that the constructed framework is correct and robust. In the next sections we will investigate some more physically interesting examples.

An optical antenna can be designed, synthesised and used to efficiently convert the freely propagating optical radiation to localised energy to enhance the interaction between light and mesoscopic structures [31]. In analogy to the radio antenna design, we investigate the interaction between a dynamic electric dipole and a parabolic nano-antenna transmitter device, as shown in Fig. 3. The shape and location of the nano parabolic transmitter are given by $(x,y,z)=(a\sin{\theta}\cos{\varphi},\,a\sin{\theta}\sin{\varphi},\,0.1a(\cos{\theta}-1)+0.3\sin^{2}{\theta})$ with $\theta$ the polar angle measured from the $z$-axis, $\varphi$ the azimuthal angle and $a$ the characteristic size of the nano-antenna. The dynamic dipole is positioned at the focal point of this antenna at $(0,\,0,\,1.35a)$ which is along the symmetry axis of the parabolic nano-antenna ($z$-axis). In this example, we chose the characteristic non-dimensional value $ka=20$. For such a $ka$, when the light wavelength is $\lambda=550\,\text{nm}$, the corresponding characteristic size of the parabolic nano-antenna would be $a=1.75\,\mu\text{m}$. Also, we assume that this antenna is well coated to be a perfect mirror (PEC) when the surrounding medium is air ($n_{\text{air}}=1$), and use the theory and numerical method demonstrated in Sec. 2.1 to perform the calculations. For this test, 4002 nodes and 2000 quadratic triangular surface elements were used on the antenna surface.

In Fig. 3(b) and (d), the distribution of the electric field magnitude in the area surrounding the nano-antenna when a linearly polarised dipole interacts with it are shown. When the direction of dipole moment, $\boldsymbol{p}$, is along the symmetry axis of the antenna, the conversion of the energy of the free radiation of the dipole is very low as shown in Fig. 3(d). On the other hand, when $\boldsymbol{p}$ is perpendicular to the symmetry axis of the parabolic nano-antenna, a strong localised energy flow is formed as shown in Fig. 3(b). For the circularly polarised dipole examples of Fig. 3(a) and (c), other than the obvious conversion of energy, note that an interesting phenomenon is the changing direction of the energy flow depending on the direction of the circular polarisation. Such a change of the direction can illustrate information of the spin angular momentum of the dipole.

Duan et al.[32] showed that classical electromagnetic calculations are capable of capturing the essential physics of a nanoplasmonic system even down to nanometer scale sized gaps and bridges. Near field optical probes, such as metal tips, are one of the key components for near field microscopes in nano-optics [4]. The sensitivity and resolution of such a device is strongly dependent on the field distribution around the tip of the probe. We simulate the interaction between a dynamic electric dipole with two Au probes. The dipole is located at the origin of the reference coordinate system which radiates light with wavelength $\lambda=800\,\text{nm}$. The shape of the Au probe is described by: $x=a\sin{\theta}\cos{\varphi}[1+2\cos^{2}{(\theta/2)}]/10$, $y=a\sin{\theta}\sin{\varphi}[1+2\cos^{2}{(\theta/2)}]/10$ and $z=a\cos{\theta}$ where $\theta$ is the polar angle measured from the $z$-axis, $\varphi$ is the azimuthal angle and $a$ is the characteristic size of the probe. As a result, the tip of the Au probe points up along the $z$-axis. In this particular example, we use two such Au probes with length of 2 $\mu$m ($a=1\,\mu\text{m}$). We rotate one probe $90^{\text{o}}$ around the $y$-axis clockwise and the other counter-clockwise to make their symmetry axis along the $x$-axis with their tips pointing to each other. We then position the probes along the $x$-axis with the separation between the tips as $100\,\text{nm}$ which implies the distance between the dipole to the tip of one probe to be $50\,\text{nm}$. With such a gap size, classical electromagnetic phenomena are dominant, whereas quantum effects [33] are deemed not significant.

We consider two scenarios: when the electric dipole moment is along the axis of the tips, $\boldsymbol{p}=(1,0,0)$ and when it is perpendicular to the axis of the tips, $\boldsymbol{p}=(0,0,1)$. In both scenarios, the surrounding medium is set to be air with $n_{\text{air}}=1.00$. When $\lambda=800\,\text{nm}$, the complex refractive index of Au is $n_{\text{Au}}=0.19+\mathrm{i}4.71$. As such, the coupling effects between the internal and external electromagnetic fields of the Au probes were solved using the framework described in Sec. 2.2.

In Figs. 4a and Fig. 4d, the distribution of the instantaneous induced surface charge (which is equivalent to the normal component of the surface electric field: $Re(\boldsymbol{E}\cdot\boldsymbol{n})$) is presented. When the electric dipole moment is along the axis of the tips, the electric polarisation of the radiation light from that dipole is parallel to the tip axis which induces axis-symmetric surface charge density with the highest amplitude at the ends of the tips, as demonstrated in Fig. 4a. However, when the electric dipole moment is perpendicular to the axis of the tips, the electric polarisation of the radiation light from that dipole is also perpendicular to the tip axis. This leads to opposite surface charges on the diametrically opposed surface points on the probes with neutral foremost tips, as demonstrated in Fig. 4d. The difference is also reflected by the energy propagation within the probes. We compared the amplitude of the Poynting vector $\boldsymbol{S}=1/2\text{Re}[\boldsymbol{E}\times(\boldsymbol{H})^{*}]$ (with ^∗ indicating the conjugate) for those two scenarios in Figs. 4b,c and Figs. 4e,f, respectively, where $\boldsymbol{H}$ was solved by the framework illustrated in Sec. 4. We can see that the electromagnetic energy penetrates deeper into the Au probes when the electric dipole moment is along the tip axis than when it is perpendicular to that axis.