1. Introduction

Microresonators are now ubiquitous in many fields of integrated photonics, ranging from optical communications [1], bio-sensing [2], spectroscopy [3], frequency metrology [4] and quantum optics [5]. Silicon microresonators have been historically appealing due to their compact foot-print, large third order nonlinearities and wide bandwidth of operation [6]. A bottleneck when they are operated in the C-band is the presence of Two Photon Absorption (TPA) and Free Carrier Absorption (FCA), which introduce additional losses to the device [7]. Irrespective of the application, these have often a detrimental role; examples are the limitation of the maximum attainable Signal to Noise Ratio (SNR) in optical links [8], gain saturation in wavelength conversion devices [9] and reduction of the heralding efficiency in parametric photon pair sources [10]. Nonetheless, there are some exceptions in which TPA and FCA in resonators are sought for all optical information processing, such as memory storage [11], logic operations [12] or unidirectional propagation [13]. TPA creates electron-hole pairs whose relaxation self-heats the resonator, so that the temperature and the carrier dynamics become coupled [14]. Their interplay may appear at the output of the resonator as time-instabilities associated to Self-Pulsing (SP), bistability, excitability, or as a nonlinear distortion in the transmitted spectrum [8,15–18]. There is a growing interest in exploiting these regimes for neuromorphic photonics, making silicon resonators promising candidates for the realization of spiking neural networks, single perceptrons, or extended reservoirs [19–23]. The design and optimization of these systems require the accurate modeling of each node, which is traditionally implemented by a set of three coupled ordinary differential equations (ODE): one for the internal energy of the resonator, one for the temperature and one for the free carrier population [24]. The linearization of this set around equilibrium points, combined with the tool of Hopf bifurcations, allows to identify the optimal regions of operation in the parameter space defined by the input power $P$ and the resonance detuning $\Delta \nu$ of the pump laser [17,23,24]. In this work, we experimentally and theoretically show that there are cases in which the commonly used ODEs fail to describe the dynamics of the silicon microresonator. This happens whenever the Newton’s law of cooling, from which the temperature equation derives, does not approximate the actual thermal dynamics. Using Finite Element Method (FEM) simulations, we show that in our planar geometry of a silicon racetrack resonator embedded in a silicon dioxide cladding, the temperature is described by a fast ($\sim 100\,\textrm{ns}$) and a slow ($\sim 1\,\mu \textrm{s}$) relaxation time. The latter can be observed when self oscillations of similar period are internally established due the interplay between thermal and free carrier nonlinearities. In our case, a SP regime of sub-MHz frequency establishes as a consequence of a relatively long ($\sim 40\,\textrm{ns}$) free carrier lifetime. We introduce a new temperature equation that reproduces such temporal features, leading to predictions very close to the experimental observations. We investigate the generalization capabilities of our model to other resonators with similar geometry but different quality factors, proposing and experimentally validating a simple method to assess the error on the predicted SP period. Finally, we discuss the conditions where Newton’s law can be applied, and the artifacts which arise on the specific heat, on the carrier lifetime or on the carrier dispersion coefficient when an inadequate model for temperature is used to fit the experimental data.

2. Observation of sub-MHz self pulsing

Our racetrack microresonator (MR) has been fabricated at the IMEC/Europractice facility within a multi-project wafer run on a $200\,\textrm{mm}$ SOI wafer of $220\,\textrm{nm}$ thickness. Silicon waveguides are embedded in a silica cladding, and have a width of $450\,\textrm{nm}$ to ensure single mode operation on both polarizations. An Add-Drop filter configuration is realized by evanescently coupling two bus waveguides to the MR of $7\,\mu \textrm{m}$ of radius. This is accomplished through a $250\,\textrm{nm}$ gap and a coupling length of $3\,\mu \textrm{m}$, which yields a measured coupling coefficient of $\kappa ^{2}=0.063(3)$. From the low-power transmission spectra recorded at the Through port, shown in Fig. 1(a), we extracted an intrinsic quality factor (Q) of $Q_{\textrm{i}} = 1.11(8)\times 10^5$ and a loaded Q of $Q_{\textrm{L}} = 6.5(2)\times 10^{3}$. We used the experimental setup shown in Fig. 1(b) to probe the resonator dynamics. A Continuous Wave (CW) tunable laser at $1555\,\textrm{nm}$ is amplified and injected at the input waveguide of the MR through a grating coupler ($\sim 3.5\,\textrm{dB}$ insertion loss). The input power is controlled by a Variable Optical Attenuator. Optionally, a C-band pulsed laser source of $40\,\textrm{ps}$ pulse width and $1\,\textrm{MHz}$ of repetition rate can be combined to the CW laser to perform pump and probe experiments. Light transmitted from the Drop port of the resonator is coupled off-chip by a second grating coupler, subsequently amplified and directed to a fast photodiode with $20\,\textrm{GHz}$ bandwidth. A tunable band-pass filter centered at the laser wavelength is used to remove the spontaneous emission of the amplifiers, hence increasing the SNR. A $40\,\textrm{GS}$ oscilloscope records the output of the photodiode. As a first experiment, we investigated the stability regions of the MR in the two-dimensional plane $(P,\Delta \nu )$. The map shown in Fig. 2(a) marks with circles the points where SP is observed, and with crosses the ones where the system has a stable CW output. This experiment is performed by resetting the initial conditions at each point $(P,\Delta \nu )$ by closing the input VOA shutter for few seconds. Self pulsing is observed at both positive and negative detunings, predominantly on the blue side of the resonance. Despite the similar MR geometry, this contrasts with the results in [16], where SP is exclusively detected at the blue side of the resonance. According to the linear stability analysis in [24], the blue shift of the SP region indicates that free carrier dispersion (FCD) is dominating over the thermo optic effect (TOE), which will find justification in Section 3. and Section 5. The map shown in Fig. 2(b) reports the frequency of SP, which we extract from time traces similar to those shown in Fig. 2(c). The general trend agrees with the results of [14,16], where it is reported that the SP frequency increases from positive to negative detunings and by increasing the input power. The time traces in Fig. 2(c) show all the distinctive features of SP induced by the interplay between FCD and TOE: the cycle starts with a free carrier bistability (narrow peak in the drop signal), which is then counteracted by a slower thermal shift pulling the resonance frequency in the opposite direction, partially restoring the initial transmittivity. When the thermal shift overwhelms FCD, the carrier concentration suddenly drops, leaving an hot cavity completely out of resonance (drop of the transmittivity in Fig. 2(c)). This is followed by a slow thermal relaxation, where the temperature of the MR decays until the initial conditions are restored and a new cycle of oscillation begins. As shown in Fig. 2(b), the average SP frequency is of the order of $\sim \,1\,\textrm{MHz}$, and can be as low as $\sim \,400\,\textrm{kHz}$ for large positive detunings. This is an order of magnitude lower than the SP frequencies reported in [16,25], despite the very similar MR geometries. In general, SP periods of $\sim 100\,\textrm{ns}$ arise when the thermal and the carrier lifetime are respectively of the order of $\tau _{\textrm{th}} \sim 80-150\,\textrm{ns}$ and $\tau _{\textrm{fc}}\sim 1-10\,\textrm{ns}$ [24]. This is because a complete cycle of oscillation requires both thermal and carrier relaxation, so its period is roughly given by $\tau _{\textrm{fc}}+\tau _{\textrm{th}}$. In order to justify a sub-MHz dynamics, a longer thermal and/or carrier lifetime has to be postulated. Since the first can be assessed from a FEM simulation, we performed a second experiment to measure the carrier lifetime.

 

