import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
from tidy3d import web
TIDY3D
LEARNING CENTER
Example Library
Two-photon absorption (TPA) and free-carrier absorption (FCA) in a Si waveguide¶
waveguide
In this notebook, we will reproduce part of the experimental results in Rong, H., Liu, A., Jones, R. et al. An all-silicon Raman laser. Nature 433, 292–294 (2005) DOI: 10.1038/nature03273. Specifically, the authors provided experimental results for the input-output power relationship for a 4.8cm long rib waveguide, with the output power following a non-linear dependence on the input power due to two-photon absorption (TPA) and free-carrier absorption (FCA) effects. All of the physical parameters are provided in the paper which makes for a great benchmark to reproduce with the first-principle TPA and FCA support in the Tidy3D solver.
First, define some fundamental parameters as they appear in the paper. It also provides the nonlinear parameters like the TPA coefficient and the FCA optical cross-section. The carrier lifetime \(\tau\) on the other hand is controlled by placing the waveguide in a p-i-n junction and controlling the voltage across the junction, which can help sweep out free carriers faster. We are going to model three different carrier lifetimes, as presented in Fig. 2 of the paper. The paper also reports a linear waveguide loss of 0.35 dB/cm. We will account for it by adding an imaginary part of the permittivity, or equivalently, a finite conductivity to the silicon medium.
wavelength = 1.55 # wavelength of interest
freq0 = td.C_0 / wavelength # corresponding frequency
wg_w = 1.5 # rib waveguide width
wg_h = 0.7 # rib waveguide height
si_h = 0.85 # silicon slab thickness
wg_l = 4.8e4 # wavelength length (4.8 cm)
si_n0 = 3.47 # refractive index of silicon
si_cond = 7.4e-5 # conductivity in S / um of silicon to account for the linear loss
tpa_beta = 0.5e-5 # TPA coefficient in um / W
fca_sigma = 1.47e-9 # FCA cross-section in um**(-2)The TPA + FCA modeling in Tidy3D follows these equations:
\[\begin{align} P_{NL} &= P_{TPA} + P_{FCA} + P_{FCPD}, \\ P_{TPA} &= -\frac{4}{3}\frac{c_0^2 \varepsilon_0^2 n_0^2 \beta}{2 i \omega} |E|^2 E, \\ P_{FCA} &= -\frac{c_0 \varepsilon_0 n_0 \sigma N_f}{i \omega} E, \\ \frac{dN_f}{dt} &= \frac{8}{3}\frac{c_0^2 \varepsilon_0^2 n_0^2 \beta}{8 q_e \hbar \omega} |E|^4 - \frac{N_f}{\tau}, \\ N_e &= N_h = N_f, \\ P_{FCPD} &= \varepsilon_0 2 n_0 \Delta n (N_f) E, \\ \Delta n (N_f) &= (c_e N_e^{e_e} + c_h N_h^{e_h}), \end{align}\]
where \(P\) denote various polarization densities, and \(N_f\) is the free carrier concentration. The factors of \(\frac{4}{3}\) and \(\frac{8}{3}\) in the equations are coming from using real fields instead of complex fields in the solver. More details can be found on this study.
We will be looking at continuous-wave operation, in which we evolve the system under a continuous-wave source until a steady state is reached both for the optical oscillations and the free-carrier distribution. In that regime, \(\frac{dN_f}{dt} = 0\), and
\[N_f = {\tau}\frac{8}{3}\frac{c_0^2 \varepsilon_0^2 n_0^2 \beta}{8 q_e \hbar \omega} |E|^4\]
This means that the free-carrier polarization density \(P_{FCA}\) is proportional to \(\sigma \tau\), which allows us to make a simplifying adjustment to be able to run the FDTD simulation within a reasonable number of time steps. Namely, the typical \(\tau\) in that system is on the order of nanoseconds, which is much larger than the optical period, and too long for the FDTD simulation. However, we can simultaneously decrease \(\tau\) by a factor and increase \(\sigma\) by the same factor, such that the steady-state concentration remains the same.
For this illustration, we also will not be simulating the whole 4.8cm waveguide from the paper. We can scale the TPA coefficient \(\beta\), making the overall nonlinearity stronger, while correspondingly shortening the waveguide. This is a valid modification to the setup since both the two-photon absorption, and the free-carrier concentration \(N_f\), and thus the free-carrier absorption, are proportional to \(\beta\).
# make the tau shorter by increasing sigma
tau_scale_factor = 1e-6 # characteristic tau will be femtoseconds instead of nanoseconds
# make the waveguide shorter by increasing beta
beta_scale_factor = 1e3Now we can define a function to generate the simulation. Note that we will not be including the doped silicon regions and the electrodes, as the authors say that they have verified that their introduction does not affect the light propagation in the linear regime.
def make_sim(
tpa: bool = True,
input_power: float = 1.0,
fca_tau: float = 1.2e-9, # seconds
):
# reduce the total waveguide length and correspondingly increase the TPA strength
wg_l_scaled = wg_l / beta_scale_factor
if tpa:
# define TPA specification
# reduce the FCA tau and correspondingly increase the cross-section
tpa_material = td.TwoPhotonAbsorption(
beta=tpa_beta * beta_scale_factor,
tau=fca_tau * tau_scale_factor,
sigma=fca_sigma / tau_scale_factor,
)
# define silicon medium
si = td.Medium(
permittivity=si_n0**2,
conductivity=si_cond,
nonlinear_spec=td.NonlinearSpec(models=[tpa_material], num_iters=5),
)
else:
si = td.Medium(
permittivity=si_n0**2,
conductivity=si_cond,
)
mat_sio2 = td.Medium(permittivity=1.44**2) # define oxide medium
# define the rib waveguide
rib = td.Structure(
geometry=td.Box.from_bounds(
rmin=[-wg_l_scaled - 10, -wg_w / 2, 0], rmax=[wg_l_scaled + 10, wg_w / 2, wg_h]
),
medium=si,
)
# define the silicon slab
slab = td.Structure(
geometry=td.Box.from_bounds(
rmin=[-wg_l_scaled - 10, -100, -si_h], rmax=[wg_l_scaled + 10, 100, 0]
),
medium=si,
)
# define the BOX layer
box = td.Structure(
geometry=td.Box.from_bounds(
rmin=[-wg_l_scaled - 10, -100, -100], rmax=[wg_l_scaled + 10, 100, -si_h]
),
medium=mat_sio2,
)
# define a mode source
source = td.ModeSource(
center=(-wg_l_scaled / 2, 0, 0),
size=(0, 4, 4),
mode_spec=td.ModeSpec(),
source_time=td.ContinuousWave(
freq0=freq0, fwidth=freq0 / 20, amplitude=np.sqrt(input_power)
), # use continuous wave
direction="+",
)
# define a flux time monitor
flux_time = td.FluxTimeMonitor(
center=(wg_l_scaled / 2, 0, 0),
size=(0, 4, 4),
name="flux_time",
interval=1,
)
# define a auxiliary field time monitor to check the auxiliary field which is related to the nonlinear polarizations
aux_field = td.AuxFieldTimeMonitor(
center=flux_time.center, size=(0, 0, 0), fields=["Nfy"], interval=2, name="aux_field"
)
# define a simulation
sim = td.Simulation(
size=(wg_l_scaled + 2, 5, 5),
structures=[box, slab, rib],
monitors=[flux_time, aux_field],
sources=[source],
run_time=2e-12,
normalize_index=None,
)
return simWe can examine the simulation and the time signal of the CW source. Note: the simulation axes do not have the physical aspect ratio as the simulation is very long along x.
sim = make_sim(tpa=False, input_power=0.5)
fig, ax = plt.subplots(3, 1, figsize=(6, 10), tight_layout=True)
sim.plot(y=0, ax=ax[0])
ax[0].set_aspect("auto")
sim.plot(z=0.2, ax=ax[1])
ax[1].set_aspect("auto")
sim.sources[0].source_time.plot(times=sim.tmesh, ax=ax[2])
plt.show()07:13:16 UTC WARNING: Monitor: 'aux_field' (simulation.monitors[1]) stores field 'Nfy', which is not used by any of the nonlinear models present in the mediums in the simulation. The resulting data will be zero.
We will first examine the waveguide mode to make sure we have what we expect.
ms = td.ModeSimulation(
center=sim.center,
size=sim.size,
structures=sim.structures,
plane=sim.sources[0],
freqs=[freq0],
mode_spec=td.ModeSpec(),
)
mode_data = web.run(ms, task_name="wg_modes")
fig, ax = plt.subplots(1, 1)
mode_data.plot_field(field_name="Ey", ax=ax)
mode_data.modes.to_dataframe()07:13:17 UTC Created task 'wg_modes' with resource_id 'mos-7d4dcf4b-0935-490e-bed1-e250996e8271' and task_type 'MODE'.
07:13:19 UTC Estimated FlexCredit cost: 0.004. For this solver type, the estimate is the final billed cost.
07:13:20 UTC status = queued
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.
07:13:30 UTC starting up solver
running solver
07:13:31 UTC status = success
07:13:34 UTC Loading results from simulation_data.hdf5
TERMS & CONDITIONS
How the process works:
Try this real URL to see an example: