Include cathodoluminescence and EELS examples in meep

[Nikolaos-Matthaiakakis]

Dear all,

I am creating this issue as I am planning to try and include examples of Cathodoluminescence (CL) and Electron Energy Loss Spectroscopy (EELS) simulations in meep. If any other users are interested in this I would be happy to work together to include this useful feature. I am currently quite busy with several projects so it might take a while to get a final result.

In the future, this could potentially be included as a built-in feature including a function for e-beam excitation, and monitors that include the relevant calculations to obtain (for example) the EELS spectra by default.

For this first post I include some relevant information on how this could be implemented based on literature and commercial FDTD packages.

1) https://pubs.acs.org/doi/abs/10.1021/acs.nanolett.5b00716 has useful information about CL FDTD simulations in the supplementary file 2) https://support.lumerical.com/hc/en-us/articles/360042706693-Electron-beam-spectroscopy this is a nice example of how CL is implemented in a commercial FDTD package. 3) https://iopscience.iop.org/article/10.1088/2632-959X/abdc3e/meta another potentially useful reference.

For EELS 1) https://pubs.acs.org/doi/abs/10.1021/ph500408e this is a useful reference on how EELS can be implemented on FDTD

For general information regarding CL and EELS this is an excellent reference https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.82.209

Thank you Nikolaos

edit: some possibly useful links https://meep.readthedocs.io/en/latest/FAQ/#how-do-i-model-a-moving-point-charge and https://github.com/NanoComp/meep/issues/1118

[Nikolaos-Matthaiakakis]

So I understand that by using

def move_source(sim):
    sim.change_sources([mp.Source(mp.ContinuousSource(frequency=1e-10),
                                  component=mp.Ex,
                                  center=mp.Vector3(-0.5*sx+dpml+v*sim.meep_time()))])

when using existing monitors this error RuntimeError: meep: allocate field components (by adding sources) before adding dft objects cannot be avoided as meep calculates the workload for the cpu cores based on the defined source etc, so the sources need to be defined from the start.

Instead I tried capturing the fields in the time domain and used fft to obtain the frequency domain response. While this works it might be more useful or practical to have an electron beam source that still supports all the monitors available in meep so that sources can be used interchangeably. Maybe this could be done by creating a source that is made through a series of dipoles that have a certain phase delay that depends on the desired electron speed. (or could the previous error be somehow fixed?)

[stevengj]

Just add source when you create the Simulation object (it will get replaced by move_sources during the run, but that way Meep will know what field components you need).

[Nikolaos-Matthaiakakis]

Thank you,

I now get the following error

TypeError: 'Source' object is not iterable

but will keep testing

[ahoenselaar]

The Simulation initializer expects a list of sources. To pass a single source:

sim = mp.Simulation(..., sources=[src])

If this does not help, please paste a complete stack trace where the issue occurs.

[Nikolaos-Matthaiakakis]

I attach part of my test code (edited to include correction)

## Reference simulation

import meep as mp
import numpy as np
import h5py
import matplotlib.pyplot as plt
import math

sx = 60
sy = 60
cell_size = mp.Vector3(sx,sy,0)

resolution=65
courant=0.5
dpml = 1
pml_layers = [mp.PML(thickness=dpml)]

v = 0.7 # velocity of point charge

## Define wavelength range 
alpha=1.0 # corresponds to 1um 

## Monitor settings
x_span=0
y_span=sy
z_span=0
#----
x_pos=-0.5*sx+dpml*2
y_pos=10
z_pos=0
#----
frq_cen=0 # random value for testing
dfrq=500 # random value for testing
nfrq=100

# set tmp source
source_tmp=[mp.Source(mp.ContinuousSource(frequency=1e-10),
                                  component=mp.Ex,
                                  center=mp.Vector3(-0.5*sx+dpml+v*0,1,0))]
