Low-overhead distribution strategy for simulation and optimization of large-area metasurfaces

The vector spherical wavefunctions are solutions of the frequency-domain Maxwell’s equations in homogenous media. Consider a medium with wavenumber $k=\omega\sqrt{\varepsilon}/c$ — the electric field $\textbf{E}(\textbf{x})$ in source-free regions then satisfies:

$\displaystyle\nabla\times\nabla\times\textbf{E}(\textbf{x})-k^{2}\textbf{E}(\textbf{x})=0$

(1)

along with the transversality condition $\nabla\cdot\textbf{E}(\textbf{x})=0$ and the transversality condition $\nabla\cdot\textbf{E}(\textbf{x})=0$.

To construct the vector spherical wave functions, we begin by considering the solutions of the scalar helmholtz equation, which is the scalar counterpart of Eq. 1:

$\displaystyle\nabla^{2}\phi(\textbf{x})+k^{2}\phi(\textbf{x})=0$

(2)

It can be shown, with a straightforward application of the separation of variables method, that the spherical wavefunctions $\phi_{l,m}(k_{b}\textbf{x})$ and $\mathcal{R}\phi_{l,m}(k_{b}\textbf{x})$ defined below are solutions of the scalar helmholtz equation:

	$\displaystyle\phi_{l,m}(k\textbf{x})=h_{l}^{(1)}(kr)P_{l}^{|m|}(\cos\theta)\exp(im\varphi)$		(3a)
	$\displaystyle\mathcal{R}\phi_{l,m}(k\textbf{x})=j_{l}(kr)P_{l}^{|m|}(\cos\theta)\exp(im\varphi)$		(3b)

