To calculate the transmission spectrum, much as in the bend example in Tutorial/Basics, we'll measure the flux spectrum at one end of the waveguide from a source at the other end, normalized by the flux from a case with no holes in the waveguide. First, we'll load the necessary modules.
import meep as mp
import numpy as np
from matplotlib import pyplot as plt
from IPython.display import Video
%matplotlib notebook
Using MPI version 3.1, 1 processes
Next, we'll define some parameters of our structure. All lengths are in units of microns ($\mu$m). The periodicity of the photonic crystal is 1 $\mu$m.
resolution = 20 # pixels/um
eps = 13 # dielectric constant of waveguide
w = 1.2 # width of waveguide
r = 0.36 # radius of holes
d = 1.4 # defect spacing (ordinary spacing = 1)
N = 3 # number of holes on either side of defect
sy = 6 # size of cell in y direction (perpendicular to wvg.)
pad = 2 # padding between last hole and PML edge
dpml = 1 # PML thickness
Given these parameters, the size of the cell in the X direction, which we'll denote sx
, is given by:
sx = 2 * (pad + dpml + N) + d - 1 # size of cell in x direction
Now, the computational cell is:
cell = mp.Vector3(sx, sy, 0)
The geometry
will consist of a single Block
for the waveguide, and 2N
cylindrical holes:
blk = mp.Block(size=mp.Vector3(mp.inf, w, mp.inf), material=mp.Medium(epsilon=eps))
geometry = [blk]
for i in range(N):
geometry.append(mp.Cylinder(r, center=mp.Vector3(d / 2 + i)))
geometry.append(mp.Cylinder(r, center=mp.Vector3(-(d / 2 + i))))
To create the holes, we have used a for
loop. Note that the geometry
objects are combined using the append
function. As usual, later objects in geometry
take precedence over earlier objects, so the Cylinder
objects will punch holes through the Block
.
The absorbing boundaries surrounding the computational cell are:
pml_layers = [mp.PML(1.0)]
Now, we'll define a couple of parameters to determine the frequency range to investigate. We already know from our previous calculations that this structure has a $H_z$-polarized band gap for frequencies in the range 0.2 to 0.3, so we'll want to cover this interval.
fcen = 0.25 # pulse center frequency
df = 0.2 # pulse frequency width
The source will now be the usual Gaussian pulse centered at fcen
, located at one edge of the cell just outside the PML, at x = - 0.5 * sx + dpml
. Ideally, we would excite exactly the fundamental mode of the waveguide, but it is good enough to just excite it with a line source. Moreover, since we are interested in the $P$ polarization (electric field in the plane), we will excite it with a $J_y$ current source (transverse to the propagation direction), which is specified as $E_y$:
src = [
mp.Source(
mp.GaussianSource(fcen, fwidth=df),
component=mp.Ey,
center=mp.Vector3(-0.5 * sx + dpml),
size=mp.Vector3(0, w),
)
]
The structure has mirror symmetry planes through the $X$ and $Y$ axes. The source breaks the mirror symmetry through the $Y$ axis, but we still have odd mirror symmetry through the $Z$ axis:
sym = [mp.Mirror(mp.Y, phase=-1)]
Note that we specify the plane by its normal, the $Y$ direction. See also Exploiting Symmetry. Putting all these objects together via the Simulation
object:
sim = mp.Simulation(
cell_size=cell,
geometry=geometry,
boundary_layers=pml_layers,
sources=src,
symmetries=sym,
resolution=resolution,
)
We need to compute the flux spectrum at the other end of the computational cell, after the holes but before the PML:
freg = mp.FluxRegion(
center=mp.Vector3(0.5 * sx - dpml - 0.5), size=mp.Vector3(0, 2 * w)
)
nfreq = 500 # number of frequencies at which to compute flux
# transmitted flux
trans = sim.add_flux(fcen, df, nfreq, freg)
Finally, we can visualize the simulation domain to ensure that everything was coded correctly.
f = plt.figure(dpi=150)
sim.plot2D(ax=f.gca())
plt.show()
----------- Initializing structure... block, center = (0,0,0) size (1e+20,1.2,1e+20) axes (1,0,0), (0,1,0), (0,0,1) cylinder, center = (0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1)
Now, we can run the simulation until the sources have finished plus some additional time to allow the fields to propagate through the structure. As in Tutorial/Basics, we'll use stop_when_fields_decayed
to increment the time in steps of 50 time units (about 13 periods) until $|E_y|^2$ has decayed by at least 1/1000 at the transmission-flux plane.
We'll also animate the fields and generate a video.
f = plt.figure(dpi=150)
animate = mp.Animate2D(f=f, fields=mp.Hz, realtime=False, normalize=True)
sim.run(
mp.during_sources(mp.at_every(0.4, animate)),
until_after_sources=mp.stop_when_fields_decayed(
50, mp.Ey, mp.Vector3(0.5 * sx - dpml - 0.5), 1e-3
),
)
plt.close()
Normalizing field data... field decay(t = 50.025000000000006): 3.95790820471984e-05 / 3.95790820471984e-05 = 1.0 field decay(t = 100.05000000000001): 6.062401894821014e-05 / 6.062401894821014e-05 = 1.0 field decay(t = 150.07500000000002): 1.7734570553698603e-05 / 6.062401894821014e-05 = 0.2925337326918046 field decay(t = 200.10000000000002): 7.468025297259177e-06 / 6.062401894821014e-05 = 0.12318591586016357 field decay(t = 250.125): 5.439633091350728e-06 / 6.062401894821014e-05 = 0.08972735865627277 field decay(t = 300.15000000000003): 4.23291165174129e-06 / 6.062401894821014e-05 = 0.06982235300759884 field decay(t = 350.175): 3.476628845824414e-06 / 6.062401894821014e-05 = 0.05734738320127583 field decay(t = 400.20000000000005): 2.847134687922034e-06 / 6.062401894821014e-05 = 0.04696380638100359 field decay(t = 450.225): 2.3289750979705154e-06 / 6.062401894821014e-05 = 0.03841670576079938 field decay(t = 500.25): 1.9177580551625897e-06 / 6.062401894821014e-05 = 0.031633634464268874 field decay(t = 550.275): 1.5662051283234017e-06 / 6.062401894821014e-05 = 0.025834729460304812 field decay(t = 600.3000000000001): 1.2900383339035116e-06 / 6.062401894821014e-05 = 0.02127932717567875 field decay(t = 650.325): 1.0534058195581604e-06 / 6.062401894821014e-05 = 0.017376047280172294 field decay(t = 700.35): 8.676272594305365e-07 / 6.062401894821014e-05 = 0.014311609069859468 field decay(t = 750.375): 7.087820739725053e-07 / 6.062401894821014e-05 = 0.011691439899060524 field decay(t = 800.4000000000001): 5.837770400453704e-07 / 6.062401894821014e-05 = 0.009629467827662155 field decay(t = 850.4250000000001): 4.808205960298484e-07 / 6.062401894821014e-05 = 0.007931189722684065 field decay(t = 900.45): 3.92681229771858e-07 / 6.062401894821014e-05 = 0.006477320979120792 field decay(t = 950.475): 3.234259653563735e-07 / 6.062401894821014e-05 = 0.005334947615938655 field decay(t = 1000.5): 2.6407774313354585e-07 / 6.062401894821014e-05 = 0.004355992026182594 field decay(t = 1050.525): 2.1750627194785985e-07 / 6.062401894821014e-05 = 0.0035877903794809614 field decay(t = 1100.55): 1.7768960210318136e-07 / 6.062401894821014e-05 = 0.0029310099393934598 field decay(t = 1150.575): 1.4635346961931978e-07 / 6.062401894821014e-05 = 0.002414116915349122 field decay(t = 1200.6000000000001): 1.1953143759707475e-07 / 6.062401894821014e-05 = 0.0019716844853058654 field decay(t = 1250.625): 9.845036591681545e-08 / 6.062401894821014e-05 = 0.001623949840754035 field decay(t = 1300.65): 8.038713471789712e-08 / 6.062401894821014e-05 = 0.0013259948138141452 field decay(t = 1350.6750000000002): 6.621009853390941e-08 / 6.062401894821014e-05 = 0.0010921430100249761 field decay(t = 1400.7): 5.408528044424668e-08 / 6.062401894821014e-05 = 0.0008921427741445288 run 0 finished at t = 1400.7 (56028 timesteps)
filename = "media/hole-wvg-cavity.mp4"
animate.to_mp4(10, filename)
Video(filename)
Generating MP4...
We can see the pulse propagating to the right, bouncing off of the holes, and also exciting a resonant mode in the cavity that sits in the center for a long time as it starts slowly leaking to the right.
Of course, the main point of this section is to get the quantitative transmission spectrum. To do this, we need to normalize our flux by running the simulation with no holes. We'll first wrap the above code in a function so that we can parameterize the simulation domain.
def sim_cavity(N=3, sy=6):
sx = 2 * (pad + dpml + N) + d - 1 # size of cell in x direction
cell = mp.Vector3(sx, sy, 0)
blk = mp.Block(size=mp.Vector3(mp.inf, w, mp.inf), material=mp.Medium(epsilon=eps))
geometry = [blk]
for i in range(N):
geometry.append(mp.Cylinder(r, center=mp.Vector3(d / 2 + i)))
geometry.append(mp.Cylinder(r, center=mp.Vector3(-(d / 2 + i))))
src = [
mp.Source(
mp.GaussianSource(fcen, fwidth=df),
component=mp.Ey,
center=mp.Vector3(-0.5 * sx + dpml),
size=mp.Vector3(0, w),
)
]
sim = mp.Simulation(
cell_size=cell,
geometry=geometry,
boundary_layers=pml_layers,
sources=src,
symmetries=sym,
resolution=resolution,
)
freg = mp.FluxRegion(
center=mp.Vector3(0.5 * sx - dpml - 0.5), size=mp.Vector3(0, 2 * w)
)
nfreq = 500
trans = sim.add_flux(fcen, df, nfreq, freg)
sim.run(
until_after_sources=mp.stop_when_fields_decayed(
50, mp.Ey, mp.Vector3(0.5 * sx - dpml - 0.5), 1e-3
)
)
freqs = mp.get_flux_freqs(trans)
psd = mp.get_fluxes(trans)
return freqs, psd
Now we can compute the fluxes for both the simple waveguide and the cavity:
freqs_wg, psd_wg = sim_cavity(N=0) # simple waveguide
freqs_cav, psd_cav = sim_cavity() # cavity
----------- Initializing structure... block, center = (0,0,0) size (1e+20,1.2,1e+20) axes (1,0,0), (0,1,0), (0,0,1) field decay(t = 50.025000000000006): 0.024304124551377593 / 0.024304124551377593 = 1.0 field decay(t = 100.05000000000001): 0.00029478773644369007 / 0.024304124551377593 = 0.012129123837418002 field decay(t = 150.07500000000002): 8.914669657976724e-14 / 0.024304124551377593 = 3.667965755825354e-12 run 0 finished at t = 150.07500000000002 (6003 timesteps) ----------- Initializing structure... block, center = (0,0,0) size (1e+20,1.2,1e+20) axes (1,0,0), (0,1,0), (0,0,1) cylinder, center = (0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) field decay(t = 50.025000000000006): 3.95790820471984e-05 / 3.95790820471984e-05 = 1.0 field decay(t = 100.05000000000001): 6.062401894821014e-05 / 6.062401894821014e-05 = 1.0 field decay(t = 150.07500000000002): 1.7734570553698603e-05 / 6.062401894821014e-05 = 0.2925337326918046 field decay(t = 200.10000000000002): 7.468025297259177e-06 / 6.062401894821014e-05 = 0.12318591586016357 field decay(t = 250.125): 5.439633091350728e-06 / 6.062401894821014e-05 = 0.08972735865627277 field decay(t = 300.15000000000003): 4.23291165174129e-06 / 6.062401894821014e-05 = 0.06982235300759884 field decay(t = 350.175): 3.476628845824414e-06 / 6.062401894821014e-05 = 0.05734738320127583 field decay(t = 400.20000000000005): 2.847134687922034e-06 / 6.062401894821014e-05 = 0.04696380638100359 field decay(t = 450.225): 2.3289750979705154e-06 / 6.062401894821014e-05 = 0.03841670576079938 field decay(t = 500.25): 1.9177580551625897e-06 / 6.062401894821014e-05 = 0.031633634464268874 field decay(t = 550.275): 1.5662051283234017e-06 / 6.062401894821014e-05 = 0.025834729460304812 field decay(t = 600.3000000000001): 1.2900383339035116e-06 / 6.062401894821014e-05 = 0.02127932717567875 field decay(t = 650.325): 1.0534058195581604e-06 / 6.062401894821014e-05 = 0.017376047280172294 field decay(t = 700.35): 8.676272594305365e-07 / 6.062401894821014e-05 = 0.014311609069859468 field decay(t = 750.375): 7.087820739725053e-07 / 6.062401894821014e-05 = 0.011691439899060524 field decay(t = 800.4000000000001): 5.837770400453704e-07 / 6.062401894821014e-05 = 0.009629467827662155 field decay(t = 850.4250000000001): 4.808205960298484e-07 / 6.062401894821014e-05 = 0.007931189722684065 field decay(t = 900.45): 3.92681229771858e-07 / 6.062401894821014e-05 = 0.006477320979120792 field decay(t = 950.475): 3.234259653563735e-07 / 6.062401894821014e-05 = 0.005334947615938655 field decay(t = 1000.5): 2.6407774313354585e-07 / 6.062401894821014e-05 = 0.004355992026182594 field decay(t = 1050.525): 2.1750627194785985e-07 / 6.062401894821014e-05 = 0.0035877903794809614 field decay(t = 1100.55): 1.7768960210318136e-07 / 6.062401894821014e-05 = 0.0029310099393934598 field decay(t = 1150.575): 1.4635346961931978e-07 / 6.062401894821014e-05 = 0.002414116915349122 field decay(t = 1200.6000000000001): 1.1953143759707475e-07 / 6.062401894821014e-05 = 0.0019716844853058654 field decay(t = 1250.625): 9.845036591681545e-08 / 6.062401894821014e-05 = 0.001623949840754035 field decay(t = 1300.65): 8.038713471789712e-08 / 6.062401894821014e-05 = 0.0013259948138141452 field decay(t = 1350.6750000000002): 6.621009853390941e-08 / 6.062401894821014e-05 = 0.0010921430100249761 field decay(t = 1400.7): 5.408528044424668e-08 / 6.062401894821014e-05 = 0.0008921427741445288 run 0 finished at t = 1400.7 (56028 timesteps)
fig = plt.figure(figsize=(11, 6), dpi=100)
ax = fig.add_subplot(111)
plt.plot(freqs_cav, np.array(psd_cav) / np.array(psd_wg), "o-")
plt.grid(True)
plt.xlabel("Frequency")
plt.ylabel("Transmission")
ax2 = fig.add_axes([0.52, 0.6, 0.2, 0.25])
plt.plot(freqs_cav, np.array(psd_cav) / np.array(psd_wg), "o-")
plt.xlim(0.23, 0.24)
plt.ylim(0, 0.8)
plt.grid(True)
The band gap is clearly visible as the range of very low transmission, and in the middle of the band gap is a sharp peak corresponding to the resonant mode trapped in the defect. The inset enlarges this peak, and shows that we didn't use quite enough frequency points to capture the whole shape although we could fit to a Lorentzian if we wanted. At the edges of the band gaps, the transmission goes up in broad Fabry-Perot resonance peaks which we will examine in more detail below. There is also some high-frequency oscillation visible at the left of the plot, which is a numerical artifact due to our pulse not having enough amplitude in that range.
The narrower the resonance peak (higher Q), the harder this sort of direct transmission simulation is to perform — because of the Fourier uncertainty principle, we need to run for a time inversely related to the frequency resolution we would like to obtain. Fortunately, there is a much better way to study high-Q resonances, as described in the next section. See also Tutorial/Basics/Modes of a Ring Resonator.
To study high-Q (long lifetime) resonant modes, it is much more efficient to excite them directly, placing a source inside the cavity, and analyze the resulting fields to obtain the frequencies and lifetimes of the modes. Here, we do precisely that for the above structure. See also Tutorial/Basics/Modes of a Ring Resonator and the Introduction.
The structure is exactly the same as above, and only the sources and analysis are different. We'll create another function that parameterizes the simulation and returns the corresponding simulation object.
This time, the new source is still a Gaussian, but is now a point source at the origin. Moreover, we are now using a magnetic current oriented in the z direction ($H_z$). This source matches the symmetry of the $H_z$-polarized resonant mode that we are looking for. If we didn't know in advance what symmetry we were looking for, we would put the source off-center in a non-symmetric location, which would excite all modes regardless of symmetry. However, in many cases the symmetry is known, and placing a symmetric source allows us to limit the number of modes we excite and also to exploit the fact that we now have two mirror symmetry planes in this problem, saving us a factor of four in computation
def sim_cavity(N=3, sy=6, fcen=0.25, df=0.2):
sx = 2 * (pad + dpml + N) + d - 1 # size of cell in x direction
cell = mp.Vector3(sx, sy, 0)
blk = mp.Block(size=mp.Vector3(mp.inf, w, mp.inf), material=mp.Medium(epsilon=eps))
geometry = [blk]
for i in range(N):
geometry.append(mp.Cylinder(r, center=mp.Vector3(d / 2 + i)))
geometry.append(mp.Cylinder(r, center=mp.Vector3(-(d / 2 + i))))
src = [mp.Source(mp.GaussianSource(fcen, fwidth=df), mp.Hz, mp.Vector3(0))]
sym = [mp.Mirror(mp.Y, phase=-1), mp.Mirror(mp.X, phase=-1)]
sim = mp.Simulation(
cell_size=cell,
geometry=geometry,
boundary_layers=pml_layers,
sources=src,
symmetries=sym,
resolution=resolution,
)
return sim
Here, you may notice a strange thing: we have specified phase=-1
for both mirror planes corresponding to odd symmetry. However, it may seem at first glance that an $H_z$ dipole at the origin has even symmetry! The subtlety here is that the magnetic field is a pseudovector, and is multiplied by −1 under mirror flips, so it is odd when it looks even and vice versa. We aren't just being pedantic here — if you don't realize the difference between vectors, such as electric fields and currents, and pseudovectors, such as magnetic fields and currents, then you will have endless confusion because the electric and magnetic fields will seem to have different symmetry. See also Exploiting Symmetry.
We'll next visualize our domain to make sure everything looks correct:
sim = sim_cavity()
f = plt.figure(dpi=100)
sim.plot2D(ax=f.gca())
plt.show()
----------- Initializing structure... block, center = (0,0,0) size (1e+20,1.2,1e+20) axes (1,0,0), (0,1,0), (0,0,1) cylinder, center = (0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1)
Finally, we can begin the time stepping and animate the fields.
Just as in Tutorial/Basics/Modes of a Ring Resonator, we use the harminv
command (which calls Harminv) to analyze the response at a point (here the $H_z$ field at the origin) for some time after the source has turned off.
h = mp.Harminv(mp.Hz, mp.Vector3(), fcen, df)
sim.run(mp.after_sources(h), until_after_sources=400)
f = plt.figure(dpi=150)
animate = mp.Animate2D(f=f, fields=mp.Hz, realtime=False, normalize=True)
sim.run(mp.at_every(1 / fcen / 20, animate), until=1 / fcen)
plt.close()
harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error harminv0:, 0.23445415345698725, -0.0003147812363599428, 372.4080827818086, 5.812143030902256, -3.763107501450195-4.429450140163555i, 4.306293489743804e-09+0.0i run 0 finished at t = 450.0 (18000 timesteps) Normalizing field data... run 1 finished at t = 454.0 (18160 timesteps)
filename = "media/hole-wvg-cavity-res.mp4"
animate.to_mp4(10, filename)
Video(filename)
Generating MP4...
Inspecting the animation, we see a single resonant mode in the gap. We can verify this by pulling Harminv's results:
f = [m.freq for m in h.modes]
Q = [m.Q for m in h.modes]
for fiter, qiter in zip(f, Q):
print(f"Resonant frequency: {fiter}, Q: {qiter}")
Resonant frequency: 0.23445415345698725, Q: 372.4080827818086
The mode has a frequency of 0.235, just as we saw in the transmission spectrum, and a $Q$ of 373 which we could have also found by fitting the transmission spectrum. This lifetime $Q$ includes two independent decay channels: light can decay from the cavity into the waveguide with lifetime $Q_w$, or it can radiate from the cavity into the surrounding air with lifetime $Q_r$, where
$$\frac{1}{Q} = \frac{1}{Q_w} + \frac{1}{Q_r}
See Chapter 10 of Photonic Crystals: Molding the Flow of Light (second edition) for more details. There are a variety of ways to separate out the two decay channels. For example, we can look at the power radiated in different directions. Here, we'll just increase the number N
of holes and see what happens — as we increase N
, $Q_w$ should increase exponentially while $Q_r$ remains roughly fixed, so that $Q$ eventually saturates at $Q_r$.
N_vec = np.arange(2, 16, 2)
f = []
Q = []
for N in N_vec:
sim = sim_cavity(N=N)
h = mp.Harminv(mp.Hz, mp.Vector3(), fcen, df)
sim.run(mp.after_sources(h), until_after_sources=400)
f.append([m.freq for m in h.modes])
Q.append([m.Q for m in h.modes])
----------- Initializing structure... block, center = (0,0,0) size (1e+20,1.2,1e+20) axes (1,0,0), (0,1,0), (0,0,1) cylinder, center = (0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error harminv0:, 0.2340278559161547, -0.0018734141261924503, 62.4602570900316, 4.626071402369734, -2.8205983802960652-3.6667098872005526i, 4.878861173642782e-09+0.0i run 0 finished at t = 450.0 (18000 timesteps) ----------- Initializing structure... block, center = (0,0,0) size (1e+20,1.2,1e+20) axes (1,0,0), (0,1,0), (0,0,1) cylinder, center = (0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (3.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-3.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error harminv0:, 0.23453222952425626, -7.30923087629776e-05, 1604.356419256047, 6.010010961462772, -3.965258312230556-4.516299178994806i, 2.5784755758812747e-09+0.0i harminv0:, 0.3203380792800332, -0.0028107745362353324, 56.98395142519764, 0.08939729348630548, -0.07611073206285417+0.04689384338623845i, 1.5030156708333963e-06+0.0i run 0 finished at t = 450.0 (18000 timesteps) ----------- Initializing structure... block, center = (0,0,0) size (1e+20,1.2,1e+20) axes (1,0,0), (0,1,0), (0,0,1) cylinder, center = (0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (3.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-3.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (4.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-4.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (5.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-5.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error harminv0:, 0.23454147739266504, -2.0401831532648347e-05, 5748.049556662018, 6.042373374596647, -3.9957221571700208+4.532601950396154i, 2.948523805875737e-09+0.0i harminv0:, 0.32588515566889864, -0.0018520908802034915, 87.97763628993552, 0.08796924097874716, -0.05533533077673982+0.06838558712335283i, 1.0559310146633343e-06+0.0i run 0 finished at t = 450.0 (18000 timesteps) ----------- Initializing structure... block, center = (0,0,0) size (1e+20,1.2,1e+20) axes (1,0,0), (0,1,0), (0,0,1) cylinder, center = (0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (3.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-3.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (4.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-4.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (5.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-5.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (6.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-6.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (7.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-7.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error harminv0:, 0.19490287264916165, -0.0007098826592777302, 137.27823190347084, 0.021882941703098496, -0.012928672972631956-0.017655383109633823i, 2.500625566886593e-06+0.0i harminv0:, 0.23454251746066898, -1.9496236307475725e-05, 6015.071672339521, 6.045187754860183, -3.9989015117156104-4.5335506715057345i, 3.135084966927066e-09+0.0i harminv0:, 0.325392318417585, -0.0028756252341051265, 56.57766432120712, 0.10346319679053909, -0.06988052223995311+0.07629774375162895i, 1.7085893425281652e-06+0.0i run 0 finished at t = 450.0 (18000 timesteps) ----------- Initializing structure... block, center = (0,0,0) size (1e+20,1.2,1e+20) axes (1,0,0), (0,1,0), (0,0,1) cylinder, center = (0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (3.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-3.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (4.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-4.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (5.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-5.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (6.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-6.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (7.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-7.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (8.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-8.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (9.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-9.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) Meep progress: 397.32500000000005/450.0 = 88.3% done in 4.0s, 0.5s to go harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error harminv0:, 0.19654202912381852, -0.0001319711534889467, 744.6401123571296, 0.010139002924627917, -0.009080829279710545-0.0045097582971114i, 3.34141550460465e-06+0.0i harminv0:, 0.23454237393497965, -1.928288781118331e-05, 6081.619522749968, 6.044742625704984, -3.998564131543015-4.533254713442919i, 2.0904183509113547e-09+0.0i harminv0:, 0.32477010784223226, -0.0009459772540947283, 171.65851844557693, 0.10284612856999271, -0.05368440033289738+0.0877229235874676i, 2.3129729952416795e-07+0.0i run 0 finished at t = 450.0 (18000 timesteps) ----------- Initializing structure... block, center = (0,0,0) size (1e+20,1.2,1e+20) axes (1,0,0), (0,1,0), (0,0,1) cylinder, center = (0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (3.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-3.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (4.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-4.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (5.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-5.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (6.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-6.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (7.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-7.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (8.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-8.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (9.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-9.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (10.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-10.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (11.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-11.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) Meep progress: 366.1/450.0 = 81.4% done in 4.0s, 0.9s to go harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error harminv0:, 0.1689088120846072, -0.0008445701020748057, 99.9969165790139, 0.00035136553844277774, -7.403109585265841e-05-0.00034347800286486707i, 4.7246283165373895e-05+0.0i harminv0:, 0.2345423978640612, -1.9291522826319817e-05, 6078.897969217604, 6.044913362697281, -3.9983665938157515-4.5336565974910785i, 1.6528012731386283e-09+0.0i harminv0:, 0.3112255725190105, -0.000515974557591914, 301.5900376672834, 0.014525642975353973, -0.013494086961301489-0.005376236688268417i, 3.824912357074424e-07+0.0i harminv0:, 0.3253120677284504, -0.0007471267911785511, 217.7087420565448, 0.10327042908975875, -0.0574294536238316+0.08582912897642073i, 1.5876314759103416e-07+0.0i run 0 finished at t = 450.0 (18000 timesteps) ----------- Initializing structure... block, center = (0,0,0) size (1e+20,1.2,1e+20) axes (1,0,0), (0,1,0), (0,0,1) cylinder, center = (0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (3.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-3.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (4.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-4.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (5.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-5.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (6.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-6.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (7.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-7.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (8.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-8.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (9.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-9.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (10.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-10.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (11.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-11.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (12.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-12.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (13.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-13.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) Meep progress: 331.45000000000005/450.0 = 73.7% done in 4.0s, 1.4s to go harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error harminv0:, 0.18634339022479032, -0.0005848109413962261, 159.31934325638528, 0.008396135658688452, -0.0015504592340709706-0.008251737402667624i, 9.722082433583947e-06+0.0i harminv0:, 0.23454234757429138, -1.930720606601693e-05, 6073.9587792330785, 6.0453078999997, -3.9988080338829404-4.533793325015162i, 1.6246192346348785e-09+0.0i harminv0:, 0.3251536244537306, -0.0007666444458027871, 212.06285797404294, 0.10960733464915637, -0.0601241424854098+0.09164526883198224i, 1.4407296748106877e-07+0.0i run 0 finished at t = 450.0 (18000 timesteps)
for fiter, qiter in zip(f, Q):
print(f"Resonant frequency: {fiter}, Q: {qiter}")
Resonant frequency: [0.2340278559161547], Q: [62.4602570900316] Resonant frequency: [0.23453222952425626, 0.3203380792800332], Q: [1604.356419256047, 56.98395142519764] Resonant frequency: [0.23454147739266504, 0.32588515566889864], Q: [5748.049556662018, 87.97763628993552] Resonant frequency: [0.19490287264916165, 0.23454251746066898, 0.325392318417585], Q: [137.27823190347084, 6015.071672339521, 56.57766432120712] Resonant frequency: [0.19654202912381852, 0.23454237393497965, 0.32477010784223226], Q: [744.6401123571296, 6081.619522749968, 171.65851844557693] Resonant frequency: [0.1689088120846072, 0.2345423978640612, 0.3112255725190105, 0.3253120677284504], Q: [99.9969165790139, 6078.897969217604, 301.5900376672834, 217.7087420565448] Resonant frequency: [0.18634339022479032, 0.23454234757429138, 0.3251536244537306], Q: [159.31934325638528, 6073.9587792330785, 212.06285797404294]
When we look at the Harminv output for larger N, something strange happens — it starts to find more modes! We'll examine why later. First, let's visualize the Q around the fundamental mode we initially detected ($\omega=0.234$).
idx = np.zeros(N_vec.shape)
Q_fund = np.zeros(N_vec.shape)
for i in range(N_vec.size):
idx[i] = np.abs(np.array(f[i]) - 0.234).argmin()
Q_fund[i] = Q[i][int(idx[i])]
plt.figure(dpi=150)
plt.semilogy(N_vec, Q_fund, "o-")
plt.grid(True)
plt.xlabel("N (# of periods)")
plt.ylabel("Q")
plt.show()
The results, shown above, are exactly what we expected: at first, an exponential increase of $Q$ with N
, and then a saturation at $Q_r≈8750$.
Now, what is this extra mode at $\omega$=0.32823? This is right around the edge of the band gap (actually, just above the edge). There are two possibilities. First, it could be a band edge state: the propagating states in the periodic waveguide go to zero group velocity as they approach the edge of the gap, corresponding to long-lived resonances in a long but finite crystal. Second, it could be a higher-order resonant mode that for a slightly larger defect will be pulled further into the gap, but is currently very delocalized. In this case, it turns out to be the latter.
To see the mode, we will simply run the simulation again with a narrow-band source, and we will also increase the y cell size sy
because it turns out that the mode is fairly spread out in that direction:
sim = sim_cavity(N=16, sy=12, fcen=0.328227374843021, df=0.1)
sim.run(until_after_sources=600)
----------- Initializing structure... block, center = (0,0,0) size (1e+20,1.2,1e+20) axes (1,0,0), (0,1,0), (0,0,1) cylinder, center = (0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-0.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-1.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-2.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (3.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-3.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (4.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-4.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (5.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-5.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (6.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-6.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (7.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-7.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (8.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-8.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (9.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-9.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (10.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-10.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (11.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-11.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (12.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-12.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (13.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-13.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (14.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-14.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (15.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) cylinder, center = (-15.7,0,0) radius 0.36, height 1e+20, axis (0, 0, 1) Meep progress: 184.5/700.0 = 26.4% done in 4.0s, 11.2s to go Meep progress: 363.225/700.0 = 51.9% done in 8.0s, 7.4s to go Meep progress: 547.7/700.0 = 78.2% done in 12.0s, 3.3s to go run 0 finished at t = 700.0 (28000 timesteps)
f = plt.figure(dpi=150)
sim.plot2D(ax=f.gca(), fields=mp.Hz)
plt.show()
From the image, the field is clearly localized around the defect in the center as opposed to being spread out evenly in the crystal like a band-edge state would be. In the defect, the pattern is higher order than the previous mode. It has an extra pair of nodes in the y direction.