Electromagnetic energy in multilayered spherical particles

TL;DR

This paper derives exact analytical formulas for electromagnetic energy distribution in multilayered spherical particles, applicable to both dielectric and magnetic media, with practical MATLAB code provided.

Contribution

It introduces a recursive transfer-matrix method for precise energy calculations in multilayered spheres, including special cases and magnetic media.

Findings

Analytic expressions for energy density inside and outside multilayered spheres.

Explicit formulas for lossless core and shell configurations.

Demonstrations of electric field enhancement in silica-gold multilayered spheres.

Abstract

We obtain exact analytic expressions for (i) the electromagnetic energy radial density within and outside a multilayered sphere and (ii) the total electromagnetic energy stored within its core and each of its shells. Explicit expressions for the special cases of lossless core and shell are also provided. The general solution is based on compact recursive transfer-matrix method and its validity includes also magnetic media. The theory is illustrated on the examples of electric field enhancement within various metallo-dielectric silica-gold multilayered spheres. User-friendly MATLAB code which includes the theoretical treatment, is available as a supplement to the paper.

Equations

$${\bf E}=\frac{ic}{\omega\varepsilon}\,\left(\mbox{\boldmath$\nabla$}\times{\bf H}\right)\ ,\qquad{\bf H}=-\frac{ic}{\omega\mu}\,\left(\mbox{\boldmath$\nabla$}\times{\bf E}\right)\ ,$$

$$\mbox{\boldmath$\nabla$}\times\left[\mbox{\boldmath$\nabla$}\times{\bf F}_{\gamma L}(k,{\bf r})\right]=k^{2}{\bf F}_{\gamma L}(k,{\bf r})$$

$$\begin{array}{rl}\mathbf{F}& {}_{ML}(k_n,\mathbf{r})=f_{ML}(k_{n}r){\mathbf{Y}}_L^{(m)}(\mathbf{r}),\\ \mathbf{F}& {}_{EL}(k_n,\mathbf{r})=\\ & \frac{1}{k_{n}r}\left\{\sqrt{l(l+1)}{f}_{EL}(k_{n}r){\mathbf{Y}}_L^{(o)}(\mathbf{r})+\frac{\mathrm{d}}{\mathrm{d}r}\left[r{f}_{EL}(k_{n}r)\right]{\mathbf{Y}}_L^{(e)}(\mathbf{r})\right\},\end{array}$$

$$\begin{split}\tilde{\bf F}_{EL}(k_{n},{\bf r})&=\frac{1}{k_{n}}\mbox{\boldmath$\nabla$}\times{\bf F}_{ML}(k_{n},{\bf r})\ ,\\ \tilde{\bf F}_{ML}(k_{n},{\bf r})&=\frac{1}{k_{n}}\mbox{\boldmath$\nabla$}\times{\bf F}_{EL}(k_{n},{\bf r})\ ,\end{split}$$

$$\begin{array}{rl}{\mathbf{E}}_{\gamma }(\mathbf{r})=& \sum _L{\mathbf{F}}_{\gamma L}(k_n,\mathbf{r})\\ =& \sum _{\gamma ,L}\left[A_{\gamma L}(n){\mathbf{J}}_{\gamma L}(k_n,\mathbf{r})+B_{\gamma L}(n){\mathbf{H}}_{\gamma L}(k_n,\mathbf{r})\right],\end{array}$$

$$f_{\gamma L}=A_{\gamma L}(n)j_l(k_{n}r)+B_{\gamma L}(n)h_l^{(1)}(k_{n}r)$$

$$C_{\gamma L}\equiv A_{\gamma L}(N+1),D_{\gamma L}\equiv B_{\gamma L}(N+1).$$

$$\left(\begin{array}{c}{A}_{\gamma L}(n+1)\\ B_{\gamma L}(n+1)\end{array}\right)=T_{\gamma l}^{+}(n)\left(\begin{array}{c}{A}_{\gamma L}(n)\\ B_{\gamma L}(n)\end{array}\right),$$

$$\left(\begin{array}{c}{A}_{\gamma L}(n)\\ B_{\gamma L}(n)\end{array}\right)=T_{\gamma l}^{-}(n)\left(\begin{array}{c}{A}_{\gamma L}(n+1)\\ B_{\gamma L}(n+1)\end{array}\right).$$

$${\left[T_{\gamma l}^{+}(n)\right]}^{-1}=T_{\gamma l}^{-}(n),{\left[T_{\gamma l}^{-}(n)\right]}^{-1}=T_{\gamma l}^{+}(n).$$

$$B_{EL}(1)=B_{ML}(1)\equiv 0;$$

$${\mathbf{E}}_0\exp \left(i\mathbf{k}\cdot \mathbf{r}\right)=\sum _L{C}_{\gamma L}{\mathbf{J}}_{\gamma L}(kr),$$

$$C_{ML}=4\pi i^l{\mathbf{E}}_0\cdot {\mathbf{Y}}_L^{(m)\ast }(\mathbf{k}),C_{EL}=4\pi i^{l-1}{\mathbf{E}}_0\cdot {\mathbf{Y}}_L^{(e)\ast }(\mathbf{k}),$$

$$\begin{array}{rl}{C}_{ML}& =i^l\sqrt{(2l+1)\pi }{E}_0{\delta }_{m,\pm 1},\\ C_{EL}& =\pm i^l\sqrt{(2l+1)\pi }{E}_0{\delta }_{m,\pm 1},\end{array}$$