where $l\in\{0,1,2\dots\}$, $m\in\{-l,-l+1,\dots l-1,l\}$, $(r,\theta,\varphi)$ are the spherical coordinates of the point x, $h_{l}^{(1)}(x)$ is the spherical hankel function of the first kind, $j_{l}(x)$ is the spherical bessel function and $P_{l}^{|m|}(x)$ is the normalized associated legendre polynomial. Note that each solution of the scalar wave function is indexed by two numbers — the orbital index $l$ and the magnetic index $m$. $\mathcal{R}$ indicates whether the wavefunction is well behaved at the origin — $\mathcal{R}\phi_{l,m}(k\textbf{x})$ evaluates to 0 at the origin while $\phi_{l,m}(\textbf{x})$ diverges at the origin. Additionally, $\phi_{l,m}(k\textbf{x})$ captures spherical waves radiating to infinity, while $\mathcal{R}\phi_{l,m}(k\textbf{x})$ captures spherical standing waves (equal superpositions of outgoing and incoming waves) — this can be best appreciated by considering the asymptotic forms of these wavefunctions at large radial coordinates ($kr\gg 1$):

		$\displaystyle\phi_{l,m}(k\textbf{x})\sim\frac{\exp(ikr)}{kr}P_{l}^{|m|}(\cos\theta)\exp(im\varphi)$		(4)
		$\displaystyle\mathcal{R}\phi_{l,m}(k\textbf{x})\sim\frac{[\exp(i(kr-l\pi/2)-\exp(-i(kr-l\pi/2))]}{2ikr}P_{l}^{|m|}(\cos\theta)\exp(im\varphi)$		(5)

For each scalar wave function, we can construct two vector spherical wave functions (corresponding to two independent polarizations) which satisfy Eq. 1 along with the transversality condition $\nabla\cdot\textbf{E}(\textbf{x})=0$. We denote the vector spherical wave functions by $\Phi_{l,m,p}(k_{b}\textbf{x})$ and $\mathcal{R}\Phi_{l,m,p}(k_{b}\textbf{x})$ (where $p\in\{0,1\}$ is the polarization index) — these are then defined by:

	$\displaystyle\begin{bmatrix}\Phi_{l,m,p=0}(k\textbf{x})\\ \mathcal{R}\Phi_{l,m,p=0}(k\textbf{x})\end{bmatrix}=\frac{1}{\sqrt{2l(l+1)}}\nabla\begin{bmatrix}\phi_{l,m}(k\textbf{x})\\ \mathcal{R}\phi_{l,m}(k\textbf{x})\end{bmatrix}\times\textbf{x}$		(6a)
	$\displaystyle\begin{bmatrix}\Phi_{l,m,p=1}(k\textbf{x})\\ \mathcal{R}\Phi_{l,m,p=1}(k\textbf{x})\end{bmatrix}=\frac{1}{k}\nabla\times\begin{bmatrix}\Phi_{l,m,p=0}(k\textbf{x})\\ \mathcal{R}\Phi_{l,m,p=0}(k\textbf{x})\end{bmatrix}$		(6b)

where $l\in\{1,2,3\dots\}$, $m\in\{-l,-l+1\dots l-1,l\}$ and $p\in\{0,1\}$. The vector spherical harmonics form a basis for the solutions of the vector helmholtz equation (Eq. 1), and an orthogonality condition can be constructed for them on the surface of (an arbitrarily chosen) sphere as delineated in Doicu .

Consider exciting the scatterer with an incident electric field $\textbf{E}_{\text{inc}}(\textbf{x})$ that is a superposition of regular vector spherical wavefunctions $\mathcal{R}\Phi_{l,m,p}(k_{b}\textbf{x})$, where $k_{b}=\omega\sqrt{\varepsilon_{b}}/c$:

$\displaystyle\textbf{E}_{\text{inc}}(\textbf{x})=\sum_{l,m,p}a_{l,m,p}\mathcal{R}\Phi_{l,m,p}(k_{b}\textbf{x})$

(7)

This incident field interacts with the scatterer, resulting in a scattered field $\textbf{E}_{\text{sca}}(\textbf{x})$ radiating away from the scatterer. Outside the sphere $\partial S$, the scattered field can be expressed as a superposition of vector spherical wavefunctions $\Phi_{l,m,p}(k_{b}\textbf{x})$:

$\displaystyle\textbf{E}_{\text{sca}}(\textbf{x})=\sum_{l,m,p}s_{l,m,p}\Phi_{l,m,p}(k_{b}\textbf{x})$

(8)

Note that the scattered fields can be entirely captured without including the contribution of regular vector spherical wave functions — this is a consequence of the fact that the scattered field is expected to propagate radially outward from the scatterer, and the regular vector spherical harmonics have a radially inward propagating component (Eq. 4). Moreover, the linearity of Maxwell’s equations implies that the coefficients $s_{l,m,p}$ must be related to $a_{l,m,p}$ via a linear transformation:

$\displaystyle s_{l,m,p}=\sum_{l^{\prime},m^{\prime},p^{\prime}}T_{l,m,p;l^{\prime},m^{\prime},p^{\prime}}a_{l^{\prime},m^{\prime},p^{\prime}}$

(9)

In practice, the expansions of the incident and scattered fields Eqs. 7 and 8 are typically truncated to a finite number of terms by ignoring contributions of basis functions with $l>l_{\text{max}}$ [this corresponds to using $2l_{\text{max}}(l_{\text{max}}+2)$ basis functions]. In this case, all the coefficients $a_{l,m,p}$ ($s_{l,m,p}$) can be collected into a vector a (s), and Eq. 9 can be expressed as a matrix equation $\textbf{s}=\textbf{T}\textbf{a}$. We will refer to T as the transition matrix of the scatterer.

There are a number of approaches to compute the transition matrix elements $T_{l,m,p;l^{\prime},m^{\prime},p^{\prime}}$ — we implement the approach based on the null field method Doicu . The null field method uses the surface integral equations to individually relate the incident and scattered fields to the fields on the scatterer surface, expand all the surface fields on the vector spherical harmonic basis, and then compute the relationship between the incident and scattered fields by eliminating the surface fields from the resulting equations. Final computation of the transition matrix method involves computing two matrices P and Q whose elements are given by the following integrals:

	$\displaystyle P_{l,m,p;l^{\prime},m^{\prime},p^{\prime}}=\int_{\partial S}\bigg{[}\mathcal{R}\Phi_{l^{\prime},m^{\prime},p^{\prime}}(k_{s}\textbf{x})\times\Phi_{l,-m,1-p}(k_{b}\textbf{x})-\frac{k_{s}}{k_{b}}\mathcal{R}\Phi_{l^{\prime},m^{\prime},1-p^{\prime}}(k_{s}\textbf{x})\times\Phi_{l,-m,p}(k_{b}\textbf{x})\bigg{]}\cdot\textbf{n}(\textbf{x})\textrm{d}^{2}\textbf{x}$		(10a)
	$\displaystyle Q_{l,m,p;l^{\prime},m^{\prime},p^{\prime}}=\int_{\partial S}\bigg{[}\mathcal{R}\Phi_{l^{\prime},m^{\prime},p^{\prime}}(k_{s}\textbf{x})\times\mathcal{R}\Phi_{l,-m,1-p}(k_{b}\textbf{x})-\frac{k_{s}}{k_{b}}\mathcal{R}\Phi_{l^{\prime},m^{\prime},1-p^{\prime}}(k_{s}\textbf{x})\times\mathcal{R}\Phi_{l,-m,p}(k_{b}\textbf{x})\bigg{]}\cdot\textbf{n}(\textbf{x})\textrm{d}^{2}\textbf{x}$		(10b)

where $k_{s}$ is the wavenumber in the scatterer material, $k_{b}$ is the wavenumber in the background medium in which the scatterers are embedded, and the transition matrix T is given by:

$\displaystyle\textbf{T}=-\textbf{P}\textbf{Q}^{-1}$

(11)

We use a numerical discretization of 0.0025 microns for evaluating these surface integrals to compute the T-matrix.

Finally, we note that the transition matrix becomes an increasingly accurate description of the scatterer on increasing $l_{\text{max}}$. However, for subwavelength scatterers, a small $l_{\text{max}}$ usually suffices to accurately describe the scattering properties of the scatterer. This is a key advantage of using the transition matrix method for simulating a collection of small scatterers, since the scattered fields from each scatterer can be described with a very small number of unknowns when expanded on the vector spherical wavefunctions. For all results in this paper, we have used $l_{\text{max}}=6$.

Finally, we describe how to simulate an ensemble of scatterers when their individual transition matrices are known. Consider $N$ scatterers with transition matrices $\textbf{T}_{1},\textbf{T}_{2}\dots\textbf{T}_{N}$ located at $\textbf{x}_{1},\textbf{x}_{2}\dots\textbf{x}_{N}$. We assume that the enclosing spheres ($\partial S$ defined in section I.2) of the scatterers do not intersect each other. This system of scatterers is excited with an incident field $\textbf{E}_{\text{inc}}(\textbf{x})$ which can be expanded into regular vector spherical wavefunctions centered at each of the scatterers:

$\displaystyle\textbf{E}_{\text{inc}}(\textbf{x})=\sum_{l,m,p}a^{(0)}_{l,m,p;i}\mathcal{R}\Phi_{l,m,n}(k_{b}(\textbf{x}-\textbf{x}_{i}))$

(12)

where $i\in\{1,2,3\dots N\}$ and $k_{b}$ is the wavenumber of the background medium. Note that the coefficients $a_{l,m,n;i}^{(0)}$ for different $i$ are not independent of each other — in particular, given these coefficients for one choice of $i$, the coefficients for all other $i$ can be constructed via an application of the translation theorem for vector spherical wavefunctions Doicu . Typically, the incident field is known analytically (e.g. plane wave field), and hence these coefficients can be analytically determined for an arbitrary $i$.

We express the scattered field as a superposition of vector spherical harmonics radiated from each scatterer:

$\displaystyle\textbf{E}_{\text{sca}}(\textbf{x})=\sum_{i=1}^{N}\sum_{l,m,p}s_{l,m,p;i}\Phi_{l,m,p}(k_{b}(\textbf{x}-\textbf{x}_{i}))$

(13)

where it is assumed that x lies outside the enclosing spheres of all scatterers. Consider now the $j^{\text{th}}$ scatterer — the total incident field at this scatterer, $\textbf{E}_{{\text{inc}},j}(\textbf{x})$ is given by the sum of $\textbf{E}_{\text{inc}}(\textbf{x})$ as well as fields radiated by the neighbouring scatterers:

$\displaystyle\textbf{E}_{{\text{inc}},j}(\textbf{x})=\textbf{E}_{\text{inc}}(\textbf{x})+\sum_{\begin{subarray}{c}i=1\\ i\neq j\end{subarray}}^{N}\sum_{l,m,p}s_{l,m,p;i}\Phi_{l,m,p}(k_{b}(\textbf{x}-\textbf{x}_{i}))$

(14)

This field can be expressed entirely as a superposition of the regular vector spherical wavefunctions centered at $\textbf{x}_{j}$ by the application of the translation theorem for the vector spherical wavefunctions. The translation theorem relates vector spherical wavefunctions defined with respect to two different origins to each other Doicu ; Egel — more explicitly:

$\displaystyle\Phi_{l,m,p}(k(\textbf{x}-\textbf{x}_{a}))=\sum_{l^{\prime},m^{\prime},p^{\prime}}\xi_{l,m,p;l^{\prime},m^{\prime},p^{\prime}}(k(\textbf{x}_{a}-\textbf{x}_{b}))\Phi_{l^{\prime},m^{\prime},p^{\prime}}(k(\textbf{x}-\textbf{x}_{b}))$

(15)

where

$\displaystyle\xi_{l,m,p;l^{\prime},m^{\prime},p^{\prime}}(\textbf{x})=\delta_{p,p^{\prime}}\alpha_{l,m;l^{\prime},m^{\prime}}(\textbf{x})+(1-\delta_{p,p^{\prime}})\beta_{l,m;l^{\prime},m^{\prime}}(\textbf{x})$

(16)

and with $(r,\varphi,\theta)$ as the spherical coordinates of the vector x:

	$\displaystyle\alpha_{l,m;l^{\prime},m^{\prime}}(\textbf{x})=\exp(i(m-m^{\prime})\varphi)\sum_{q=|l-l^{\prime}|}^{l+l^{\prime}}a_{5}(l,m|l^{\prime},m^{\prime}|q)h_{q}^{(1)}(r)P_{q}^{|m-m^{\prime}|}(\cos\theta)$		(17a)
	$\displaystyle\beta_{l,m;l^{\prime},m^{\prime}}(\textbf{x})=\exp(i(m-m^{\prime})\varphi)\sum_{q=|l-l^{\prime}|+1}^{l+l^{\prime}}b_{5}(l,m|l^{\prime},m^{\prime}|q)h_{q}^{(1)}(r)P_{q}^{|m-m^{\prime}|}(\cos\theta)$		(17b)

where:

	$\displaystyle\alpha(l,m|l^{\prime},m^{\prime}|p)=$	$\displaystyle i^{|m-m^{\prime}|-|m|-|m^{\prime}|+l^{\prime}-l+p}(-1)^{m-m^{\prime}}$		(18a)
		$\displaystyle\times[l(l+1)+l^{\prime}(l^{\prime}+1)-p(p+1)]\sqrt{2p+1}$
		$\displaystyle\times\sqrt{\frac{(2l+1)(2l^{\prime}+1)}{2l(l+1)(l^{\prime}+1)}}\begin{pmatrix}l&l^{\prime}&p\\ m&-m^{\prime}&m^{\prime}-m\end{pmatrix}\begin{pmatrix}l&l^{\prime}&p\\ 0&0&0\end{pmatrix}$
	$\displaystyle\beta(l,m|l^{\prime},m^{\prime}|p)=$	$\displaystyle i^{|m-m^{\prime}|-|m|-|m^{\prime}|+l^{\prime}-l+p}(-1)^{m-m^{\prime}}$		(18b)
		$\displaystyle\times\sqrt{(l+l^{\prime}+1+p)(l+l^{\prime}+1-p)(p+l-l^{\prime})(p-l+l^{\prime})(2p+1)}$
		$\displaystyle\times\sqrt{\frac{(2l+1)(2l^{\prime}+1)}{2l(l+1)(l^{\prime}+1)}}\begin{pmatrix}l&l^{\prime}&p\\ m&-m^{\prime}&m^{\prime}-m\end{pmatrix}\begin{pmatrix}l&l^{\prime}&p\\ 0&0&0\end{pmatrix}$

with

$\displaystyle\begin{pmatrix}j_{1}&j_{2}&j_{3}\\ m_{1}&m_{2}&m_{3}\end{pmatrix}$

(19)

is the Wigner-3j symbol Luscombe . With this addition theorem, $\textbf{E}_{\text{inc},j}(\textbf{x})$ can be rewritten as:

$\displaystyle\textbf{E}_{\text{inc},j}(\textbf{x})=\sum_{l,m,p}a_{l,m,p;j}\mathcal{R}\Phi_{l,m,p}(k_{b}(\textbf{x}-\textbf{x}_{j}))$

(20)

where

$\displaystyle a_{l,m,p;j}=a^{(0)}_{l,m,p;j}+\sum_{\begin{subarray}{c}i=1\\ i\neq j\end{subarray}}^{N}\sum_{l^{\prime},m^{\prime},p^{\prime}}\xi_{l^{\prime},m^{\prime},p^{\prime};l,m,p}(k_{b}(\textbf{x}_{j}-\textbf{x}_{i}))s_{l^{\prime},m^{\prime},p^{\prime};i}$

(21)

Finally, the coefficients $a_{l,m,p;j}$ can be related to the coefficients $s_{l,m,p;j}$ via the transition matrix of the $j^{\text{th}}$ scatterer: $\textbf{s}_{j}=\textbf{T}_{j}\textbf{a}_{j}$, where $\textbf{s}_{j}$ and $\textbf{a}_{j}$ are column vectors of the coefficients $s_{l,m,n;j}$ and $a_{l,m,n;j}$. Denoting by $\boldsymbol{\xi}(\textbf{x}_{i},\textbf{x}_{j})$ as the matrix of the translation coefficients $\xi_{l,m,p;l^{\prime},m^{\prime},p^{\prime}}(k_{b}(\textbf{x}_{i}-\textbf{x}_{j}))$ with the unprimed indices corresponding to different rows and the primed indices corresponding to different columns, this equation can be compactly written as:

$\displaystyle\textbf{T}_{j}^{-1}\textbf{s}_{j}-\sum_{\begin{subarray}{c}i=1\\ i\neq j\end{subarray}}^{N}\boldsymbol{\xi}^{T}(\textbf{x}_{j},\textbf{x}_{i})\textbf{s}_{i}=\textbf{a}^{(0)}_{j}$

(22)

Such an equation can be written for each scatterer — all of these equations collected together provide a completely determined system of linear equations that can be solved to obtain the scattering coefficients $s_{l,m,p;j}$. The system of equations can be written in a matrix form:

$\displaystyle\underbrace{\begin{bmatrix}\textbf{T}_{1}^{-1}&\boldsymbol{\xi}^{T}(\textbf{x}_{1},\textbf{x}_{2})&\boldsymbol{\xi}^{T}(\textbf{x}_{1},\textbf{x}_{3})&\dots&\boldsymbol{\xi}^{T}(\textbf{x}_{1},\textbf{x}_{N})\\ \boldsymbol{\xi}^{T}(\textbf{x}_{2},\textbf{x}_{1})&\textbf{T}_{2}^{-1}&\boldsymbol{\xi}^{T}(\textbf{x}_{2},\textbf{x}_{3})&\dots&\boldsymbol{\xi}^{T}(\textbf{x}_{2},\textbf{x}_{N})\\ \boldsymbol{\xi}^{T}(\textbf{x}_{3},\textbf{x}_{1})&\boldsymbol{\xi}^{T}(\textbf{x}_{3},\textbf{x}_{2})&\textbf{T}_{3}^{-1}&\dots&\boldsymbol{\xi}^{T}(\textbf{x}_{3},\textbf{x}_{N})\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ \boldsymbol{\xi}^{T}(\textbf{x}_{N},\textbf{x}_{1})&\boldsymbol{\xi}^{T}(\textbf{x}_{N},\textbf{x}_{2})&\boldsymbol{\xi}^{T}(\textbf{x}_{N},\textbf{x}_{3})&\dots&\textbf{T}_{N}^{-1}\end{bmatrix}}_{\boldsymbol{\Omega}}\underbrace{\begin{bmatrix}\textbf{s}_{1}\\ \textbf{s}_{2}\\ \textbf{s}_{3}\\ \vdots\\ \textbf{s}_{N}\end{bmatrix}}_{\textbf{s}}=\underbrace{\begin{bmatrix}\textbf{a}^{(0)}_{1}\\ \textbf{a}^{(0)}_{2}\\ \textbf{a}^{(0)}_{3}\\ \vdots\\ \textbf{a}^{(0)}_{N}\end{bmatrix}}_{\textbf{a}^{(0)}}$

(23)

The matrix $\boldsymbol{\Omega}$ in the above system of equations is referred to as the ‘maxwell operator‘, since it is equivalent to expressing the frequency domain maxwell’s equations on the vector spherical wavefunction basis. To summarize, the problem of solving the Maxwell’s equations has been reduced to the problem of solving the linear system in Eq. 23. Once this system of equations is solved to obtain $s_{l,m,n;j}$, the electric fields at the desired points can be computed using Eq. 13, which can be concisely expressed as:

$\displaystyle\textbf{E}_{\text{sca}}(\textbf{x})=\sum_{i=1}^{N}\textbf{s}_{i}^{T}\Phi(\textbf{x}-\textbf{x}_{i})$

(24)

where $\Phi(\textbf{x})$ is a column vector of the vector spherical harmonics $\Phi_{l,m,p}(\textbf{x})$ at position x.

From a numerical perspective, the linear systems Eq. 23 is severely ill conditioned due to the transition matrices being close to singular for small scatterers. This issue can be mitigated using the following block diagonal preconditioner M:

$\displaystyle\textbf{M}=\begin{bmatrix}\textbf{T}_{1}&&&\\ &\textbf{T}_{2}&&\\ &&\textbf{T}_{3}&&\\ &&&\ddots&\\ &&&&\textbf{T}_{N}\\ \end{bmatrix}$

(25)

Instead of solving $\boldsymbol{\Omega}\textbf{s}=\textbf{a}^{(0)}$, we solve $\textbf{M}\boldsymbol{\Omega}\textbf{s}=\textbf{M}\textbf{a}^{(0)}$ (a system of equations similar to those found in markkanen2017fast ). Validation simulations for our T-matrix simulation method implementation are shown in Fig. S2 — we see good agreement between the T-matrix method, FDTD, and FDFD simulations.

In order to implement the Nyquist-sampling-based-parellelization scheme described above, it is necessary to be able to simulate the response of the metasurface to a collection of jinc sources using the T-matrix method. This requires the ability to expand the jinc sources on the vector spherical wavefunctions (Eq. 23). In this section, we mathematically develop such an expansion.

Consider a jinc source at $z=z_{0}$ propagating in the $+z$ direction and centered at $\textbf{x}^{t}_{0}=(x_{0},y_{0})$ in the transverse plane. The transverse electric field at $z=z_{0}$ is given by:

$\displaystyle\textbf{E}^{T}(\textbf{x}^{T},z=z_{0})=\textbf{A}^{T}\frac{j_{1}(k_{0}|\textbf{x}^{T}-\textbf{x}^{T}_{0}|)}{k_{0}|\textbf{x}^{T}-\textbf{x}^{T}_{0}|}$

(27)

where $j_{1}(\cdot)$ is the first order spherical bessel function and $\textbf{A}^{T}=A_{x}\hat{x}+A_{y}\hat{y}$ is the transverse polarization of the jinc field. Alternatively, in the fourier representation,

$\displaystyle\textbf{E}^{T}(x,y,z=z_{0})=\frac{\textbf{A}^{T}}{2\pi}\int_{|\textbf{k}^{T}|<k_{0}}\exp[i\textbf{k}^{T}\cdot(\textbf{x}^{T}-\textbf{x}^{T}_{0})]d^{2}\textbf{k}^{T}$

(28)

Propagating each plane-wave component, we obtain:

$\displaystyle\textbf{E}(x,y,z)=\frac{1}{2\pi}\int_{|\textbf{k}^{T}|<k_{0}}\textbf{A}(\textbf{k}^{T})\exp(i\textbf{k}^{T}\cdot(\textbf{x}^{T}-\textbf{x}^{T}_{0}))\exp(ik_{z}(z-z_{0})){d}^{2}\textbf{k}^{T}$

(29)

where $k_{z}=\sqrt{k_{0}^{2}-\textbf{k}^{T}\cdot\textbf{k}^{T}}$ and

$\displaystyle\textbf{A}(\textbf{k}^{T})=\textbf{A}^{T}-\hat{z}\frac{\textbf{k}^{T}\cdot\textbf{A}^{T}}{k_{z}(\textbf{k}^{T})}$

(30)

Changing the integration variable to $\alpha,\beta$ ($0\leq\alpha\leq 2\pi$ and $0\leq\beta\leq\pi/2$) where $\textbf{k}^{T}(\beta,\alpha)=k_{0}\sin\beta(\hat{x}\cos\alpha+\hat{y}\sin\alpha)$ and therefore $k_{z}(\textbf{k}^{T})\equiv k_{z}(\beta,\alpha)=k_{0}\cos\beta$, $\textbf{A}(\textbf{k}^{T})\equiv\textbf{A}(\beta,\alpha)=\hat{x}A_{x}+\hat{y}A_{y}-\hat{z}\tan\beta(A_{x}\cos\alpha+A_{y}\sin\alpha)$ and $d^{2}\textbf{k}^{T}=k_{0}^{2}\sin\beta\cos\beta d\alpha d\beta$. Moreover, each plane wave can be expanded into a sum of vector spherical harmonic wavefunctions Doicu :

$\displaystyle\textbf{A}(\beta,\alpha)\exp(i\textbf{k}^{T}(\beta,\alpha)\cdot(\textbf{x}^{T}-\textbf{x}^{T}_{0}))\exp(ik_{z}(\beta,\alpha)(z-z_{0}))=\exp[i\textbf{k}(\beta,\alpha)\cdot(\textbf{x}_{0}^{\prime}-\textbf{x}_{0})]\sum_{l,m,p}a_{l,m,p}(\beta,\alpha)\mathcal{R}\Phi_{l,m,p}(k_{0}(\textbf{x}-\textbf{x}_{0}^{\prime}))$

(31)

where

	$\displaystyle a_{l,m,p=0}(\beta,\alpha)=4i^{m}\textbf{A}(\beta,\alpha)\cdot\textbf{m}_{l,m}^{*}(\beta,\alpha)=-\frac{4i^{l}\exp(-im\alpha)}{\sqrt{2l(l+1)}}\big{[}im\pi^{|m|}_{l}(\beta)(\hat{\beta}\cdot\textbf{A}(\beta,\alpha))+\tau^{|m|}_{l}(\beta)(\hat{\alpha}\cdot\textbf{A}(\beta,\alpha))\big{]}$		(32a)
	$\displaystyle a_{l,m,p=1}(\beta,\alpha)=-4i^{l+1}\textbf{A}(\beta,\alpha)\cdot\textbf{n}_{l,m}^{*}(\beta,\alpha)=-\frac{4i^{l+1}\exp(-im\alpha)}{\sqrt{2l(l+1)}}\big{[}\tau^{|m|}_{l}(\beta)(\hat{\beta}\cdot\textbf{A}(\beta,\alpha))-im\pi^{|m|}_{l}(\beta)(\hat{\alpha}\cdot\textbf{A}(\beta,\alpha))\big{]}$		(32b)

where $\textbf{m}_{l,m}(\cdot,\cdot)$ and $\textbf{n}_{l,m}(\cdot,\cdot)$ are vector spherical harmonics with orbital index $l$ and azimuthal index $m$. Therefore,

$\displaystyle\textbf{E}(\textbf{x})=\sum_{l,m,p}a_{l,m,p}\mathcal{R}\Phi_{l,m,p}(k_{0}(\textbf{x}-\textbf{x}_{0}^{\prime}))$

(33)

where $a_{l,m,p}$ are given by:

$\displaystyle a_{l,m,p}=\frac{1}{2\pi}\int_{\beta=0}^{\pi/2}\int_{\alpha=0}^{2\pi}a_{l,m,p}(\beta,\alpha)\exp[i\textbf{k}(\beta,\alpha)\cdot(\textbf{x}_{0}^{\prime}-\textbf{x}_{0})]\sin\beta\cos\beta d\alpha d\beta$

(34)

In the remainder of this section, we will evaluate the integral over $\alpha$ analytically, and the resulting expression can then be numerically integrated over $\beta$. To do so, it will be convenient to define a function $\Gamma_{m}(\xi,\eta,\rho)$ by:

$\displaystyle\Gamma_{m}(\xi,\eta,\rho)=\frac{1}{2\pi}\int_{0}^{2\pi}\exp(-im\alpha)\exp[i\rho\cos(\alpha-\xi)]\cos(\alpha-\eta)d\alpha$

(35)

$\Gamma_{m}(\xi,\eta,\rho)$ can be evaluated analytically: making a change of variables to $\alpha^{\prime}=\alpha-\xi+\pi/2$, we obtain:

	$\displaystyle\Gamma_{m}(\xi,\eta,\rho)$	$\displaystyle=\frac{i^{m}\exp(-im\xi)}{2\pi}\int_{\pi/2-\xi}^{5\pi/2-\xi}\exp(-im\alpha^{\prime})\exp(i\rho\sin\alpha^{\prime})\sin(\alpha^{\prime}+\xi-\eta)d\alpha^{\prime}$
		$\displaystyle=\frac{i^{m-1}\exp(-im\xi)}{4\pi}\int_{0}^{2\pi}\exp(-im\alpha^{\prime})\exp(i\rho\sin\alpha^{\prime})(\exp(i(\alpha^{\prime}+\xi-\eta))-\exp(-i(\alpha^{\prime}+\xi-\eta)))d\alpha^{\prime}$
		$\displaystyle=\frac{i^{m-1}\exp(-im\xi)}{2}\bigg{[}\exp\{i(\xi-\eta)\}J_{m-1}(\rho)-\exp\{-i(\xi-\eta)\}J_{m+1}(\rho)\bigg{]}$		(36)

wherein we have used the identity:

$\displaystyle J_{m}(\rho)=\frac{1}{2\pi}\int_{0}^{2\pi}\exp[i\rho\sin\theta-m\theta]d\theta$

(37)

Additionally, note that $\hat{\beta}=(\hat{x}\cos\alpha+\hat{y}\sin\alpha)\cos\beta-\hat{z}\sin\beta$ and $\hat{\alpha}=-\hat{x}\sin\alpha+\hat{y}\cos\alpha$. Therefore:

	$\displaystyle\hat{\beta}\cdot\textbf{A}(\beta,\alpha)=\sec\beta(A_{x}\cos\alpha+A_{y}\sin\alpha)$		(38a)
	$\displaystyle\hat{\alpha}\cdot\textbf{A}(\beta,\alpha)=-A_{x}\sin\alpha+A_{y}\cos\alpha$		(38b)

Finally, let $\textbf{x}_{0}^{\prime}-\textbf{x}_{0}\equiv(r_{0},\theta_{0},\varphi_{0})$, and therefore $\textbf{k}(\beta,\alpha)\cdot(\textbf{x}_{0}^{\prime}-\textbf{x}_{0})=k_{0}r_{0}(\cos\beta\cos\theta_{0}+\sin\beta\sin\theta_{0}\cos(\alpha-\varphi_{0}))$.