## Sensitivity Analysis via Perturbation Theory¶

For a given mode of the ring resonator, it is often useful to know how sensitive the resonant frequency $\omega$ is to small changes in the ring radius $r$ by computing its derivative $\partial\omega/\partial r$. The gradient is also a useful quantity for shape optimization because it can be paired with fast iterative methods such as BFGS to find local optima. The "brute-force" approach for computing the gradient is via a finite-difference approximation requiring two simulations of the (1) unperturbed [$\omega(r)$] and (2) perturbed [$\omega(r+\Delta r)$] structures. Since each simulation is potentially costly, an alternative approach based on semi analytics is to use perturbation theory to obtain the gradient from the fields of the unperturbed structure. This involves a single simulation and is often more computationally efficient than the brute-force approach although some care is required to set up the calculation properly. (More generally, adjoint methods can be used to compute any number of derivatives with a single additional simulation.)

Pertubation theory for Maxwell equations involving high index-contrast dielectric interfaces is reviewed in Chapter 2 of Photonics Crystals: Molding the Flow of Light, 2nd Edition (2008). The formula (equation 30 on p.19) for the frequency shift $\Delta \omega$ resulting from the displacement of a block of $\varepsilon_1$-material towards $\varepsilon_2$-material by a distance $\Delta h$ (perpendicular to the boundary) is:

$$\Delta\omega = -\frac{\omega}{2} \frac{ \iint d^2 \vec{r} \big[ (\varepsilon_1 - \varepsilon_2) |\vec{E}_{\parallel}(\vec{r})|^2 - \big(\frac{1}{\varepsilon_1} - \frac{1}{\varepsilon_2}\big)|\varepsilon\vec{E}_{\perp}|^2\big] \Delta h}{\int d^3\vec{r} \varepsilon(\vec{r})|\vec{E}(\vec{r})|^2} + O(\Delta h^2)$$

In this expression, $\vec{E}_{\parallel}(\vec{r})$ is the component of $\vec{E}$ that is parallel to the surface, and $\varepsilon\vec{E}_{\perp}$ is the component of $\varepsilon\vec{E}$ that is perpendicular to the surface. These two components are guaranteed to be continuous across an interface between two isotropic dielectric materials. In this demonstration, $\partial\omega/\partial r$ is computed using this formula and the results are validated by comparing with the finite-difference approximation: $[\omega(r+\Delta r)-\omega(r)]/\Delta r$.

There are three parts to the calculation: (1) find the resonant frequency of the unperturbed geometry using a broadband Gaussian source, (2) find the resonant mode profile of the unperturbed geometry using a narrowband source and from these fields compute the gradient via the perturbation-theory formula, and (3) find the resonant frequency of the perturbed geometry and from this compute the gradient using a finite-difference approximation. The perturbation is applied only to the inner and outer ring radii. The ring width is constant. There are two types of modes which are computed in separate simulations using different source polarizations: parallel ($E_z$) and perpendicular ($H_z$) to the interface.

In :
import meep as mp
import numpy as np

resolution = 100        # pixels/um

perpendicular = False   # perpendicular (Hz) or parallel (Ez) source?

if perpendicular:
src_cmpt = mp.Hz
fcen = 0.21         # pulse center frequency
else:
src_cmpt = mp.Ez
fcen = 0.17         # pulse center frequency

n = 3.4                 # index of waveguide
w = 1                   # ring width
r = 1                   # inner radius of ring
dpml = 2                # thickness of PML
m = 5                   # angular dependence

pml_layers = [mp.PML(dpml)]

sr = r + w + pad + dpml        # radial size (cell is from 0 to sr)
dimensions = mp.CYLINDRICAL    # coordinate system is (r,phi,z) instead of (x,y,z)
cell = mp.Vector3(sr)

geometry = [mp.Block(center=mp.Vector3(r + (w / 2)),
size=mp.Vector3(w, mp.inf, mp.inf),
material=mp.Medium(index=n))]