Fig. 1. (a) Low power transmission spectra of the MR, collected at the Through port. Highlighted in red is the resonance order where the SP regime is analyzed. An enlarged view of this resonance is shown on the right panel. Here, the blue region marks the range of laser detuning $\Delta \nu$ investigated in this work. (b) The experimental setup implemented to measure the SP of the MR and the carrier lifetime. TLS = Tunable Laser Source, PLS = Pulsed Laser Source, FPC = Fiber Polarization Controller, EDFA = Erbium Doped Optical Amplifier, VOA = Variable Optical Attenuator, BS = Beam Splitter, DUT = Device Under Test, BPF = Band Pass Filter, PD = Photodiode, OSC = Oscilloscope.

Download Full Size | PDF

 

Fig. 2. (a) Stability map of the MR in the $\Delta \nu$ = laser frequency detuning, P = input power plane. Red crosses indicate the points where the MR, after an initial transient, shows stable output. Green circles indicate the points where the MR is SP. (b) Map of the oscillation frequency of the self pulsing regime. (c) Examples of time traces recorded at the output of the Drop port of the MR. The maximum of the intensity is normalized to one. The labels (1), (2) and (3) refer to the values of $(P,\Delta \nu )$ associated to each trace, which are indicated in panel (b).

Download Full Size | PDF

3. Determination of the free carrier lifetime

The carrier lifetime is determined with a pump and probe technique applied to a spiral waveguide of length $7.1\,\textrm{mm}$, located on the same die of the MR. As shown in Fig. 1(b), a strong pulsed pump ($40\,\textrm{ps}$) with $\sim 1.5\,\textrm{W}$ peak power and a wavelength of $1550\,\textrm{nm}$ is combined with the CW probe and injected into the spiral. The high intensity of the pump creates an initial free carrier population $\Delta N_0$ through TPA, which recombines via trap states located at the Si/SiO$_2$ interface [26]. During this transient regime, the CW probe experiences an extra-loss due to FCA, whose temporal profile follows the instantaneous free carrier population $\Delta N(t)$ along the waveguide. The time-dependent transmittivity of the spiral (see inset in Fig. 3(a)) is used to extract the instantaneous carrier lifetime $\tau _{\textrm{fc}}$, using the approach detailed in [27]. Figure 3(b) shows $\tau _{\textrm{fc}}$ as a function of the carrier concentration $\Delta N$ inside the waveguide.

The value of $\tau _{\textrm{fc}}$ is bounded between $20\,\textrm{ns}$ and $60\,\textrm{ns}$. As shown in Fig. 3(b), at the very early stage, carriers recombine with a characteristic lifetime of $\tau _{\textrm{slow}}\sim 60\,\textrm{ns}$, which is significantly larger than the one in the late stage ($\tau _{\textrm{fast}}\sim 25\, \textrm{ns}$), in which the carrier concentration dropped by nearly two orders of magnitude. Even though the instantaneous lifetime is larger than the ones traditionally reported for SOI waveguides [28–30], this is not entirely surprising. Indeed, as reported in [27], the free carrier lifetime appears to evolve nonlinearly with the carrier concentration, and ultra-long lifetimes of the order of $300\,\textrm{ns}$ have been observed. It is worth to mention that the waveguide geometry and the fabrication facility adopted in [27] is the same of our work.

 