$$T_{Ml}^{-}(n)=-i\left(\begin{array}{cc}\tilde{\eta }{\zeta }_l^{\prime }(x){\psi }_l(\tilde{x})-\tilde{\mu }{\zeta }_l(x){\psi }_l^{\prime }(\tilde{x})& \tilde{\eta }{\zeta }_l^{\prime }(x){\zeta }_l(\tilde{x})-\tilde{\mu }{\zeta }_l(x){\zeta }_l^{\prime }(\tilde{x})\\ -\tilde{\eta }{\psi }_l^{\prime }(x){\psi }_l(\tilde{x})+\tilde{\mu }{\psi }_l(x){\psi }_l^{\prime }(\tilde{x})& -\tilde{\eta }{\psi }_l^{\prime }(x){\zeta }_l(\tilde{x})+\tilde{\mu }{\psi }_l(x){\zeta }_l^{\prime }(\tilde{x})\end{array}\right),$$

$$T_{El}^{-}(n)=-i\left(\begin{array}{cc}\tilde{\mu }{\zeta }_l^{\prime }(x){\psi }_l(\tilde{x})-\tilde{\eta }{\zeta }_l(x){\psi }_l^{\prime }(\tilde{x})& \tilde{\mu }{\zeta }_l^{\prime }(x){\zeta }_l(\tilde{x})-\tilde{\eta }{\zeta }_l(x){\zeta }_l^{\prime }(\tilde{x})\\ -\tilde{\mu }{\psi }_l^{\prime }(x){\psi }_l(\tilde{x})+\tilde{\eta }{\psi }_l(x){\psi }_l^{\prime }(\tilde{x})& -\tilde{\mu }{\psi }_l^{\prime }(x){\zeta }_l(\tilde{x})+\tilde{\eta }{\psi }_l(x){\zeta }_l^{\prime }(\tilde{x})\end{array}\right),$$

$$T_{Ml}^{+}(n)=-i\left(\begin{array}{cc}{\zeta }_l^{\prime }(\tilde{x}){\psi }_l(x)/\tilde{\eta }-{\zeta }_l(\tilde{x}){\psi }_l^{\prime }(x)/\tilde{\mu }& {\zeta }_l^{\prime }(\tilde{x}){\zeta }_l(x)/\tilde{\eta }-{\zeta }_l(\tilde{x}){\zeta }_l^{\prime }(x)/\tilde{\mu }\\ -{\psi }_l^{\prime }(\tilde{x}){\psi }_l(x)/\tilde{\eta }+{\psi }_l(\tilde{x}){\psi }_l^{\prime }(x)/\tilde{\mu }& -{\psi }_l^{\prime }(\tilde{x}){\zeta }_l(x)/\tilde{\eta }+{\psi }_l(\tilde{x}){\zeta }_l^{\prime }(x)/\tilde{\mu }\end{array}\right),$$

$$T_{El}^{+}(n)=-i\left(\begin{array}{cc}{\zeta }_l^{\prime }(\tilde{x}){\psi }_l(x)/\tilde{\mu }-{\zeta }_l(\tilde{x}){\psi }_l^{\prime }(x)/\tilde{\eta }& {\zeta }_l^{\prime }(\tilde{x}){\zeta }_l(x)/\tilde{\mu }-{\zeta }_l(\tilde{x}){\zeta }_l^{\prime }(x)/\tilde{\eta }\\ -{\psi }_l^{\prime }(\tilde{x}){\psi }_l(x)/\tilde{\mu }+{\psi }_l(\tilde{x}){\psi }_l^{\prime }(x)/\tilde{\eta }& -{\psi }_l^{\prime }(\tilde{x}){\zeta }_l(x)/\tilde{\mu }+{\psi }_l(\tilde{x}){\zeta }_l^{\prime }(x)/\tilde{\eta }\end{array}\right),$$

$$x_n=k_n{r}_n,{\tilde{\eta }}_n={\eta }_n/{\eta }_{n+1},{\tilde{x}}_n=x_n/{\tilde{\eta }}_n,{\tilde{\mu }}_n={\mu }_n/{\mu }_{n+1}.$$

$${\mathcal{T}}_{\gamma l}(n)=\prod _{j=n-1}^1{T}_{\gamma l}^{+}(j),{\mathcal{M}}_{\gamma l}(n)=\prod _{j=n}^N{T}_{\gamma l}^{-}(j).$$

$${\left[{\mathcal{T}}_{\gamma l}(N+1)\right]}^{-1}={\mathcal{M}}_{\gamma l}(1),{\left[{\mathcal{M}}_{\gamma l}(1)\right]}^{-1}={\mathcal{T}}_{\gamma l}(N+1).$$

$$D_{\gamma L}/C_{\gamma L}={\mathcal{T}}_{21;\gamma l}(N+1)/{\mathcal{T}}_{11;\gamma l}(N+1).$$

$$W=\sum _{n=1}^N{W}_n=\sum _{n=1}^N{\int }_{r_{n-1}}^{r_n}{w}_n(r)r^2\mathrm{d}r,$$

$$w_n(r)=\frac{1}{4}\oint \left[G_e\left({\epsilon }_n\right){\left|\mathbf{E}(\mathbf{r})\right|}^2+G_m\left({\mu }_n\right){\left|\mathbf{H}(\mathbf{r})\right|}^2\right]\mathrm{d}\Omega .$$