# find resonant frequency of unperturbed geometry using broadband source

df = 0.2*fcen           # pulse width (in frequency)

sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
component=src_cmpt,
center=mp.Vector3(r+0.1))]

sim = mp.Simulation(cell_size=cell,
geometry=geometry,
boundary_layers=pml_layers,
resolution=resolution,
sources=sources,
dimensions=dimensions,
m=m)

h = mp.Harminv(src_cmpt, mp.Vector3(r+0.1), fcen, df)
sim.run(mp.after_sources(h),
until_after_sources=100)

frq_unperturbed = h.modes.freq

sim.reset_meep()

# unperturbed geometry with narrowband source centered at resonant frequency

fcen = frq_unperturbed
df = 0.05*fcen

sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
component=src_cmpt,
center=mp.Vector3(r+0.1))]

sim = mp.Simulation(cell_size=cell,
geometry=geometry,
boundary_layers=pml_layers,
resolution=resolution,
sources=sources,
dimensions=dimensions,
m=m)

sim.run(until_after_sources=100)

deps = 1 - n**2
deps_inv = 1 - 1/n**2

if perpendicular:
para_integral = deps*2*np.pi*(r*abs(sim.get_field_point(mp.Ep, mp.Vector3(r)))**2 - (r+w)*abs(sim.get_field_point(mp.Ep, mp.Vector3(r+w)))**2)
perp_integral = deps_inv*2*np.pi*(-r*abs(sim.get_field_point(mp.Dr, mp.Vector3(r)))**2 + (r+w)*abs(sim.get_field_point(mp.Dr, mp.Vector3(r+w)))**2)
numerator_integral = para_integral + perp_integral
else:
numerator_integral = deps*2*np.pi*(r*abs(sim.get_field_point(mp.Ez, mp.Vector3(r)))**2 - (r+w)*abs(sim.get_field_point(mp.Ez, mp.Vector3(r+w)))**2)

denominator_integral = sim.electric_energy_in_box(center=mp.Vector3(0.5*sr), size=mp.Vector3(sr))

perturb_theory_dw_dR = -frq_unperturbed * numerator_integral / (4 * denominator_integral)

sim.reset_meep()

# perturbed geometry with narrowband source

dr = 0.01

sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
component=src_cmpt,
center=mp.Vector3(r + dr + 0.1))]

geometry = [mp.Block(center=mp.Vector3(r + dr + (w / 2)),
size=mp.Vector3(w, mp.inf, mp.inf),
material=mp.Medium(index=n))]

sim = mp.Simulation(cell_size=cell,
geometry=geometry,
boundary_layers=pml_layers,
resolution=resolution,
sources=sources,
dimensions=dimensions,
m=m)

h = mp.Harminv(src_cmpt, mp.Vector3(r+dr+0.1), fcen, df)
sim.run(mp.after_sources(h),
until_after_sources=100)

frq_perturbed = h.modes.freq

finite_diff_dw_dR = (frq_perturbed - frq_unperturbed) / dr

print("dwdR:, {} (pert. theory), {} (finite diff.)".format(perturb_theory_dw_dR,finite_diff_dw_dR))