Fig. 3. (a) Temporal decay of the carrier concentration $\Delta N$ subsequent to the absorption of a pump pulse at $t=0$. The two red dashed lines indicate the expected profile for a purely exponential decay, i.e., $\log (\Delta N)=\log (\Delta N_0)-\frac {t}{\tau _{\textrm{fc}}}$. The line indicated as $\tau _{\textrm{slow}}\sim 60\,\textrm{ns}$ is obtained by evaluating the average local slope for $0\leq t\leq 50\,\textrm{ns}$, while the time interval $250\leq t\leq 300\,\textrm{ns}$ is used to evaluate the slope of the line labeled as $\tau _{\textrm{fast}}\sim 25\,\textrm{ns}$. The inset shows the time-dependent transmittivity of the CW probe beam. Dips mark the time arrival of the pump pulses, which locally induce extra-losses due to free carrier absorption. The transmission has been normalized by taking into account the linear propagation loss. (b) Measured (black dots) carrier lifetime $\tau _{\textrm{fc}}$ as a function of the carrier concentration. This is extracted starting from curves similar to the one in panel (a), and by using the method outlined in [27]. The red dashed line is a fit of the lifetime obtained by implementing the Shockley Read Hall expression in Eq. (1) and the parameters listed in Table 2 of Appendix B.

Download Full Size | PDF

Trap-assisted recombination leads to electron-hole lifetimes which depend on the amount of excess carriers [31–33]. When the density of traps is small compared to the equilibrium concentration of the majority carriers, the following Shockley Read Hall (SRH) expression holds [32]:

4. Evidences of inadequacy of the standard coupled equations

The temporal dynamics of silicon MR is commonly modeled by a set of three, coupled ordinary differential equations, which govern the internal energy $U_{\textrm{int}}=|a|^2$ of the resonator, the free carrier population $\Delta N$, and the differential temperature with the cladding $\Delta T$ [14]:

As a first step, we replaced the constant free carrier lifetime $\tau _{\textrm{fc}}$ in Eq. (3) with the expression in Eq. (1). We then numerically integrated Eqs. (2)–(4) by changing $\tau _{\textrm{th}}$ until the simulated SP period matches the one of the experiment. Arbitrarily, we used as a reference the experimental SP trace for $P = 7\,\textrm{dBm}$ and $\Delta \nu = -10\,\textrm{GHz}$. The best value of $\tau _{\textrm{th}}$ is found by repetitively minimizing the least squares error between the experimental and the simulated Drop signal by using a Particle Swarm Optimizer (PSO) [34]. The best match is found with $\tau _{\textrm{th}} = 270(10)\,\textrm{ns}$. Remarkably, we had to increase the FCD coefficient up to $\sigma _{\textrm{FCD}} = -9.5\times 10^{-27}\,\textrm{m}^3$ in order to observe SP, which is $\sim 2-3$ times higher than the values reported in literature [35]. Using the set of parameters listed Table 1 of Appendix A, we generated the stability and the SP frequency maps shown in Figs. 4(a) and (b). A comparison between some simulated and experimental traces is shown in Fig. 4(c) (labelled as "Model 1"). Clearly, the overlap (fraction of matching points) with the map in Fig. 2(a) is very poor. In contrast with the experiment, the SP region is totally confined at negative detunings, and the frequency trend in Fig. 4(b) does not follow the one in Fig. 2(b). Moreover, as emerges from Fig. 4(c), the predicted duty cycle is almost half the one which we measured. We also performed more refined analysis, in which, for each pair of $(\sigma _{\textrm{FCD}},\tau _{\textrm{th}})$, we generated the same maps in Figs. 4(a) and (b), and we evaluated the overlap with the ones in Figs. 2(a) and (b), but we never obtained overlaps higher than $60\%$. According to [24], the SP region shifts towards positive detunings if we increase the relative weight $q$ between FCD and TOE. Using the same notation of [24], this scales as $q\propto \frac {c_p}{\tau _{\textrm{th}}\left ( \frac {dn}{dT}\right )}\sqrt {\sigma _{\textrm{FCD}}\tau _{\textrm{fc}}}$, where $\left ( \frac {dn}{dT}\right )$ is the thermo-optic coefficient and $c_p$ the specific heat. Since $\tau _{\textrm{fc}}$ has been measured, and the continuity equation governing free carrier evolution is well established [36], we focused our attention to the temperature equation. We let the specific heat and the thermal relaxation time to be free parameters, and we independently fit each experimental trace in the plane $(P,\Delta \nu )$. This produced the values of $c_p(P,\Delta \nu )$ and $\tau _{\textrm{th}}(P,\Delta \nu )$ shown in Fig. 4(d), while the related temporal traces are shown in Fig. 4(c) (labelled as "Model 2"). Several considerations now arise, with the most striking being the fact that, on average, the specific heat should be $\sim 7-10$ times higher than its tabulated reference value $c_{\textrm{p,ref}}=700\,\textrm{J}/\textrm{(kgK)}$ . Furthermore, $\tau _{\textrm{th}}$ lies between $350-650\,\textrm{ns}$, and as $c_p$, shows a clear correlation with the SP frequency in Fig. 2(b). In particular, the longer the SP period, the higher is the value of $c_p$ and $\tau _{\textrm{th}}$ required to fit the data. This is shown in Fig. 4(d), which evidences a linear correlation between the variables, indicating that the $q$ parameter is constant in the plane $(P,\Delta \nu )$. Using the average values $\tau _{\textrm{th}} = 450$ and $c_p = 8c_{\textrm{p,ref}}$, we obtain $q\sim 5q_0$, where $q_0$ is the value of $q$ generating the maps in Figs. 4(a) and (b). This confirms that the original model was underestimating the role of free carriers on the overall resonance shift. The agreement with the experimental traces in Fig. 4(c) is excellent, but it is reached with the nonphysical assumption that $c_p$ should be almost one order of magnitude larger than the one reported in literature. Similarly, the thermal decay constant is $\sim 3.5$ times higher than the ones reported in [16,25], despite the MR geometry and the material platform are almost the same. This suggests that they may be artifacts which fix the inadequacy of the temperature equation.

 

