Efficiently computing the emission of a quantum well for non-periodic structures

[oskooi]

Tutorial/An Efficient Approach for Large Numbers of Dipole Emitters described a procedure to simulate the emission of a quantum well (QW) (i.e., spatially incoherent planar source) within a periodic structure (e.g., 1D diffraction grating) using a cosine Fourier series expansion and line sources.

It would be useful to demonstrate how a similar approach can be used to simulate the emission of a QW in a non-periodic structure such as an isolated disc. Alternatively, it would be useful to demonstrate that arranging individual dipole sources non uniformly requires fewer dipole sources (i.e., simulations) to obtain convergence compared to dipole sources arranged uniformly over the same area.

[stevengj]

You can still use a cosine series for a non-periodic case. You can also use any orthonormal basis, I believe — if it's periodic, it might actually more sense to use both cosines and sines (i.e. a full Fourier series), now that I think about it.

With an isolated object such as disc, the trick is to get an orthonormal basis defined over the emitter region. A natural choice might be to use a spectral basis: eigenfunctions of $-\nabla^2$, sorted by eigenvalue (in increasing order), over the domain given by the emitter region with Neumann boundary conditions. (You could even modify this to handle a non-uniform distribution of sources by including a weighting factor in the Laplacian operator). For example, if you do this for a line source you get back the cosine series. For a disc region you would get a tensor product of Bessel functions (in the radial direction) multiplied by a $e^{im\phi}$ Fourier series (in the angular direction). For a general emitter region one could evaluate an eigenbasis computationally with eigs on a finite-difference Laplacian.

[oskooi]

Here is a plot of the extraction efficiency per dipole for a collection of dipoles within a disc. In this example, there are 11 dipoles equally spaced within the QW. The extraction efficiency for all dipoles is the sum of the extraction efficiencies of the individual dipoles (~93%). A total of 92 simulations in Cylindrical (2D) coordinates were necessary because each nonaxisymmetric dipole involves a Fourier-series expansion of the fields with azimuthal dependence $e^{im\phi}$ for $m = 0, ..., M$ where $M$ is the cutoff. The second plot shows the number of simulations ($M + 1$) per dipole. The script used to generate these results is here.

The extraction efficiency for the $n$ th dipole at $r_n$ is computed using this expression:

where $P_{rad}$ in the numerator is the radiated flux (obtained by integrating over the solid angle subtended by $\theta$ in [0, $\pi/2$]), $s(r_n)$ is a weighting function necessary to ensure equal flux contribution of the dipole at $r_n$ relative to the dipole at $r_0 \approx 0$, and $P_{dipole}$ in the denominator is the flux emitted by the dipole. Note that the relative contribution of each dipole to the total flux is proportional to its position $r_n$. As such, dipoles close to the middle of the disc contribute less than those positioned closer to the edge of the disc.

The problem involves developing a procedure that can reproduce the first plot using a fewer number of simulations than 92.

[stevengj]

If we ignore the non-smoothness caused by the discretization, in principle you can do better than equally-spaced points and the trapezoidal rule by using Gaussian quadrature, in particular Gauss–Jacobi quadrature with $\alpha = 0$ and $\beta = 1$.

Using Julia to generate the points and weights, I can do:

using FastGaussQuadrature
y, w = gaussjacobi(11, 0, 1)

However, these are not quite right, because they are integrating for $y \in [-1, +1]$, whereas we want integrals from $x\in [0,1]$, or eventually $r \in [0,R]$. We do this in two steps. First, we do a change of variables $x = (1 + y)/2$ to get us to $[0,1]$:

$$\int_{-1}^{+1} f(y) (1 + y) dy = 4\int_{0}^{1} f(2x - 1) x dx$$

So we need to divide the weights by 4:

x = (y .+ 1) ./ 2
w /= 4
@show x
@show w

which gives:

x = [0.025273620397520347, 0.08304161344740518, 0.16917510037718142, 0.27779671510903203, 0.4015027202328608, 0.531862386910416, 0.6599918420853348, 0.7771593929561622, 0.8753807748555569, 0.9479645488728194, 0.9899817195383196]
w = [0.0002659168233111355, 0.0015109802401199453, 0.004165905862920175, 0.008160386678458337, 0.012897840168074392, 0.017382968789545754, 0.020469775747015875, 0.021164855721004412, 0.018905012014296743, 0.01373452278321788, 0.006341835172035366]

which look like:

These are points and weights for computing:

$$\sum_{k} f(x_k) w_k \approx \int_0^1 f(x) x , dx$$

If you want $\int_0^R g(r) r , dr$, then then you should make the change of variables $r = x R$, so your desired integral is:

$$\int_0^R g(r) r , dr = R^2 \int_0^1 g(xR) x , dx \approx R^2 \sum_{k} g(x_k R) w_k$$

In principle, this should give you something that is more accurate than using equally spaced points and a trapezoidal rule … assuming that the functions can be treated as smooth, despite the discretization.

Here is the same calculation but for 6 points:

x = [0.07305432868025885, 0.2307661379699455, 0.44132848122844986, 0.6630153097188457, 0.8519214003315156, 0.9706835728402151]
w = [0.00873830181360952, 0.043955165550508976, 0.09866115089065525, 0.14079255378819888, 0.13554249723151862, 0.07231033072550873]