-----------
Initializing structure...
time for choose_chunkdivision = 0.000187874 s
Working in Cylindrical dimensions.
Computational cell is 8 x 0 x 0.01 with resolution 100
block, center = (1.5,0,0)
size (1,1e+20,1e+20)
axes (1,0,0), (0,1,0), (0,0,1)
dielectric constant epsilon diagonal = (11.56,11.56,11.56)
time for set_epsilon = 0.00455093 s
-----------
Meep: using complex fields.
Meep progress: 63.24/394.1176452636719 = 16.0% done in 4.0s, 20.9s to go
on time step 12648 (time=63.24), 0.000316286 s/step
Meep progress: 127.035/394.1176452636719 = 32.2% done in 8.0s, 16.8s to go
on time step 25410 (time=127.05), 0.00031345 s/step
Meep progress: 190.945/394.1176452636719 = 48.4% done in 12.0s, 12.8s to go
on time step 38195 (time=190.975), 0.00031288 s/step
Meep progress: 254.66/394.1176452636719 = 64.6% done in 16.0s, 8.8s to go
on time step 50941 (time=254.705), 0.000313832 s/step
Meep progress: 316.92/394.1176452636719 = 80.4% done in 20.0s, 4.9s to go
on time step 63396 (time=316.98), 0.000321182 s/step
Meep progress: 376.46/394.1176452636719 = 95.5% done in 24.0s, 1.1s to go
on time step 75307 (time=376.535), 0.000335853 s/step
harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error
harminv0:, 0.17578499227611236, -4.816700722801396e-05, 1824.7448034707427, 0.24613200861807683, -0.14560581011759283-0.19844372937023874i, 1.1797259436409042e-10+0.0i
run 0 finished at t = 394.12 (78824 timesteps)
-----------
Initializing structure...
time for choose_chunkdivision = 0.000282049 s
Working in Cylindrical dimensions.
Computational cell is 8 x 0 x 0.01 with resolution 100
block, center = (1.5,0,0)
size (1,1e+20,1e+20)
axes (1,0,0), (0,1,0), (0,0,1)
dielectric constant epsilon diagonal = (11.56,11.56,11.56)
time for set_epsilon = 0.00452399 s
-----------
Meep: using complex fields.
Meep progress: 64.11/1237.7535400390625 = 5.2% done in 4.0s, 73.2s to go
on time step 12822 (time=64.11), 0.00031198 s/step
Meep progress: 129.395/1237.7535400390625 = 10.5% done in 8.0s, 68.5s to go
on time step 25882 (time=129.41), 0.0003063 s/step
Meep progress: 194.32/1237.7535400390625 = 15.7% done in 12.0s, 64.4s to go
on time step 38871 (time=194.355), 0.000307977 s/step
Meep progress: 260.045/1237.7535400390625 = 21.0% done in 16.0s, 60.2s to go
on time step 52019 (time=260.095), 0.00030423 s/step
Meep progress: 325.825/1237.7535400390625 = 26.3% done in 20.0s, 56.0s to go
on time step 65178 (time=325.89), 0.000303988 s/step
Meep progress: 390.665/1237.7535400390625 = 31.6% done in 24.0s, 52.0s to go
on time step 78150 (time=390.75), 0.00030838 s/step
Meep progress: 455.2/1237.7535400390625 = 36.8% done in 28.0s, 48.1s to go
on time step 91060 (time=455.3), 0.000309849 s/step
Meep progress: 519.66/1237.7535400390625 = 42.0% done in 32.0s, 44.2s to go
on time step 103958 (time=519.79), 0.000310148 s/step
Meep progress: 584.73/1237.7535400390625 = 47.2% done in 36.0s, 40.2s to go
on time step 116975 (time=584.875), 0.000307307 s/step
Meep progress: 649.64/1237.7535400390625 = 52.5% done in 40.0s, 36.2s to go
on time step 129957 (time=649.785), 0.000308119 s/step
Meep progress: 714.945/1237.7535400390625 = 57.8% done in 44.0s, 32.2s to go
on time step 143025 (time=715.125), 0.000306096 s/step
Meep progress: 780.855/1237.7535400390625 = 63.1% done in 48.0s, 28.1s to go
on time step 156209 (time=781.045), 0.000303408 s/step
Meep progress: 846.395/1237.7535400390625 = 68.4% done in 52.0s, 24.0s to go
on time step 169322 (time=846.61), 0.000305051 s/step
Meep progress: 911.86/1237.7535400390625 = 73.7% done in 56.0s, 20.0s to go
on time step 182419 (time=912.095), 0.000305427 s/step
Meep progress: 977.405/1237.7535400390625 = 79.0% done in 60.0s, 16.0s to go
on time step 195532 (time=977.66), 0.000305057 s/step
Meep progress: 1042.795/1237.7535400390625 = 84.2% done in 64.0s, 12.0s to go
on time step 208616 (time=1043.08), 0.000305738 s/step
Meep progress: 1107.3500000000001/1237.7535400390625 = 89.5% done in 68.0s, 8.0s to go
on time step 221526 (time=1107.63), 0.000309854 s/step
Meep progress: 1172.115/1237.7535400390625 = 94.7% done in 72.0s, 4.0s to go
on time step 234482 (time=1172.41), 0.000308753 s/step
Meep progress: 1237.6100000000001/1237.7535400390625 = 100.0% done in 76.0s, 0.0s to go
run 0 finished at t = 1237.755 (247551 timesteps)
-----------
Initializing structure...
time for choose_chunkdivision = 0.000117064 s
Working in Cylindrical dimensions.
Computational cell is 8 x 0 x 0.01 with resolution 100
block, center = (1.51,0,0)
size (1,1e+20,1e+20)
axes (1,0,0), (0,1,0), (0,0,1)
dielectric constant epsilon diagonal = (11.56,11.56,11.56)
time for set_epsilon = 0.0026741 s
-----------
Meep: using complex fields.
Meep progress: 64.245/1237.7535400390625 = 5.2% done in 4.0s, 73.1s to go
on time step 12849 (time=64.245), 0.000311334 s/step
Meep progress: 128.08/1237.7535400390625 = 10.3% done in 8.0s, 69.3s to go
on time step 25621 (time=128.105), 0.000313206 s/step
Meep progress: 192.4/1237.7535400390625 = 15.5% done in 12.0s, 65.2s to go
on time step 38487 (time=192.435), 0.0003109 s/step
Meep progress: 256.76/1237.7535400390625 = 20.7% done in 16.0s, 61.1s to go
on time step 51365 (time=256.825), 0.00031061 s/step
Meep progress: 321.26/1237.7535400390625 = 26.0% done in 20.0s, 57.1s to go
on time step 64271 (time=321.355), 0.000309956 s/step
Meep progress: 385.66/1237.7535400390625 = 31.2% done in 24.0s, 53.0s to go
on time step 77154 (time=385.77), 0.000310499 s/step
Meep progress: 449.995/1237.7535400390625 = 36.4% done in 28.0s, 49.0s to go
on time step 90026 (time=450.13), 0.000310771 s/step
Meep progress: 514.22/1237.7535400390625 = 41.5% done in 32.0s, 45.0s to go
on time step 102876 (time=514.38), 0.000311289 s/step
Meep progress: 578.95/1237.7535400390625 = 46.8% done in 36.0s, 41.0s to go
on time step 115826 (time=579.13), 0.000308899 s/step
Meep progress: 642.9250000000001/1237.7535400390625 = 51.9% done in 40.0s, 37.0s to go
on time step 128622 (time=643.11), 0.000312613 s/step
Meep progress: 706.915/1237.7535400390625 = 57.1% done in 44.0s, 33.0s to go
on time step 141424 (time=707.12), 0.00031246 s/step
Meep progress: 770.58/1237.7535400390625 = 62.3% done in 48.0s, 29.1s to go
on time step 154164 (time=770.82), 0.000313974 s/step
Meep progress: 834.445/1237.7535400390625 = 67.4% done in 52.0s, 25.1s to go
on time step 166942 (time=834.71), 0.00031304 s/step
Meep progress: 898.195/1237.7535400390625 = 72.6% done in 56.0s, 21.2s to go
on time step 179696 (time=898.48), 0.000313635 s/step
Meep progress: 961.405/1237.7535400390625 = 77.7% done in 60.0s, 17.2s to go
on time step 192338 (time=961.69), 0.000316413 s/step
Meep progress: 1025.53/1237.7535400390625 = 82.9% done in 64.0s, 13.2s to go
on time step 205173 (time=1025.87), 0.000311651 s/step
Meep progress: 1089.78/1237.7535400390625 = 88.0% done in 68.0s, 9.2s to go
on time step 218020 (time=1090.1), 0.000311357 s/step
Meep progress: 1151.805/1237.7535400390625 = 93.1% done in 72.0s, 5.4s to go
on time step 230424 (time=1152.12), 0.000322493 s/step
Meep progress: 1210.695/1237.7535400390625 = 97.8% done in 76.0s, 1.7s to go
on time step 242207 (time=1211.04), 0.000339497 s/step
harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error
harminv0:, 0.17493286736703903, -4.814897914072506e-05, 1816.5791932551942, 1.216470973025523, 1.2069307717323075-0.15205176901081904i, 2.316660865172475e-11+0.0i
run 0 finished at t = 1237.755 (247551 timesteps)
dwdR:, -0.08544696397218933 (pert. theory), -0.08521249090733263 (finite diff.)