$$\begin{array}{rl}\oint {\left|\mathbf{E}\right|}^2\mathrm{d}\Omega =\sum _{l=1}^{\infty }\sum _{m=-l}^l& [|f_{Mlm}{|}^2+\frac{l+1}{2l+1}|f_{E,l-1,m}{|}^2\\ & +\frac{l}{2l+1}|f_{E,l+1,m}{|}^2],\end{array}$$

$$\begin{array}{rl}\oint {\left|\mathbf{H}\right|}^2\mathrm{d}\Omega =\frac{|{\epsilon }_n|}{|{\mu }_n|}\sum _{l=1}^{\infty }\sum _{m=-l}^l& [|f_{Elm}{|}^2+\frac{l+1}{2l+1}|f_{M,l-1,m}{|}^2\\ & +\frac{l}{2l+1}|f_{M,l+1,m}{|}^2].\end{array}$$

$$\begin{array}{rl}{A}_{\gamma L}(n)=& {\mathcal{M}}_{11;\gamma l}(n)C_{\gamma L}+{\mathcal{M}}_{12;\gamma l}(n)D_{\gamma L}\\ =& C_{\gamma L}\left[{\mathcal{M}}_{11;\gamma l}(n)+{\mathcal{M}}_{12;\gamma l}(n)\frac{{\mathcal{T}}_{21;\gamma l}(N+1)}{{\mathcal{T}}_{11;\gamma l}(N+1)}\right],\end{array}$$

$$\begin{array}{rl}{B}_{\gamma L}(n)=& {\mathcal{M}}_{21;\gamma l}(n)C_{\gamma L}+{\mathcal{M}}_{22;\gamma l}(n)D_{\gamma L}\\ =& C_{\gamma L}\left[{\mathcal{M}}_{21;\gamma l}(n)+{\mathcal{M}}_{22;\gamma l}(n)\frac{{\mathcal{T}}_{21;\gamma l}(N+1)}{{\mathcal{T}}_{11;\gamma l}(N+1)}\right].\end{array}$$

$$A_{\gamma L}(n)=C_{\gamma L}{\bar{A}}_{\gamma l}(n),B_{\gamma L}(n)=C_{\gamma L}{\bar{B}}_{\gamma l}(n).$$

$$f_{\gamma L}=C_{\gamma L}{\bar{f}}_{\gamma l},$$

Full text

††journal: josaa

1 Introduction

Multilayered spherical particles of various sizes and material composition are an important part of modern science and technology due to exceptionally adjustable and extraordinary electromagnetic properties. In this regard, the interaction of the electromagnetic wave with a multilayered spherical particle under plane-wave [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] or general beam [12, 13, 14, 15, 16] illuminations represents a problem of a long-standing interest. The solution of this problem implies the definition of the electromagnetic field within or outside of a sphere which allows to get its absorption, scattering, extinction or other important characteristics. Among these properties, the cycle- and orientation-averaged electric $|{\bf E}|^{2}$ and magnetic $|{\bf H}|^{2}$ fields (in general, electromagnetic energy) within a particular layer (shell) or in the vicinity of a multilayered sphere is of great importance since it defines performance and suitability of a multilayered sphere for a large number of intriguing applications: nonlinear optics [17, 18, 19], lasing [20, 21, 22], heating [23, 24, 25], photocatalysis [26], fluorescence enhancement [27, 28, 29, 30], plasmon-enhanced upconversion [31, 32], energy harvesting and storing [33, 34, 35, 36], surface-enhanced Raman spectroscopy [37, 38, 39], biology and medicine [40, 41, 42, 43, 44].

Thus, a specific attention is drawn to theoretical considerations of the electromagnetic energy within and in a proximity of multilayered spherical particles. This fundamental problem has been thoroughly studied for homogeneous spheres: exact analytic expressions are reported for the electromagnetic energy in dielectric [45, 46], magnetic [47] and chiral [48] spheres. These solutions have been extended for two-layered [20, 49] and three-layered [39] spheres by using the recursive relations, which makes corresponding analytic representation quite cumbersome and difficult to generalize for $N>3$ shells. Alternatively, semi-analytic approach might be used for estimating orientation-averaged local fields [37], which involves analytic representation of electromagnetic field within the spherical particle and its consequent numeric integration within the volume of a particle [50]. However, given that various optimization algorithms [51, 52] are used to find the optimal design of multilayered sphere for a particular application, the development of closed-form analytic expressions becomes highly desirable. Here, we fulfill this need and present rigorous, and, quite importantly, compact analytic solution for (i) the electromagnetic energy radial density within and outside a multilayered sphere and (ii) the total electromagnetic energy within its core and each of its shells.

The paper is organized as follows. In Sec. 2, we provide a brief overview of the recursive transfer-matrix solution of the electromagnetic light scattering from a multilayered sphere, which has been proposed and thoroughly discussed in Ref. [27]; in Sec. 3, within the framework of this formalism, we derive a solution for the electromagnetic energy and its density within the multilayered sphere; in Sec. 4 we provide explicit expressions for specific cases of the electromagnetic energy stored within a core of multilayered sphere, in a lossless shell, and in surrounding medium close to a sphere. Discussion of numerical results for silica-gold multilayered nanospheres is given in Sec. 5. Finally, we draw conclusive remarks in Sec. 6 and provide useful relations and derivations in the Appendix.