Fig. 4. (a) Simulated stability map obtained by integrating Eqs. (2)–(4) and by using the set of parameters listed in Table 1 of Appendix A. (b) Simulated frequency of the SP cycles. (c) Examples of experimental (black) and simulated time traces at the Drop port of the resonator, corresponding to points (1), (2) and (3) in the $(P,\Delta \nu )$ plane shown in panel (b). Red lines use Eqs. (2)–(4) with $c_p=c_{\textrm{p,ref}} = 700\,\textrm{J/(kgK)}$ and $\tau _{\textrm{th}}=275\,\textrm{ns}$. Blue lines use Eqs. (2)–(4) and the values of $c_p$ and $\tau _{\textrm{th}}$ shown in panel (d). (d) Normalized specific heat $c_p/c_{\textrm{p,ref}}$ as a function of the thermal lifetime $\tau _{\textrm{th}}$. The color of the dots labels the position in the plane $(P,\Delta \nu )$. This color code is shown in the upper left inset.

Download Full Size | PDF

5. Temperature equation: FEM vs ODE predictions

The ODE for temperature in Eq. (4) is derived from the Netwon’s law of cooling, which has a single exponential as autonomous solution [37]. However, we can use a FEM solver (COMSOL Multiphysics [38]) to show that Newton’s law does not model the temperature dynamics of our MR. As an example, in Fig. 5(a) we report the average temperature of the MR when a heat source supplying a costant power of $2\,\textrm{mW}$ for $20\,\mu \textrm{s}$ is placed in the waveguide core. These conditions emulate what happens during SP, where the absorbed power origins from TPA and FCA.

The inset in Fig. 5(a) indicates that the temperature profile in the surroundings of the MR is not spatially uniform. Simulations reveal that the external temperature distribution is evolving in time, since heat progressively propagates toward the outer regions of the silica cladding and beneath the silicon substrate. The region which embeds the silicon waveguide does not act as a thermal bath at a constant temperature, which is one of the assumptions behind Newton’s law. As a direct consequence, as shown in Fig. 5(a), the raising(falling) edge of $\Delta T$ is not a single exponential. At this point, one could fix the discrepancy by directly coupling Eqs. (2) and (3) to the FEM solver, and integrate the temperature profile inside the MR to lump the spatial variables, i.e., $\Delta T(x,y,z,t) \rightarrow \Delta T(t)$. This is certainly possible, but on one hand it is time consuming, and on the other it is not suited for a linear stability analysis. Alternatively, one can construct an ODE which mimics the temporal solution given by the FEM. The starting point is the observation that a double exponential decay of the form $\Delta T = C_s\exp ^{-\frac {t}{\tau _s}} + C_f\exp ^{-\frac {t}{\tau _f}}$ (solid blue line in Fig. 5(a)) well matches the FEM solution. Here, we call $\tau _s\gg \tau _f$ respectively the slow and the fast thermal lifetimes, while $C_1$ and $C_2$ are constants to be determined from the initial conditions. The associated ODE is $\frac {d^2\Delta T}{dt^2}+\left ( \frac {1}{\tau _{\textrm{eff}}} \right )\frac {d\Delta T}{dt}+\left ( \frac {1}{\tau _s\tau _f} \right ) \Delta T = 0$, where $\tau _{\textrm{eff}}=\left (\frac {1}{\tau _s}+\frac {1}{\tau _f}\right )^{-1}$. This equation is homogeneous, so we must provide the driving terms related to the absorbed power. Heuristically, we propose the form:

 

Fig. 5. (a) Temperature evolution $\Delta T$ (black) of the MR when a constant power of $2\,\textrm{mW}$ is applied through the whole cavity from $0$ to $20\,\mu \textrm{s}$. After $20\,\mu \textrm{s}$, the power is switched off. The temperature is obtained by solving the heat transfer equation with a commercial FEM solver (COMSOL Multiphysics [38]), and by spatially averaging the temperature profile inside the core of the waveguide. The inset shows a cross section of the instantaneous distribution of $\Delta T$ at $t\sim 21\,\mu \textrm{s}$ (axisymmetric view). The red line is a single exponential fit of the raising edge of $\Delta T$. The blue line fits the falling edge with a double exponential. (b) Temperature variation of the MR predicted by the FEM solver (black squares) through some cycles of SP for $(P=8\,\textrm{dBm},\Delta \nu =0\,\textrm{GHz})$. The driving power $P_{\textrm{abs}}$ is obtained by solving Eqs. (2)–(4) with $c_p = 6.9c_{\textrm{p,ref}}$ and $\tau _{\textrm{th}} = 417\,\textrm{ns}$. The red dashed line (ODE1 (1)) is obtained by solving Eqs. (2)–(4) using the same driving power $P_{\textrm{abs}}$, but by forcing $c_p=c_{\textrm{p,ref}}$. The value of the thermal lifetime $\tau _{\textrm{th}}=40\,\textrm{ns}$ minimizes the least square error with the FEM trace. The solid red line (ODE1 (2)) is obtained by solving Eqs. (2)–(4) with $c_p=6.9c_{\textrm{p,ref}}$ and $\tau _{\textrm{th}} = 417\,\textrm{ns}$. The solid blue line (ODE2) is obtained by integrating Eq. (7) and by minimizing the least square error with the FEM trace. This yielded $(\kappa = 8.87\times 10^{-4}\,\textrm{ns}^{-1} ,\tau _f = 55.8\,\textrm{ns} ,\tau _s = 1200\,\textrm{ns} ,\eta = 0.31)$.

