Tutorial 5: Part 2 -- silicon slab transmission¶
In this part, we move on to treat more complex material dispersion. As an example, we illustrate how to use our dispersion fitting tool to fit the dispersion relation of crystalline silicon based on tabulated data, and subsequently convert it to Tidy3D components for use in computing silicon slab transmission.
# standard python imports
import numpy as np
import matplotlib.pyplot as plt
import tidy3d as td
from tidy3d import web
# fitter plugins
from tidy3d.plugins.dispersion import FastDispersionFitter as Fitter
Dispersion fitter¶
We load the dispersion data for crystalline silicon at room temperature directly from a URL. Considering that we are interested in the solar cell applications, we narrow the wavelength range to [0.4 um,1.2 um] to perform the fitting.
fitter = Fitter.from_url('https://refractiveindex.info/data_csv.php?datafile=database/data-nk/main/Si/Green-2008.yml')
fitter = fitter.copy(update={'wvl_range' : (0.4, 1.2)}) #um
Let's visulize the complete dispersion data.
plt.rcParams.update({'font.size': 18})
fig, ax = plt.subplots(1,figsize=(5,4))
ax.scatter(fitter.wvl_um, fitter.n_data, label='n', color='crimson', edgecolors='black', linewidth=0.5)
ax.scatter(fitter.wvl_um, fitter.k_data, label='k', color='blueviolet', edgecolors='black', linewidth=0.5)
ax.set_xlabel('Wavelength ($\mu m$)')
ax.legend()
plt.show()
The fitting is then performed. More details on setting up the fitter can be found here.
medium, rms_error = fitter.fit(max_num_poles=2, tolerance_rms=5e-2)
medium = medium.copy(update={'name':"fitted silicon"})
Output()
Let’s visualize how well the fit captures the dispersion data in the wavelength range of interest. As you can see, this fit looks great and should be sufficient for our simulation.
fig, ax = plt.subplots(1,figsize=(6,5))
fitter.plot(medium,ax=ax)
plt.show()
Other basic setups¶
Now let us setup the rest of the basic components.
# Wavelength and frequency range
lambda_range = (0.4, 1.2)
lam0 = np.sum(lambda_range)/2
freq_range = (td.constants.C_0/lambda_range[1], td.constants.C_0/lambda_range[0])
# frequencies and wavelengths of monitor
Nfreq = 333
monitor_freqs = np.linspace(freq_range[0], freq_range[1], Nfreq)
monitor_lambdas = td.constants.C_0 / monitor_freqs
# central frequency, frequency pulse width and total running time
freq0 = monitor_freqs[Nfreq // 2]
freqw = 0.3 * (freq_range[1] - freq_range[0])
t_stop = 100 / freq0
# Thicknesses of slabs
t_slabs = [1.0] # um
# medium of slabs
mat_slabs = [medium]
# space between slabs and sources and PML
spacing = 1 * lambda_range[-1]
# simulation size
sim_size = Lx, Ly, Lz = (1, 1, 4*spacing + sum(t_slabs))
Create Simulation¶
Now we set everything else up (structures, sources, monitors, simulation) to run the example.
# First, we define the multilayer stack structure (a single layer in this example).
slabs = []
slab_position = -Lz/2 + 2*spacing
count_slab = 0
for t, mat in zip(t_slabs, mat_slabs):
slab = td.Structure(
geometry=td.Box(
center=(0, 0, slab_position + t/2),
size=(td.inf, td.inf, t),
),
medium=mat,
name='slab'+ str(count_slab),
)
slabs.append(slab)
slab_position += t
count_slab += 1
# Here we define the planewave source, placed just in advance (towards negative z) of the slab
source = td.PlaneWave(
source_time = td.GaussianPulse(
freq0=freq0,
fwidth=freqw
),
size=(td.inf, td.inf, 0),
center=(0, 0, -Lz/2+spacing*1.5),
direction='+',
pol_angle=0,
name='planewave',
)
# We are interested in measuring the transmitted flux, so we set it to be an oversized plane.
monitor = td.FluxMonitor(
center = (0, 0, Lz/2 - spacing),
size = (td.inf, td.inf, 0),
freqs = monitor_freqs,
name='trans',
)
Now it is time to define the simulation object.
sim = td.Simulation(
center = (0, 0, 0),
size = sim_size,
grid_spec=td.GridSpec.auto(min_steps_per_wvl=30),
structures = slabs,
sources = [source],
monitors = [monitor],
run_time = t_stop,
boundary_spec = td.BoundarySpec.pml(z=True),
)
Plot The Structure¶
plt.rcParams.update({'font.size': 12})
fig, ax = plt.subplots(1, tight_layout=True, figsize=(4, 6))
sim.plot(y=0, freq=freq0, ax=ax);
ax.set_title('')
ax.set_xlabel('x ($\mu$m)')
ax.set_ylabel('z ($\mu$m)')
plt.show()
Run the simulation¶
We will submit the simulation to run as a new project.
sim_data = web.run(sim, task_name='lecture05_silicon', path='data/data_silicon.hdf5')
17:24:06 PST Created task 'lecture05_silicon' with task_id 'fdve-2e9adbd0-f29f-4e73-869d-e07360b8963e' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-2e9adbd0-f29 f-4e73-869d-e07360b8963e'.
Output()
17:24:07 PST status = queued
Output()
17:24:10 PST status = preprocess
17:24:14 PST Maximum FlexCredit cost: 0.156. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
17:24:15 PST running solver
To cancel the simulation, use 'web.abort(task_id)' or 'web.delete(task_id)' or abort/delete the task in the web UI. Terminating the Python script will not stop the job running on the cloud.
Output()
17:24:38 PST early shutoff detected at 76%, exiting.
status = postprocess
Output()
17:24:41 PST status = success
View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-2e9adbd0-f29 f-4e73-869d-e07360b8963e'.
Output()
17:24:42 PST loading simulation from data/data_silicon.hdf5
Postprocess and Plot¶
Once the simulation has completed, we can download the results and load them into the simulation object.
Now, we compare the transmitted flux computed from a transfer matrix method (TMM) code to that computed from the FDTD simulation.
# import TMM package
import tmm
# prepare list of thicknesses including air boundaries
d_list = [np.inf] + t_slabs + [np.inf]
# loop through wavelength and record TMM computed transmission
transmission_tmm = []
for i, lam in enumerate(monitor_lambdas):
# create list of refractive index at this wavelength including outer material (air)
n_list = [1, np.sqrt(medium.eps_model(td.C_0/lam)), 1]
# get transmission at normal incidence
T = tmm.coh_tmm('s', n_list, d_list, 0, lam)['T']
transmission_tmm.append(T)
transmission = sim_data['trans'].flux
fig, ax = plt.subplots(1,figsize=(5,4))
ax.plot(monitor_lambdas, transmission_tmm, '-', label='TMM')
ax.plot(monitor_lambdas, transmission, 'k--', label='Tidy3D')
ax.set_xlabel('Wavelength ($\mu m$)')
ax.set_ylabel('Transmission')
ax.set_xlim(fitter.wvl_range)
ax.legend()
plt.show()