2 Recursive Transfer-Matrix Method

Consider a multilayered sphere with $N$ concentric shells as shown in Fig. 1. The sphere core counts as shell with number $n=1$ and the host medium is the $n=N+1$ shell. Each shell is assumed to be homogeneous and isotropic with scalar permittivity $\varepsilon_{n}$ and permeability $\mu_{n}$. The outer radius of the $n$-th shell is denoted by $r_{n}$. Spherical coordinates are centered at the sphere origin.

We assume that the multilayered sphere is illuminated with a harmonic electomagnetic wave having vacuum wavelength $\lambda$. Corresponding wave vector in the $n$-th shell is $k_{n}=\eta_{n}\omega/c=2\pi\eta_{n}/\lambda$, where $c$ is the speed of light in vacuum, $\omega$ is frequency, and $\eta_{n}=\sqrt{\varepsilon_{n}\mu_{n}}$ is the refractive index. Electromagnetic fields in any shell are described by the stationary macroscopic Maxwell’s equations (in Gaussian units, with time dependence $e^{-i\omega t}$ assumed and suppressed throughout the paper):

[TABLE]

where the permittivity $\varepsilon$ and permeability $\mu$ are scalars.

Following the notation of [27], the basis of normalized (the normalization here refers to angular integration) transverse vector multipole fields \mbox{\boldmath\nabla}\cdot{\bf F}_{\gamma L}\equiv 0 that satisfy the vector Helmholtz equation

[TABLE]

for $n$-th shell, $1\leq n\leq N+1$, can be formed as [27]:

[TABLE]

where

[TABLE]

with $\tilde{f}_{EL}=f_{ML}$ and $\tilde{f}_{ML}=f_{EL}$. Here $L=lm$ is a composite angular momentum index; ${\bf Y}^{(m)}_{L}$, ${\bf Y}^{(o)}_{L}$ and ${\bf Y}^{(e)}_{L}$ are magnetic, longitudinal, and electric vector spherical harmonics of degree $l$ and order $m$ (see Appendix 6.1 for their definition), and $f_{\gamma L}$ is a suitable linear combination of spherical Bessel functions. Provided that the multipole fields in (2) represent ${\bf E}$, the respective subscripts $M$ and $E$ denote the magnetic, or transverse electric (TE), and electric, or transverse magnetic (TM), polarizations [53].

In the respective cases that $f_{\gamma l}=j_{l}$ and $f_{\gamma l}=h_{l}^{(1)}$, where $j_{l}$ and $h_{l}^{(1)}$ being the spherical Bessel functions of the first and third kind, correspondingly, the multipoles ${\bf F}_{\gamma L}$ will be denoted as ${\bf J}_{\gamma L}$ and ${\bf H}_{\gamma L}$. General solution for the electric field in the $n$-th shell, $1\leq n\leq N+1$, reads then as [27]:

[TABLE]

with corresponding

[TABLE]

to be determined. The expansion of magnetic field H is related to that of the electric field E by the stationary macroscopic Maxwell’s equations (1) on using relations (3).

Expansion coefficients $A_{\gamma L}(n)$ and $B_{\gamma L}(n)$ are determined by matching fields across the shell interfaces, i.e. requiring that the tangential components of ${\bf E}$ and ${\bf H}$ are continuous. For surrounding medium, i.e. the $(N+1)$-th shell, the expansion coefficients will occasionally be written as

[TABLE]

Coefficients $A_{\gamma L}(n)$ and $B_{\gamma L}(n)$ can be found via transfer matrix solution in terms of $2\times 2$ lowering (backward) $T_{\gamma l}^{-}(n)$ and raising (forward) $T_{\gamma l}^{+}(n)$ transfer matrices [27]. The raising transfer matrices translate the expansion coefficients $A_{\gamma L}(n)$ and $B_{\gamma L}(n)$ from the $n$-th shell into the coefficients $A_{\gamma L}(n+1)$ and $B_{\gamma L}(n+1)$ in the $(n+1)$-th shell:

[TABLE]

whereas the lowering transfer matrices translate the coefficients $A(n+1)$ and $B(n+1)$ in the $(n+1)$-th shell into the coefficients $A(n)$ and $B(n)$ in the $n$-th shell:

[TABLE]

It can be easily seen from (5) and (6) that the lowering and raising transfer matrices are related as:

[TABLE]

Provided that the coefficients $A_{\gamma L}(n+1)$ and $B_{\gamma L}(n+1)$ are known, the coefficients $A_{\gamma L}(n)$ and $B_{\gamma L}(n)$ can be unambiguously determined, and vice versa. The constituent transfer matrices $T_{\gamma l}^{+}$ and $T_{\gamma l}^{-}$ can be viewed as analagous to ladder operators of quantum mechanics.

To determine the expansion coefficients at any shell, mixed boundary conditions are imposed which fix two of the coefficients $A_{\gamma L}(n)$ and $B_{\gamma L}(n^{\prime})$ for each given $\gamma$ and $L$, where in general $n\neq n^{\prime}$:

The regularity condition of the solution at the sphere origin:

[TABLE] 2. 2.

For a given frequency $\omega$, the $A_{\gamma L}(N+1)=C_{\gamma L}$ coefficients are equal to the expansion coefficients of an incident electromagnetic field in spherical coordinates.