Download Full Size | PDF

In a last step, we learned the parameters $(\kappa ,\eta ,\tau _s,\tau _f)$ using the experimental time response of the MR. This is accomplished by running a multi-objective PSO optimizer simultaneously minimizing the least squares error between the MR output predicted by Eqs. (2), (3),and (7) and three arbitrarily chosen experimental traces in the plane $(P,\Delta \nu )$. The result of the optimization yield $(\kappa = 7.45(5)\times 10^{-4}\,\textrm{ns}^{-1} ,\tau _f = 71(2)\,\textrm{ns} ,\tau _s = 1030(20)\,\textrm{ns} ,\eta = 0.29(3))$, where the errors have been taken as the standard deviation over $10$ independent runs of the algorithm. We used this set to generate maps of stability, SP frequency and duty cycle, which we report in Figs. 6(a) and (b). The improvement obtained by adopting the new model Eq. (7) is substantial. The overlap between the experimental and the simulated stability map raised from $41\%$ to $89\%$; the relative error on the SP frequency decreased from $13(7)\%$ to $7(8)\%$ and the one on the duty cycle from $60\%$ to $30(10)\%$. As can be appreciated from the time traces in Fig. 6(c), all the characteristic features of SP discussed in Section 2 are maintained by the new model given by Eq. (7). Worth to note, the improved agreement between the simulation and the experiments is obtained by solely correcting the temperature equation, with all the other parameters extracted from the experiment or taken from tabulated values (including $c_p=c_{\textrm{p,ref}}$).

 

Fig. 6. (a) Simulated stability map which uses Eqs. (2) and (3) for the energy and the carrier concentration and Eq. (7) for the temperature. The values of $(\kappa ,\tau _f,\tau _s,\eta )$ are indicated in the main text. (b) Upper panel: comparison between the simulated (lower bars) and the experimental (upper bars) SP frequencies. The position in the plane $(P,\Delta \nu )$ is indicated by the same color code of Fig. 4(d). Lower panel: same representation used in the upper panel, but relative to the duty cycle. This is defined as DC = $T_{\textrm{th}}\nu _{\textrm{SP}}$, where $T_{\textrm{th}}$ is the total time in which the Drop trace lies above the threshold $0.45\cdot \textrm{(max[Drop]-min[Drop])}$, while $\nu _{\textrm{sp}}$ is the SP frequency. (c) Examples of experimental (black) and simulated (red) Drop traces relative to the points (1), (2) and (3) in panel (a).

Download Full Size | PDF

6. Accuracy of the temperature equation

Equation (7) fixes the discrepancies between the model and the experiment by introducing an inertial term proportional to $\Delta T$ to increase the relative weight of FCD with respect to TOE . However, Eq. (7) is not isomorphic to the heat transfer equation (where we implicitly assumed that the spatial coordinates have been dropped by integrating the temperature profile over the MR volume). In this sense, it is an approximation whose range of validity has to be quantified. A natural question is whether or not it is possible to train the model on a specific SP frequency range and use it to predict the SP behaviour of a MR having different oscillation period. We addressed this question by training Eq. (7) to reproduce the solution of the FEM solver when the heat source is a square wave with $33\%$ duty-cycle, $2\,\textrm{mW}$ of amplitude, $100\%$ extinction and $T_{\textrm{ref}} = 1\,\mu \textrm{s}$ of period. These conditions well approximate the temporal profile of $P_{\textrm{abs}}$ during the SP regime of a MR. We then changed the period of the square wave to some other value $T_{\textrm{sp}} \neq T_{\textrm{ref}}$, and used the FEM solver to evaluate the maximum temperature excursion $\Delta T_{\textrm{max}}$ for each SP cycle (initial transient effects are discarded in this analysis). For example, in Fig. 5(b), $\Delta T_{\textrm{max}} \sim 6\,\textrm{K}$ for $T_{\textrm{sp}} = 0.8\,\mu \textrm{s}$. Note that $\Delta T_{\textrm{max}}$ is a function of $T_{\textrm{sp}}$. We assume that the temperature variations predicted by FEM represent what occurs in a MR with a SP period of $T_{\textrm{sp}}$. Using the same driving power $P_{\textrm{abs}}$ in Eq. (7), we obtained the approximated ODE solution, from which we evaluated the time $T_{\textrm{ODE}}$ required to reach the same temperature variation $\Delta T_{\textrm{max}}$. The error on the predicted SP period by Eq. (7) is then determined by $T_{\textrm{sp}}-T_{\textrm{ODE}}$. This error is a function of the SP period and is shown as a solid red line in Fig. 7. This approach is then experimentally validated. We measured the SP period $T_{\textrm{exp}}$ of a MR with an identical perimeter as the one described in Section 1., but with $Q=20000$, and evaluated the difference with the SP period $T_{\textrm{sim}}$ predicted by the full model given by Eqs. (2), (3), and (7). To this purpose, we used the set of parameters $(\kappa ,\eta ,\tau _s,\tau _f)$ learned from the fits in Fig. 6(c), which are related to a MR with $Q=6500$. The errors on the SP period $T_{\textrm{exp}}-T_{\textrm{sim}}$ as a function of $T_{\textrm{sp}}$ are shown in Fig. 7 (black diamonds). As a comparison, we also included the errors on the SP period $T_{\textrm{exp}}-T_{\textrm{sim}}$ relative to the MR with $Q=6500$ (black squares). As can be seen, the errors of the model with respect to the experiment or to the full FEM simulation are similar, which allows us to evaluate the reliability of the approximated temperature model. In general, Eq. (7) tends to underestimate the SP period, with an error that super-linearly increases with the latter. As expected, the relative error is at minimum where the model has been trained (i.e., $1\,\mu \textrm{s}$). If we set the desired accuracy threshold to $10\%$, the model loses accuracy as long as the period exceeds $2\,\mu \textrm{s}$ or gets lower than $500\,\textrm{ns}$. This interval is sufficiently large to accommodate all the SP frequencies reported in Fig. 2(b).

 