[stevengj]

However, using Gaussian quadrature in $r$ still requires you to go to high $m$ for each point at the larger radii. That's why I think it would be better to use an orthonormal basis where you sort by radial and angular frequency simultaneously.

The basic principle, following similar math as Yao et al. (2022) is as follows. Suppose you have dipoles uniformly distributed in some region $\mathscr{R}$, and you have a Maxwell-solve operator $M$ that maps a source current vector $j$ in this region to the field vector $e = M j$. You are computing power = $\Re[e^*j]$ or maybe some quadratic function of the fields. This then corresponds to a quadratic function $j^* H j$ of the currents, where $H$ is Hermitian (e.g. $H = (M + M^*)/2$ for $\Re[e^*j]$). Moreover, you want the expectation value of this quadratic function over random currents $j$. Employing the trace trick, we have:

$E[j^* H j] = E[\text{tr}(j^* H j)] = E[\text{tr}(H j j^*)] = \text{tr}(H \, E[j j^*]) = \text{tr}(H)$

since the current correlation matrix is $E[j j^*] = I$ for spatially uncorrelated currents, with unit amplitude.

Now, suppose that instead, we write the current $j$ in some basis $j = B c$ for basis coefficients $c$, and we want to simply average over uncorrelated coefficients: $E[c c^*] = I$. Then, we have:

$E[\text{tr}(c^* B^* H Bc)] = E[\text{tr}(H B c c^* B^*)] = \text{tr}(H B E[c c^*] B^*)] = \text{tr}(H B B^*)$

So for this to work, we just need $B B^* = I$, i.e. an orthonormal basis.

(More generally, if you had a non-uniform distribution of currents $E[j j^*] = W$, e.g. $W$ is a diagonal matrix of the mean-square current at each point, then you would want a weighed orthogonal matrix $B^* W^{-1} B = I \Longleftrightarrow B B^* = W$.)

[stevengj]

Let's take the case of a disc with radius $R$, and suppose that we only need scalar basis functions for simplicity (e.g. corresponding to currents $J_z$ only). One simple way to get orthonormal functions is to consider eigenfunctions of $-\nabla^2$, say with Neumann (zero-slope) boundary conditions.

In cylindrical coordinates, these eigenfunctions are simply

$$u_{m,n}(r,\phi) = e^{i m\phi} J_m(k_{n,m} r)$$

where we need $J_m'(k_n R) = 0$, so we choose $k_{n,m} R$ to be the $n$-th root of $J_m'(x)$, which you can find tabulated in various places (the root of the slope of the Bessel functions), e.g. it is given by scipy.special.jnp_zeros.

I believe the eigenvalue is then $k_{n,m}^2$, so this is your sort parameter.

(More generally, if we have a weighted distribution of currents with mean-square current $w(x)$, then a corresponding weighted-orthogonal basis could be constructed from eigenfunctions of $-\frac{1}{w} \nabla^2$, which could be computed numerically.)

[stevengj]

Eventually, we want a plot of the accuracy vs. the number of simulations. Starting with the scalar ($J_z$) case seems easiest.

A simple proxy for this is just plotting the contribution vs. the number of basis functions, since it should be rapidly decaying.

[oskooi]

Here are the results for the extraction efficiency (EE) of a QW in the middle of a disc surrounded by vacuum for various dipole arrangements. The configuration involving six nonuniformly spaced dipoles (with weights provided above) with EE obtained using Gaussian quadrature is on the first row with EE of 50.90%. The extraction efficiency for uniformly spaced dipoles (obtained using the trapezoidal rule) seems to be converging to a similar EE value with increasing number of dipoles.

These results seem to indicate that Gaussian quadrature can indeed be used to significantly reduce the number of dipoles while maintaining the accuracy of the EE (in this case reducing the number of dipoles by 15 or a 71% reduction in the number of single-dipole simulations).

Dipole Arrangement	Extraction Efficiency (%)
6 dipoles (nonuniform, Gauss. quad.)	50.90
6 dipoles (nonuniform, trapez. rule)	50.89
6 dipoles (uniform, trapez. rule)	51.07
16 dipoles (uniform, trapez. rule)	50.98
21 dipoles (uniform, trapez. rule)	50.95

[stevengj]

To get the same points and weights in Python, do:

import numpy as np
import scipy.special as ss

n = 11
y, w = ss.roots_jacobi(n, 0, 1)
x = (y + 1) / 2
w = w / 16

(for the n-point rule)

[stevengj]

What you want is a plot error vs number of points $n$ for both Gaussian quadrature and the trapezoidal rule, compared to an "exact" result. For the "exact" result, just use the trapezoidal rule with a large $n$ (I would say at least twice the number of pixels).

I think that the best we can hope for is that Gaussian quadrature is significantly more accurate than trapezoidal for some range of moderately small $n$ values. For large $n$, the non-smoothness of the function (due to the spatial discretization), will spoil the convergence of Gaussian quadrature and make it no better than trapezoidal (and maybe worse). (What this range is will depend on the spatial resolution — Gaussian quadrature should get better at higher resolution where things get smoother.)

Because of the non-smoothness, the gains here may be marginal at typical spatial resolutions. (And if the gains are on the same order as the discretization error, they are not useful.) I still think that bigger gains should arise from the orthonormal-basis approach.