For general (e.g. focused or Gaussian) beams numerical integration is, as a rule, required to arrive at the expansion coefficients [12, 13, 15, 16]. Nevertheless, closed analytic expressions for the $C_{\gamma L}=0$’s are known for two important incident fields. The most familiar example, which we shall examine in detail below, is furnished by incident plane electromagnetic wave ${\bf E}({\bf r})={\bf E}_{0}\exp\left({i{\bf k}\cdot{\bf r}}\right)$, which expansion in vector spherical wave functions reads as

[TABLE]

where

[TABLE]

and asterisk denotes complex conjugate. In particular, for the plane wave incident along the $z$-axis and polarized along the $x$-axis (i.e. parallel to $\hat{\bf e}_{\vartheta}$ for $\varphi=0$):

[TABLE]

where $\delta_{mm^{\prime}}$ is the Kronecker delta function. For an incident dipole field see Ref. [27].

Irrespective of the incident field, one has the following explicit expressions for the constituent backward and forward transfer matrices [27]:

{strip}

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\psi_{l}(x)=xj_{l}(x)$ and $\zeta_{l}(x)=xh^{(1)}_{l}(x)$ are the Riccati-Bessel functions, prime denotes the derivative with respect to the argument in parentheses, and

[TABLE]

For the sake of clarity, the $n$-subscript has been suppressed in (11) – (14). The above relations for $T^{-}_{\gamma l}(n)$ and $T^{+}_{\gamma l}(n)$ are general and valid for any homogeneous and isotropic medium, including magnetic materials with $\mu_{n}\neq 1$.

It occurs that the formalism becomes compact if one introduces the composite transfer matrices ${\cal T}_{\gamma l}(n)$ and ${\cal M}_{\gamma l}(n)$ defined as ordered (from the left to the right) products of the constituent forward and backward $2\times 2$ matrices:

[TABLE]

Composite matrices ${\cal T}_{\gamma l}(n)$ and ${\cal M}_{\gamma l}(n)$ transfer expansion coefficients from the sphere core to the $n$-th shell, and from the surrounding medium to the $n$-th shell, respectively. Note that ${\cal T}_{\gamma l}(n)$ is defined for $2\leq n\leq N+1$, while ${\cal M}_{\gamma l}(n)$ is defined for $1\leq n\leq N$. Analogously to (7), the following relations can be applied:

[TABLE]

Note that the regularity condition given by (8), unambiguously determines the $m$-independent ratio $D_{\gamma L}/C_{\gamma L}$ [53]:

[TABLE]

Here ${\cal T}_{ij;\gamma l}(n)$ denotes the $(i,j)$-th element of the $2\times 2$ matrix ${\cal T}_{\gamma l}(n)$.

Thus far, the electromagnetic field anywhere inside and outside a multilayered sphere is unambiguously determined from a pair of expansion coefficients $A(n)$ and $B(n)$ for the respective $n$-th shell (including the host medium denoted as $(N+1)$-th shell).

3 Energy within a multilayered sphere

The energy here will have the usual meaning of instant power integrated over a cycle of harmonic excitation [54, 55, 56, 57]. The total electromagnetic energy $W$ within the multilayered sphere can be obtained by integrating the electromagnetic energy radial density $w_{n}(r)$ within each $n$-th shell and, consequently, summing up the total electromagnetic energies $W_{n}$ stored in each $n$-th shell:

[TABLE]

where $r_{0}=0$ for the core, and

[TABLE]

Here $G_{e}\left(\varepsilon_{n}\right)={\rm Re}\left(\varepsilon_{n}\right)$ and $G_{m}\left(\mu_{n}\right)={\rm Re}\left(\mu_{n}\right)$ in the non-dispersive case [46]. For dispersive and absorbing metallic shells $G_{m}$ remains the same, while $G_{e}$ can be described by Loudon’s formula [54]: $G_{e}=[{\rm Re}(\varepsilon_{n})+2\omega{\rm Im}(\varepsilon_{n})/\Gamma_{n}]$, with $\Gamma_{n}$ being the free electron damping constant in the Drude formula.

Angular integration of $\left|{\bf E}\right|^{2}$ and $\left|{\bf H}\right|^{2}$ in (18) in any shell can be performed as follows (see Appendix 6.2 for details):

[TABLE]

[TABLE]

The spherical Bessel functions of the order $l\pm 1$ here originate from eliminating the radial derivation in ${\bf F}_{EL}$ of Eq. (2) by the identity (41) of Appendix 6.2.

Below we present exact analytic expressions for each of the above quantities, i.e. $w_{n}(r)$, $W_{n}$, and $W$. There are two important steps. The first is the summation over $m$ in (19) and (20), which we perform by (i) factorizing $f_{\gamma L}$, and then (ii) using sum rules for $C_{\gamma L}$. The second is a radial integration in (17), which we implement by using Lommel’s integration formulas of Appendix 6.3.

3.1 Factorization of expansion coefficients

Expansion coefficients $A_{\gamma L}(n)$ and $B_{\gamma L}(n)$ can be reformulated via forward ${\cal T}_{\gamma l}(n)$ and backward ${\cal M}_{\gamma l}(n)$ composite transfer matrices as follows:

[TABLE]

[TABLE]