Fig. 7. Predicted absolute error (red line) on the SP period in case Eq. (7) is trained to reproduce the temperature evolution of the MR driven by a square wave with $2\,\textrm{mW}$ amplitude, $33\%$ duty-cycle and $1\,\mu \textrm{s}$ period. The relative error is shown with a red dashed line. Black squares indicate the error between the simulated and the experimental period for a MR of $Q=6500$. Black diamonds are relative to a MR of $Q=20000$, which naturally exhibits longer SP periods. Only the experimental data points associated to a duty-cycle in the range $[20-30]\%$ are shown for a fair comparison with the simulation.

Download Full Size | PDF

7. General remarks on the different modeling approaches

In the previous section, we have seen that the model trained on a given MR does not generalize well to other MR with SP periods that are $\sim 2-3$ times higher(lower). It is then expected that different MR have to be modeled by different sets $(\kappa ,\eta ,\tau _s,\tau _f)$. This holds also for MR which share the same geometry but have different Q. This is of little use for the engineering of such dynamical systems, where we would like to know in advance their performances before fabrication and testing. An universal modeling of these systems can only be obtained if Eqs. (2) and (3) are coupled to the FEM solution of the temperature equation. In order to reduce the computational effort, a possible strategy could be to use the FEM solver to obtain accurate solutions for only very few points in the phase space $(P,\Delta \nu )$, which are then used to train Eq. (7) to predict the time response for all the other points. We expect that Newton’s law will nevertheless produce accurate solutions for SP periods of the order of $\tau _f$, since the slow thermal lifetime will appear as a simple constant background. However, as shown in Fig. 5(b), this will not fix the need of increasing $c_p$ to reproduce the correct temperature excursion. In principle, we could also leave $c_p=c_{\textrm{p,ref}}$, but we will have to artificially increase $\tau _{\textrm{fc}}$ and/or $\sigma _{\textrm{FCD}}$ in order to restore the relative weight between FCD and TOE. Similarly, a single exponential solution works well if the driving input power is changing much faster than $\tau _f$, as it happens for OOK signals transmitted at Gpbs data rates [8,39]. Indeed, these are effectively seen as CW signals from the point of view of the temperature evolution, so $\tau _{\textrm{th}}$ can be tailored to match the steady state temperature observed in the experiment or predicted by the FEM simulation.

There could be also some cases for which the single exponential solution can be safely used without any sort of corrections and for every SP period, which are the ones where Newton’s law of cooling holds. An example could be an air cladded silicon microdisk/microtoroid suspended on a silica pedestal, as the ones in [14,40]. In this case, the substrate acts as a thermal bath with constant temperature, which justifies the validity of Newton’s law.

8. Conclusions

In this work, we refined the coupled equations describing the dynamics of silicon microresonators under thermal and free carrier nonlinearities. We showed that for MR buried in a silicon dioxide cladding, the Newton’s law of cooling, which was traditionally used to model the waveguide temperature, can be inadequate. In particular, it leads to an overestimation of the thermo optic effect over free carrier dispersion. The immediate consequence, which we support by experimental evidences, is that these equations wrongly predicts the positions of the stability regions in the $(P,\Delta \nu )$ plane. Using inputs from FEM simulations, we proposed and implemented an alternative equation for the temperature evolution, which overcomes most of the limitations of Newton’s law. The new equation introduces a slow and a fast thermal relaxation time, which can be both observed when, as in our case, the latter is comparable to the SP period. This decisively improved the agreement with the experiment of both the stability maps and the shape of the self pulsing waveforms. We envisage that FEM simulations could be used as inputs for training this ODE model for temperature, which simultaneously speeds up computation and allows to handle the linear stability analysis. As a last step, we discussed those conditions where Newton’s law can still lead to accurate predictions, but at the cost of introducing artifacts on the value of the specific heat, the carrier lifetime or the free carrier dispersion coefficient. We believe that this work completes the comprehension of a well established phenomena in silicon resonators, and can be helpful for the design and performance assessment of telecommunication devices or dynamical nodes for neuromorphic computing.

Appendix A: values of the terms in the coupled equations

Table 1 lists the definitions and the values of all the terms which appear in the set of coupled Eqs. (2)–(4), with the exception of $\tau _{\textrm{fc}}$ and $\tau _{\textrm{th}}$, that are defined in the main text. The normalized resonance shift induced by free carriers is modeled as $\delta _{\textrm{fc}} = -\frac {\sigma _{\textrm{FCD}}\Delta N}{n_0}$ to produce the maps in Figs. 4(a) and (b). The refined expression $\delta _{\textrm{fc}} = -\frac {\sigma _{\textrm{FCD,1}}\Delta N+\sigma _{\textrm{FCD,2}}\Delta N^{0.8}}{n_0}$ [35] is used to produce the maps in Figs. 6(a) and (b).

Table 1. Definition, description and values of the parameters which appear in Eqs. (2)–(4) in the main text.

View Table | View all tables in this article

Appendix B: values of the terms in the Shockley Read Hall expression for $\tau _{\textrm{fc}}$

