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.

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 = \frac{1}{\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$.

Now 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.

We 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.

We will first examine the waveguide mode to make sure we have what we expect.

Indeed the mode looks like the fundamental waveguide mode, and the mode area of ~1.6 μm$^2$ matches the value given in the paper. Since we shortened the waveguide by a factor of 1000, we need to amplify the linear loss by the same factor to account for it. Here the loss of ~350 dB/cm is what we expect (0.35 dB/cm x 1000).

Now we will run FDTD simulations with three different values for the carrier lifetime, over an array of input powers.

We can examine the flux-time dependence of the longest-lifetime simulation at the highest input power.

Similarly we can plot the auxiliary field profile, which is related to the nonlinear polarizations, and see if it has reached steady state. Here we plot Nfy since the dominant field is in the $y$ direction. Different components don't interact nonlinearly as a modeling approximation

The output power is just the time-average of the time-dependent flux once steady-state is reached.

Here we turn off the warnings by setting the logging level to "ERROR". This is to prevent many repeated warnings generated from the parameter sweep. In this case, we will get warnings about the field decay not meeting the shutoff condition. Since we inject a CW source, the field will never decay so this is expected. In general, however, we highly advise against turning the warnings off since they can indicate real issues regarding the simulation setup.

Now we can compute this output power for all values of $\tau$ and the input power, and recreate Fig. 2 from the paper. The simulation results here closely match the experimental measurements.