Note that in this representation, coefficients $A_{\gamma L}$ and $B_{\gamma L}$ are products of $m$-dependent $C_{\gamma L}$ and $m$-independent expression in square brackets. Due to the linearity of the equations and spherical symmetry of the problem (the latter being reflected in $m$-independent transfer matrix elements), each of the expansion coefficients $A_{\gamma L}$ and $B_{\gamma L}$ can be factorized as a product of a $m$-dependent factor resulting from $m$-dependence of the expansion coefficients $C_{\gamma L}$, and a $m$-independent factor coming from the transfer matrices:

[TABLE]

Hence each $f_{\gamma L}$ factorizes as

[TABLE]

where $m$-independent $\bar{f}_{\gamma l}$ is explicitly defined as

[TABLE]

Note that in (19) and (20) one has to consider the radial solutions as indexed by $\gamma$, $l$, and $l\pm 1$,

[TABLE]

Since the spherical Bessel functions of the order $l\pm 1$ in (19) and (20) originate from eliminating the radial derivation in ${\bf F}_{EL}$ of (2) according to the identity (41) of Appendix 6.2, they are multiplied by the expansions coefficients $\bar{A}_{\gamma l}$ and $\bar{B}_{\gamma l}$ carrying the index $l$.

After the factorization given by (23), integral in (19) reads as:

[TABLE]

Note that (20) is factorized in a similar manner.

3.2 Sum rules

One can eliminate the $m$-dependence in (26) recalling (10) for $C_{\gamma lm}$ and using the sum rules [58, 59] for magnetic ${\bf Y}^{(m)}_{L}({\rm r})$ and electric ${\bf Y}^{(e)}_{L}({\rm r})$ vector spherical harmonics:

[TABLE]

where $\otimes$ denotes the tensor product. Thus for plane wave incidence the $m$-dependence in (26) yields

[TABLE]

which is also applicable to a factorized representation of (20). For an incident dipole field see Ref. [27].

3.3 Electromagnetic energy radial density

After the factorization of (19) and (20), and subsequent summation over the magnetic number $m$, we end up with the following expressions for the electric

[TABLE]

and for magnetic

[TABLE]

components of the electromagnetic field.

Thus, the electromagnetic energy radial density in (18) is explicitly defined with analytic expressions in (27) and (28).

3.4 Total electromagnetic energy

Finally, the radial integration of (27) and (28) in (17) can be performed by using Lommel’s integration formulas (see Appendix 6.3 for details):

[TABLE]

[TABLE]

Here $x=k_{n}r$, and purely imaginary functions

[TABLE]

are cancelled by purely imaginary $x^{2}-x^{*2}=4i{\rm Re}(x){\rm Im}(x)$ in the denominator, which results in purely real integrals in (29) and (30).

Substitution of (29) and (30) into (17) and (18) yields in explicit expression for the total electromagnetic energy $W_{n}$ stored within each shell of the multilayered sphere. The above relations are general and valid within any shell including $(N+1)$-th layer being a surrounding medium.

3.5 Normalized electromagnetic energy

In some cases, it is of practical use to estimate normalized electromagnetic energy instead of its absolute value. For example, the electromagnetic energy enhancement determines the performance of the spherical particle in surface-enhanced Raman spectroscopy [37, 39], plasmon-enhanced upconversion [32], or might be important in other cases [60, 61, 62]. To get the corresponding enhancement factor, one could compare the energy stored within the $n$-th layer compared to the energy stored in a lossless host medium of the same volume. Given that in a homogeneous medium $|{\bf H}|^{2}=(|\varepsilon_{h}|/|\mu_{h}|)|{\bf E}|^{2}$, the angularly integrated electromagnetic energy density (18) of the incident wave reduces in the lossless host medium to

[TABLE]

where $w_{0}^{(e)}$ and $w_{0}^{(m)}$ represent electric and magnetic components of the $w_{0}$, respectively. These quantities can be used to normalize the electromagnetic energy radial density $w_{n}(r)$ in the presence of a general multilayered sphere at the radial distance $r$ from the sphere origin.

On using (17), the total electromagnetic energy stored within the shell with thickness $\left(r_{n}-r_{n-1}\right)$ characterized by a lossless $\varepsilon_{h}$, is defined as:

[TABLE]

This quantity can be used to normalize the total electromagnetic energy $W_{n}$ stored within $n$-th shell.

3.6 Convergence criterion

For the completeness of the developed theory, it is insightful to provide general remarks on the summation over $l$ in Eqs. (27)–(30). Numerical implementation of these equations requires the truncation to some finite number $l_{\rm max}$, which can be defined for a particular value of the size parameter $x=kr$ with a widely used Wiscombe criterion [63]:

[TABLE]

However, this criterion may vary for near- and far-fields [64].

For large values of the size parameter, one could face convergence issues, because the theoretical treatment for multilayered spherical particles inevitably involves calculation of the difference of the products of the Riccati-Bessel functions. The most successful way to mitigate these issues is to factorize Riccati-Bessel functions with their logarithmic derivatives as shown in Ref. [8].

4 Special Cases

Although presented formalism is rigorous and valid for a general multilayered sphere, it is insightful to discuss some special cases and provide corresponding explicit expressions, which might be more handy to use.

4.1 Core region

For core region, with $n=1$, (21) reduces on using (15) to

[TABLE]

Note that $B_{\gamma L}=0$ due to regularity condition given in (8). Thus, on using (32) and Eqs. (23)–(25):

[TABLE]

We emphasize that $\bar{f}_{\gamma l}$ and $\bar{f}_{\gamma l\pm 1}$ have the same denominator according to (24) and (25).

Thus, the total electromagnetic energy stored inside a core of a general multilayered sphere extending from $r=0$ to $r=r_{1}$ reads as:

[TABLE]

where $\bar{F}^{(1)}_{l}=2i{\rm Im}\left[x_{1}j_{l+1}(x_{1})j_{l}^{*}(x_{1})\right]$.

Known results for the electromagnetic energy within a non-magnetic [46] and magnetic [47] homogeneous sphere can be recovered from (34) by considering the special case of $N=1$.

4.2 Lossless shell

The lossless shell case can be obtained by taking the limit ${\rm Im}(\eta)\to 0$, which yields in $x=x^{*}$, and, as a consequence, vanishing of the denominator in (29) and (30). After applying l’Hôpital’s rule (see Appendix 6.3 for details), (29) and (30) read as:

[TABLE]

[TABLE]

where $r_{0}=0$, and $m$-independent parameter

[TABLE]

Of note, Eqs.(35)–(37) are also valid for a lossless core, though, appropriate expressions for $\bar{f}_{\gamma l}$ and $\bar{f}_{\gamma l+1}$ from (33) have to be used.

4.3 Electromagnetic field in the vicinity of a sphere

Finally, for many applications it is of interest to get the angular averaged electric or magnetic field intensity outside a multilayered particle. Assuming the plane wave incidence and the corresponding plane wave expansion in vector spherical wave functions given by (9), the electric field outside a spherical particle is defined as:

[TABLE]

where $f_{\gamma L}$ in ${\bf F}_{\gamma L}$ are given by [27]:

[TABLE]

Here the $m$-dependent coefficients $C_{\gamma L}$ are given by (10) and (16) for $m$-independent ratio $D_{\gamma L}/C_{\gamma L}$ has been used.

One can now proceed as before when arriving at (27), to yield formally in the same formula, but with $\bar{f}_{\gamma l}$ given in (39). Corresponding enhancement of the average electric field intensity at the distance $r>r_{N}$ from the center of a spherical particle can be obtained after normalization (27) to $4\pi|E_{0}|^{2}$. The magnetic field enhancement near a multilayered sphere can be obtained analogously.

5 Discussion

For the sake of illustration, we apply theory developed here to metallo-dielectric multilayered nanospheres. Appropriately designed by well-developed fabrication procedures [65, 66], such nanoparticles attract significant attention since they may extraordinarily absorb [67, 68, 69], scatter [70] or transmit [71, 72, 73] electromagnetic irradiation, and serve as a platform for photonic bandgap structures [74, 75, 76, 77] or hyperbolic media [78].

Figure 2 shows the normalized electric energy radial density $w_{n}^{(e)}(r)/w^{(e)}_{0}$ (which corresponds to the normalized electric term in (18)) for widely used silica-gold nanospheres. $w_{n}^{(e)}(r)/w^{(e)}_{0}$ is shown for a number of ${\rm SiO}_{2}$-Au spheres of different composition, consisting of $n=1,2,3,5$ or 7 layers with varying thickness. It can be seen that our theory makes it convenient to investigate the electromagnetic field localization within the multilayered spheres. In a particular case of our study, one can observe wavelength-dependent features of the electric field localization within the various layers, depending on the particle composition: $w_{n}^{(e)}(r)/w^{(e)}_{0}$ acquires maximum values at different $\lambda$, from visible to near-IR, and within different layers. Taking as a benchmark the immediate exterior of a homogeneous Au sphere, significant energy density enhancements can be observed both inside and outside of multilayered particles.

6 Conclusion

We have presented a self-consistent rigorous theory for the electromagnetic energy within a general multilayered sphere, which is applicable to a general illumination. Our main focus was on a plane wave illumination for which we obtained exact analytic expressions for (i) the electromagnetic energy radial density within and outside a multilayered sphere and (ii) the total electromagnetic energy stored within its core and each of its shells. Other types of excitation [27, 12, 13, 15, 16, 83] require substitution of corresponding expressions of expansion coefficients $C_{\gamma L}$ for those in (10).

The reported formalism is valid for a wide range of sphere sizes and materials, including magnetic materials. Multilayered spheres from anisotropic [84, 85] or chiral [48] materials can be also considered with the presented formalism after modification.

The theory developed here could have numerous applications. The most straightforward are heating and nonlinear optics applications, which require the determination of $|{\bf E}|^{2}$ or its higher powers. Although only the electric part of the electromagnetic energy has been traditionally extensively considered in the literature, the recent development of all-dielectric nanophotonics also paves a way to a variety of exciting phenomena based on the manipulation of $|{\bf H}|^{2}$ [86, 87, 88], which also can be realized with the multilayered spherical particles of the appropriate composition [89]. In any case, a proper understanding of energy density distribution may provide a valuable insight in many other situations, some of which are the subject of future investigation.

Corresponding MATLAB code, which includes the theoretical treatment reported in this paper, is freely available on GitLab [90]. It is straightforward to modify the code to obtain separate contributions for field intensities, and to determine the stored and dissipated energies for dispersive and absorbing media [54, 55, 56, 57] (cf. also Appendix C of Ref. [27]). It is worthwhile to emphasize that a numerical implementation of the developed theory is of low computational cost due to the utilization of explicit analytic expressions for the electromagnetic energy density.

Appendix

6.1 Vector spherical harmonics