Trap assisted recombination at the Si/SiO$_2$ interface of the silicon waveguide is described using the well known Shockley Read Hall expression in Eq. (1) [32]. Here, it is assumed that the recombination is dominated by a single trap level, lying at energy $E_t$ inside the silicon energy gap $E_g$. The density of traps and the electron(hole) capture cross-section at the interface lumps into a phenomenological characteristic lifetime $\tau _n$($\tau _p$). The terms $(a,b,\tau _0)$ appearing in Eqs. (2)–(4) have the following definitions:

(8)$$a = \frac{\tau_p+\tau_n}{\tau_p(n_1+n_0)+\tau_n(p_1+p_0)} \quad\quad b = \frac{1}{n_0+p_0} \quad \quad \tau_0 = \tau_p\frac{n_0+n_1}{n_0+p_0}+\tau_n\frac{p_0+p_1}{n_0+p_0}$$

in which all the above quantities are defined in Table 2.

Table 2. Definition, description and values of the parameters which appear in Eq. (8).

View Table | View all tables in this article

Disclosures

The authors declare no conflicts of interest.

References

1. F. Testa, S. Tondini, F. Gambini, P. Velha, A. Bianchi, C. Kopp, M. Hofbauer, C. L. Manganelli, N. Zecevic, S. Faralli, G. Peres, R. Enne, A. Serrano, B. Goll, G. Fontana, A. Chalyan, J. M. Lee, P. Pintus, G. Chiaretti, H. Zimmermann, L. Pavesi, C. J. Oton, and S. Stracca, “Integrated reconfigurable silicon photonics switch matrix in iris project: Technological achievements and experimental results,” J. Lightwave Technol. 37(2), 345–355 (2019). [CrossRef]  

2. P. Steglich, M. Hülsemann, B. Dietzel, and A. Mai, “Optical biosensors based on silicon-on-insulator ring resonators: A review,” Molecules 24(3), 519 (2019). [CrossRef]  

3. M.-G. Suh, Q.-F. Yang, K. Y. Yang, X. Yi, and K. J. Vahala, “Microresonator soliton dual-comb spectroscopy,” Science 354(6312), 600–603 (2016). [CrossRef]  

4. S. B. Papp, K. Beha, P. Del’Haye, F. Quinlan, H. Lee, K. J. Vahala, and S. A. Diddams, “Microresonator frequency comb optical clock,” Optica 1(1), 10–14 (2014). [CrossRef]  

5. D. Llewellyn, Y. Ding, I. I. Faruque, S. Paesani, D. Bacco, R. Santagati, Y.-J. Qian, Y. Li, Y.-F. Xiao, M. Huber, M. Malik, G. F. Sinclair, X. Zhou, K. Rottwitt, J. L. O’Brien, J. G. Rarity, Q. Gong, L. K. Oxenlowe, J. Wang, and M. G. Thompson, “Chip-to-chip quantum teleportation and multi-photon entanglement in silicon,” Nat. Phys. 16(2), 148–153 (2020). [CrossRef]  

6. M. Borghi, C. Castellan, S. Signorini, A. Trenti, and L. Pavesi, “Nonlinear silicon photonics,” J. Opt. 19(9), 093002 (2017). [CrossRef]  

7. M. Borselli, T. J. Johnson, and O. Painter, “Beyond the rayleigh scattering limit in high-q silicon microdisks: theory and experiment,” Opt. Express 13(5), 1515–1530 (2005). [CrossRef]  

8. M. de Cea, A. H. Atabaki, and R. J. Ram, “Power handling of silicon microring modulators,” Opt. Express 27(17), 24274–24285 (2019). [CrossRef]  

9. Z. Wang, H. Liu, Q. Sun, N. Huang, S. Li, and J. Han, “The influence of thermal and free carrier dispersion effects on all-optical wavelength conversion in a silicon racetrack-shaped microring resonator,” Laser Phys. 26(7), 075403 (2016). [CrossRef]  

10. G. F. Sinclair, N. A. Tyler, D. Sahin, J. Barreto, and M. G. Thompson, “Temperature dependence of the kerr nonlinearity and two-photon absorption in a silicon waveguide at 1.55 μ m,” Phys. Rev. Appl. 11(4), 044084 (2019). [CrossRef]  

11. V. R. Almeida and M. Lipson, “Optical bistability on a silicon chip,” Opt. Lett. 29(20), 2387–2389 (2004). [CrossRef]  

12. Q. Xu and M. Lipson, “All-optical logic based on silicon micro-ring resonators,” Opt. Express 15(3), 924–929 (2007). [CrossRef]  

13. L. Fan, J. Wang, L. T. Varghese, H. Shen, B. Niu, Y. Xuan, A. M. Weiner, and M. Qi, “An all-silicon passive optical diode,” Science 335(6067), 447–450 (2012). [CrossRef]  

14. T. J. Johnson, M. Borselli, and O. Painter, “Self-induced optical modulation of the transmission through a high-q silicon microdisk resonator,” Opt. Express 14(2), 817–831 (2006). [CrossRef]  

15. L.-W. Luo, G. S. Wiederhecker, K. Preston, and M. Lipson, “Power insensitive silicon microring resonators,” Opt. Lett. 37(4), 590–592 (2012). [CrossRef]  

16. L. Zhang, Y. Fei, Y. Cao, X. Lei, and S. Chen, “Experimental observations of thermo-optical bistability and self-pulsation in silicon microring resonators,” J. Opt. Soc. Am. B 31(2), 201–206 (2014). [CrossRef]  

17. L. Zhang, Y. Fei, T. Cao, Y. Cao, Q. Xu, and S. Chen, “Multibistability and self-pulsation in nonlinear high-q silicon microring resonators considering thermo-optical effect,” Phys. Rev. A 87(5), 053805 (2013). [CrossRef]  