# set tmp source
sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    sources=source_tmp,
                    default_material=mp.Medium(index=1),
                    Courant=courant,
                    boundary_layers=pml_layers,
                    eps_averaging=False
                    )
dx=alpha/resolution
dt=courant*dx

## Flux monitor
box_x1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(x_pos,y_pos,z_pos),size=mp.Vector3(x_span,y_span,z_span)))

## move source function
def move_source(sim,counter):   
    sim.change_sources([mp.Source(mp.ContinuousSource(frequency=1e-10),
                                  component=mp.Ex,
                                  center=mp.Vector3(-0.5*sx+dpml+v*sim.meep_time(),1,0))])

## Time domain monitor   
field_vals_y = []
def get_slice(sim):
    field_vals_y.append(sim.get_array(size=mp.Vector3(x_span,y_span,z_span), center=mp.Vector3(x_pos,y_pos,z_pos), component=mp.Ey))

and the error

-----------
Initializing structure...
time for choose_chunkdivision = 0.000635147 s
Working in 2D dimensions.
Computational cell is 60 x 60 x 0 with resolution 65
time for set_epsilon = 5.34708 s
-----------

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-4-27cecf080310> in <module>
      1 ## Run simulation
----> 2 sim.run(move_source,get_slice,
      3     until = sx/v)

~/m_python/envs/pmp2/lib/python3.9/site-packages/meep/simulation.py in run(self, *step_funcs, **kwargs)
   3634 
   3635         if self.fields is None:
-> 3636             self.init_sim()
   3637 
   3638         self._evaluate_dft_objects()

~/m_python/envs/pmp2/lib/python3.9/site-packages/meep/simulation.py in init_sim(self)
   1935             self.fields.use_bloch(py_v3_to_vec(self.dimensions, v, self.is_cylindrical))
   1936 
-> 1937         self.add_sources()
   1938 
   1939         for hook in self.init_sim_hooks:

~/m_python/envs/pmp2/lib/python3.9/site-packages/meep/simulation.py in add_sources(self)
   2390         if self.fields is None:
   2391             self.init_sim() # in case only some processes have IndexedSources
-> 2392         for s in self.sources:
   2393             self.add_source(s)
   2394         self.fields.require_source_components() # needed by IndexedSource objects

TypeError: 'Source' object is not iterable

[ahoenselaar]

# set tmp source
sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    sources=source_tmp,
                    default_material=mp.Medium(index=1),
                    Courant=courant,
                    boundary_layers=pml_layers,
                    eps_averaging=False
                    )

where the sources= keyword argument expects a list of sources but the Source object source_tmp is passed in instead. Please try sources=[source_tmp].

[Nikolaos-Matthaiakakis]

Thank you, that was a silly oversight from me. It now works well. I will try to get a proper example running and will upload it here in case someone wants to run cathodoluminescence simulations.

[Nikolaos-Matthaiakakis]

I include a simple example of cathodoluminescence (CL) simulation in meep. I have not tested for accuracy in detail but from a quick test it agrees reasonably well with a BEM code that I often use for this. Either way this code could be a good starting point for people interested in performing CL simulations in meep. I might update this comment for future improvements.

## Example of cathodoluminescence simulation of an infinitely long gold nanowire (x=40nm,y=250nm,z=inf) 
#  - by Nikolaos Matthaiakakis 31/10/21

import meep as mp
import numpy as np
import h5py
import matplotlib.pyplot as plt
import math

## Define cell 
alpha=1.0 # corresponds to 1um 
sx = 5
sy = 5
cell_size = mp.Vector3(sx,sy)
dpml = 1
pml_layers = [mp.PML(thickness=dpml)]
resolution=350

## Calculate dx and dt
courant=0.5
dx=alpha/resolution
dt=courant*dx

## Define geometry
geom_x=0.04
geom_y=0.25
# Define material
from meep.materials import Au
material_set=Au