The definition of the vector spherical harmonics from (2) varies in the literature (compare, for instance Bohren [7], Jackson [53], Chew [91], Mishchenko [92]). Here we provide a “combined” representation for magnetic, longitudinal, and electric vector spherical harmonics in spherical $(\varphi,\vartheta,r)$ coordinates:

[TABLE]

where $\hat{\bf e}_{\varphi}$, $\hat{\bf e}_{\vartheta}$, and $\hat{\bf e}_{r}$ are corresponding unit vectors, and $P^{m}_{l}(x)$ are the associated Legendre functions of the first kind [93] of degree $l$ and order $m$:

[TABLE]

6.2 Angular integration of multipole fields

Recalling the fact that vector spherical harmonics are orthonormal, i.e.

[TABLE]

one obtains for the multipole field ${\bf F}_{ML}$ of (2):

[TABLE]

After substituting from [93, Eq. (10.1.22)] into [93, Eq. (10.1.20)] in the recurrence relations for the spherical Bessel functions,

[TABLE]

and the radial derivative in the multipole field ${\bf F}_{EL}$ can be easily eliminated:

[TABLE]

On using the orthonormality of ${\bf Y}^{(o,m,e)}_{L}$, together with [93, Eq. (10.1.19)] applied to the contribution resulting from $\sim{\bf Y}^{(o)}_{L}$ term, (42) is transformed into:

[TABLE]

(40) and (43) are used to perform angular integration of (18), which yields in (19) and (20).

6.3 Lommel’s integration formulas

Consider the defining spherical Bessel equation of given order $l$,

[TABLE]

where $r$ is real in our case, and $k$ is, in general, a complex parameter.

Much the same as in the case of the Lommel formula for the cylindrical Bessel equation [94, Chap. V], the corresponding Lommel formula for two arbitrary solutions $\mathcal{F}_{l}(\rho r)$ and $\mathcal{G}_{l}(\sigma r)$ (here, in general, $\sigma\neq\rho$ are arbitrary complex numbers) of the spherical Bessel equation (44) with different parameter values follows straightforwardly from a fact that

[TABLE]

The latter expression can be recast as

[TABLE]

On integrating both sides one immediately arrives at

[TABLE]

where $\mathcal{C}$ is an integration constant.

In the limit $\rho\to\sigma$, the lhs of (45) goes to zero, whereas the rhs seems to be a nonzero function of $r$. After a closer inspection one finds that the square bracket on the rhs reduces to the Wronskian $W_{r}\{\mathcal{F}_{l}(\sigma r),\mathcal{G}_{l}(\sigma r)\}$. Any solution $\mathcal{F}_{l}$ and $\mathcal{G}_{l}$ of (44) of order $l$ can be expressed as a linear combination of spherical Bessel functions $j_{l}$ and $n_{l}$ with, in general complex, coefficients [cf. (4)]. Subsequently, the Wronskian $W_{r}\{\mathcal{F}_{l},\mathcal{G}_{l}\}$ breaks down into a sum of terms proportional to $W_{r}\{j_{l}(\sigma r),j_{l}(\sigma r)\}$, $W_{r}\{n_{l}(\sigma r),n_{l}(\sigma r)\}$, and $W_{r}\{j_{l}(\sigma r),n_{l}(\sigma r)\}$. The first two factors are identically zero. The last one is proportional to $1/r^{2}$ [93, Eq. (10.1.6)], which cancels the factor $r^{2}$ in front of the square bracket in (45). This means that the term with the square bracket on the rhs of (45) reduces to, in general, nonzero constant $\tilde{\mathcal{C}}$ the limit $\rho\to\sigma$. By taking the integration constant $\mathcal{C}$ in (45) to be just the opposite of the constant $\tilde{\mathcal{C}}$, both sides of (45) go to zero to the limit $\rho\to\sigma$ and the equality is preserved. The point of crucial importance is that for any specific pair of spherical Bessel functions $\mathcal{F}_{l}$ and $\mathcal{G}_{l}$ there is a unique constant $\mathcal{C}$ which ensures equality in (45) including the limit $\rho\to\sigma$. One can thus readily apply l’Hôpital rule to investigate the limit $\rho\to\sigma$ of the expression

[TABLE]

Since the constant $\mathcal{C}$ naturally disappears after applying l’Hôpital’s rule, or when performing definite integrals, it does not impact final expressions and will be omitted below.

If one wants to consider the special case when $\mathcal{G}_{l}=\mathcal{F}_{l}^{*}$, then, as the result of complex conjugation, also the argument of Bessel functions forming $\mathcal{F}_{l}$ becomes complex conjugated. An appropriate application of the Lommel’s first integral (46) to this special case is therefore

[TABLE]

because $\mathcal{F}^{*}_{l}(\sigma r)$ satisfies the Bessel equation (44) with $k=\sigma^{*}$. After using the recurrence relation for Bessel functions [93, Eq. (10.1.22)],

[TABLE]

one arrives at:

[TABLE]

In the special case of $\sigma^{*}=\rho^{*}$:

[TABLE]

where $x=\rho r$. The latter expression in square brackets corresponds to our (31).

In the real limit $\rho^{*}\to\rho$, which corresponds to a lossless core or a lossless shell in the current work, an application of l’Hôpital’s rule and recurrence relations [93, Eqs. (10.1.21-22)] for the Bessel functions in (48) yield:

{strip}

[TABLE]

The latter expression in square brackets corresponds to our (37).