18. N. Takemura, M. Takiguchi, and M. Notomi, “Designs toward synchronization of optical limit cycles with coupled silicon photonic crystal microcavities,” Opt. Express 28(19), 27657–27675 (2020). [CrossRef]  

19. S. Li, S. Dev, S. Ohlendorf, K. Jamshidi, and S. Pachnicke, “Photonic reservoir computing enabled by silicon micro-rings,” in Asia Communications and Photonics Conference, (Optical Society of America, 2019), pp. M4C–4.

20. M. Catuneanu, R. Hamerly, N. Annavarapu, S. Sabouri, and K. Jamshidi, “Nonlinear activation function generation based on silicon microring resonators for integrated photonic neural networks,” in The European Conference on Lasers and Electro-Optics, (Optical Society of America, 2019).

21. F. Denis-Le Coarer, M. Sciamanna, A. Katumba, M. Freiberger, J. Dambre, P. Bienstman, and D. Rontani, “All-optical reservoir computing on a photonic chip using silicon-based ring resonators,” IEEE J. Sel. Top. Quantum Electron. 24(6), 1–8 (2018). [CrossRef]  

22. C. Mesaritakis, V. Papataxiarhis, and D. Syvridis, “Micro ring resonators as building blocks for an all-optical high-speed reservoir-computing bit-pattern-recognition system,” J. Opt. Soc. Am. B 30(11), 3048–3055 (2013). [CrossRef]  

23. T. Van Vaerenbergh, M. Fiers, P. Mechet, T. Spuesens, R. Kumar, G. Morthier, B. Schrauwen, J. Dambre, and P. Bienstman, “Cascadable excitability in microrings,” Opt. Express 20(18), 20292–20308 (2012). [CrossRef]  

24. T. Van Vaerenbergh, M. Fiers, J. Dambre, and P. Bienstman, “Simplified description of self-pulsation and excitability by thermal and free-carrier effects in semiconductor microcavities,” Phys. Rev. A 86(6), 063808 (2012). [CrossRef]  

25. M. Mancinelli, M. Borghi, F. Ramiro-Manzano, J. Fedeli, and L. Pavesi, “Chaotic dynamics in coupled resonator sequences,” Opt. Express 22(12), 14505–14516 (2014). [CrossRef]  

26. D. Dimitropoulos, R. Jhaveri, R. Claps, J. Woo, and B. Jalali, “Lifetime of photogenerated carriers in silicon-on-insulator rib waveguides,” Appl. Phys. Lett. 86(7), 071115 (2005). [CrossRef]  

27. I. Aldaya, A. Gil-Molina, J. L. Pita, L. H. Gabrielli, H. L. Fragnito, and P. Dainese, “Nonlinear carrier dynamics in silicon nano-waveguides,” Optica 4(10), 1219–1227 (2017). [CrossRef]  

28. W. H. Pernice, M. Li, and H. X. Tang, “Time-domain measurement of optical transport in silicon micro-ring resonators,” Opt. Express 18(17), 18438–18452 (2010). [CrossRef]  

29. V. R. Almeida, C. A. Barrios, R. R. Panepucci, M. Lipson, M. A. Foster, D. G. Ouzounov, and A. L. Gaeta, “All-optical switching on a silicon chip,” Opt. Lett. 29(24), 2867–2869 (2004). [CrossRef]  

30. J. Van Campenhout, W. M. Green, X. Liu, S. Assefa, R. M. Osgood, and Y. A. Vlasov, “Silicon-nitride surface passivation of submicrometer silicon waveguides for low-power optical switches,” Opt. Lett. 34(10), 1534–1536 (2009). [CrossRef]  

31. J. Linnros, “Carrier lifetime measurements using free carrier absorption transients. i. principle and injection dependence,” J. Appl. Phys. 84(1), 275–283 (1998). [CrossRef]  

32. W. Shockley and W. Read Jr, “Statistics of the recombinations of holes and electrons,” Phys. Rev. 87(5), 835–842 (1952). [CrossRef]  

33. K. Nomura and J. Blakemore, “Decay of excess carriers in semiconductors. ii,” Phys. Rev. 121(3), 734–740 (1961). [CrossRef]  

34. B. Birge, “Psot-a particle swarm optimization toolbox for use with matlab,” in Proceedings of the 2003 IEEE Swarm Intelligence Symposium. SIS’03 (Cat. No. 03EX706), (IEEE, 2003), pp. 182–186.

35. Q. Lin, O. J. Painter, and G. P. Agrawal, “Nonlinear optical phenomena in silicon waveguides: modeling and applications,” Opt. Express 15(25), 16604–16644 (2007). [CrossRef]  

36. S. M. Sze and K. K. Ng, Physics of semiconductor devices (John wiley & sons, 2006).

37. T. L. Bergman, F. P. Incropera, D. P. DeWitt, and A. S. Lavine, Fundamentals of heat and mass transfer (John Wiley & Sons, 2011).

38. https://www.comsol.eu.

39. X. Zheng, Y. Luo, G. Li, I. Shubin, H. Thacker, J. Yao, K. Raj, J. E. Cunningham, and A. V. Krishnamoorthy, “Enhanced optical bistability from self-heating due to free carrier absorption in substrate removed silicon ring modulators,” Opt. Express 20(10), 11478–11486 (2012). [CrossRef]  

40. C. Baker, S. Stapfner, D. Parrain, S. Ducci, G. Leo, E. M. Weig, and I. Favero, “Optical instability and self-pulsing in silicon nitride whispering gallery resonators,” Opt. Express 20(27), 29076–29089 (2012). [CrossRef]