## Monitor settings
#----span
x_span=0
y_span=sy-2*dpml
z_span=0
#----position
x_pos=-0.5*sx+dpml+0.2
y_pos=0
z_pos=0
#----wavelength
wvl_max=1.8 # random value for testing
wvl_min=0.05 # random value for testing
nfrq=200 # number of frequency points
#---- Convert monitor frequency
frq_min = 1/wvl_max
frq_max = 1/wvl_min
frq_cen = 0.5*(frq_min+frq_max)
dfrq = frq_max-frq_min

## Create moving charge source
v = 0.7 # velocity of point charge
source_y=geom_y/2+dx # source position in y axis
def move_source(sim): # Create moving dipole source function
      sim.change_sources([mp.Source(mp.ContinuousSource(frequency=1e-10),
                                    component=mp.Ex,
                                    center=mp.Vector3(-0.5*sx+dpml+v*sim.meep_time(),source_y,0))])
# create a temporary source for simulation initialization
source_tmp=[mp.Source(mp.ContinuousSource(frequency=1e-10),
                                  component=mp.Ex,
                                  center=mp.Vector3(-0.5*sx+dpml+dx+v*0,source_y,0))]

## Reference simulation
# A reference simulation is included to capture the far fields produced by the source in free space
# due to the initialization and discritized movement of the source, and to remove them from the results of
# the main simulation
# create simulation
sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    sources=source_tmp,
                    default_material=mp.Medium(index=1),
                    Courant=courant,
                    boundary_layers=pml_layers,
                    eps_averaging=False
                    )

## Add monitors
box_x1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(x_pos,y_pos,z_pos),size=mp.Vector3(x_span,y_span,z_span))) # Add flux monitor  
field_vals_y = []
field_vals_hz = []
def get_slice(sim):# Define time domain monitor
    field_vals_y.append(sim.get_array(center=mp.Vector3(x=x_pos,y=y_pos),size=mp.Vector3(y=y_span), component=mp.Ey))
    field_vals_hz.append(sim.get_array(size=mp.Vector3(sx-2*dpml,sy-2*dpml,0), center=mp.Vector3(0,0), component=mp.Hz))

## Run reference simulation 
sim.run(move_source,mp.at_every(50*dt,get_slice),
    until = sx/v)

## Flux data reference
field_map_ref=np.array(field_vals_hz)
box_x1_data = sim.get_flux_data(box_x1)
freqs = np.asarray(mp.get_flux_freqs(box_x1))

## Main simulation

# Reset simulation
sim.reset_meep()
# reset values from time domain monitor
field_vals_y = []
field_vals_hz = []

## Create geometry
geometry = [mp.Block(mp.Vector3(geom_x,geom_y,mp.inf),
                     center=mp.Vector3(0,0,0),
                     material=material_set)] 

## create simulation including geometry
sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    sources=source_tmp,
                    default_material=mp.Medium(index=1),
                    Courant=courant,
                    boundary_layers=pml_layers,
                    geometry=geometry,
                    eps_averaging=False
                    )

## Add flux monitor
box_x1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(x_pos,y_pos,z_pos),size=mp.Vector3(x_span,y_span,z_span)))
sim.load_minus_flux_data(box_x1, box_x1_data)# Remove reference data

## Run main simulation
sim.run(move_source,
    mp.at_every(50*dt,get_slice),
    until = sx/v)

## Simulation region preview
figure1=plt.figure(dpi=100)
volume=mp.Block(size=mp.Vector3(sx,sy,0), center=mp.Vector3())
sim.plot2D(output_plane=volume)
# plot and save simulation region
figure1.tight_layout(pad=3.0)
plt.savefig("simulation.png")

# Obtain flux data
box_x1_flux = np.asarray(mp.get_fluxes(box_x1))
# Plot flux data
plt.plot(1/freqs,np.abs(box_x1_flux))
plt.xlabel('wavelength (um)')
plt.ylabel('Intensity (au)')
plt.savefig("spectra.png")

Simulation region

Time domain Hz component (reference data removed)

CL spectra