There are three things to note. First, there is a built-in function electric_energy_in_box which calculates the integral of $\vec{E}\cdot\vec{D}/2 = \varepsilon|E|^2/2$. This is exactly the integral in the denominator, divided by 2. In cylindrical coordinates $(r,\phi,z)$, the integrand is multiplied by the circumference $2\pi r$, or equivalently the integral is over an annular volume. Second, for the case involving the $H_z$ source, both parallel ($E_{\parallel}=E_{\phi}$) and perpendicular ($\varepsilon E_{\perp}=D_r$) fields are present which must be included in the numerator as separate terms. Field values anywhere in the grid obtained with get_field_point are automatically interpolated; i.e., no additional post-processing is necessary. Third, when comparing the perturbation-theory result to the finite-difference approximation, there are two convergence parameters: the resolution and $\Delta r$. In principle, to demonstrate agreement with perturbation theory, the limit of the resolution should be taken to infinity and then the limit of $\Delta r$ to 0. In practice, this can be obtained by doubling the resolution at a given $\Delta r$ until it is sufficiently converged, then halving $\Delta r$ and repeating.

For an $E_z$ source (parallel to the interface), at a resolution of 100 the results are -0.0854469639 (perturbation theory) and -0.0852124909 (finite difference). Doubling the resolution to 200, the results are -0.0854460732 (perturbation theory) and -0.0852115350 (finite difference). Both results have converged to at least five digits. The relative error at resolution 200 is 0.3%. The mode has a $\omega$ of 0.175 and $Q$ of 1800.

For an $H_z$ source (perpendicular to the interface), at a resolution of 100 the results are -0.0805038572 (perturbation theory) and -0.0798087331 (finite difference). Doubling the resolution to 200, the results are -0.0802028346 (perturbation theory) and -0.0798088015 (finite difference). Both results have converged to at least three digits. The relative error at resolution 200 is 0.5%. The error is larger in this case due to the presence of the discontinuous fields at the dielectric interface which are handled by subpixel smoothing. The mode has a $\omega$ of 0.208 and $Q$ of 1200.

Finally, as reference, the same calculation can be set up in Cartesian coordinates as a 2d simulation. There is one major difference: the mode produced by a point source in 2d is actually the $\cos(m\phi)$ mode, not $\exp(im\phi)$, or equivalently it is the superposition of the $\pm m$ modes. This means that computing the numerator integral does not involve just multiplying the field at a single point on the surface by $2\pi r$ — rather, it is the integral of $\cos^2(m\phi)$ which gives a factor of 1/2. (For non-circular shapes in 2d, the surface integral must be computed numerically.) The simulation script is in examples/perturbation_theory_2d.py. The results are comparable to the cylindrical coordinate case (a 1d calculation) but the 2d simulation is much slower and less accurate at the same grid resolution.