On the modeling of thermal and free carrier nonlinearities in Silicon On Insulator microring resonators

The temporal dynamics of integrated silicon resonators has been modeled using a set of equations coupling the internal energy, the temperature and the free carrier population. Owing to its simplicity, Newton's law of cooling is the traditional choice for describing the thermal evolution of such systems. In this work, we theoretically and experimentally prove that this can be inadequate in monolithic planar devices, leading to inaccurate predictions. A new equation, that we train to reproduce the correct temperature behaviour, is introduced to fix the discrepancies with the experimental results. We discuss the limitations and the range of validity of our refined model, identifying those cases where Netwon's law provides, nevertheless, accurate solutions. Our modeling describes the phenomena underlying thermal and free carrier instabilities, and is a valuable tool for the engineering of photonic systems which relay on resonator dynamical states, such as all optical spiking neural networks or reservoirs for neuromorphic computing.


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 becomes 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,16,17,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,20,21,22,23]. The design and optimization of these systems require the accurate modeling of each node, which is traditionally implemented using 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 ∆ν of the pump laser [24,23,17]. 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 (∼ 100 ns) and a slow (∼ 1 µ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 (∼ 40 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. realized by evanescently coupling two bus waveguides to the MR of 7 µm of radius. This is accomplished through a 250 nm gap and a coupling length of 3 µm, which yields a measured coupling coefficient of κ 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 i = 1.11(8) × 10 5 and a loaded Q of Q L = 6.5(2) × 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 nm is amplified and injected at the input waveguide of the MR through a grating coupler (∼ 3.5 dB insertion loss). The input power is controlled by a Variable Optical Attenuator. Optionally, a C-band pulsed laser source of 40 ps pulse width and 1 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 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 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, ∆ν). 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, ∆ν) 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 resonace. Despite the similar MR geometry, this contrasts with the results in [16], where SP is exclusively detected at the blue side of the resonace. 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 ∼ 1 MHz, and can be as low as ∼ 400 kHz for large positive detunings. This is an order of magnitude lower than the SP frequencies reported in [25,16], despite the very similar MR geometries. In general, SP periods of ∼ 100 ns arise when the thermal and the carrier lifetime are respectively of the order of τ th ∼ 80 − 150 ns and τ fc ∼ 1 − 10 ns [24]. This is because a complete cycle of oscillation requires both thermal and carrier relaxation, so its period is roughly given by τ fc + τ 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.

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 mm, located on the same die of the MR. As shown in Fig.1(b), a strong pulsed pump (40 ps) with ∼ 1.5 W peak power and a wavelength of 1550 nm is combined with the CW probe and injected into the spiral. The high intensity of the pump creates an initial free carrier population ∆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 ∆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 τ fc , using the approach detailed in [27]. Figure 3(b) shows τ fc as a function of the carrier concentration ∆N inside the waveguide. The two red dashed lines indicate the expected profile for a purely exponential decay, i.e., log(∆N ) = log(∆N 0 ) − t τfc . The line indicated as τ slow ∼ 60 ns is obtained by evaluating the average local slope for 0 ≤ t ≤ 50 ns, while the time interval 250 ≤ t ≤ 300 ns is used to evaluate the slope of the line labeled as τ fast ∼ 25 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 τ 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.
The value of τ fc is bounded between 20 ns and 60 ns. As shown in Fig.3(b), at the very early stage, carrier recombines with a characteristic lifetime of τ slow ∼ 60 ns, which is significantly larger than the one in the late stage (τ fast ∼ 25ns), in which the carrier concentration has dropped by nearly two orders of magnitude. Even though the instantaneous lifetime is larger than the ones traditionally reported for SOI waveguides [28,29,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 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. Trap-assisted recombination leads to electron-hole lifetimes which depend on the amount of excess carriers [31,32,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]: where R(t) is the carrier recombination rate, τ 0 is the lifetime in the low carrier injection regime and (a, b) are constants which depend on material properties, such as the doping level, the electron-hole capture cross-section and the trap energy level. The definition and the numerical values of these terms is reported in Table 2 of Appendix B. Equation (1) is used to model the points in Fig.3(b), indicating a good agreement with the experiment. In principle, if the density of traps is of the same order of magnitude of the equilibrium concentration of the majority carriers, the phenomenon of trapping makes electrons and holes to recombine at different rates, so one has to separately describe the evolution of the two types of carrier [27,33]. As a consequence, τ fc is not solely determined by the carrier concentration, but also by the initial excess carrier density. An experimental signature of the presence of trapping is the dependence of the carrier lifetime on the pump power that is used to create the initial excess population [27]. Through numerical simulations, we verified that in the concentration range 10 15 − 10 18 cm −3 , which is of concern in our work, the phenomenon of trapping has a negligible impact on τ fc .

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 int = |a| 2 of the resonator, the free carrier population ∆N , and its differential temperature with the cladding ∆T [14]: where the meaning and the values of the various symbols are given in Table 1 of Appendix A. From Eq.(2), the power P d at the Drop port is obtained by P d = γ e U int and P abs = 2γU int . As a first step, we replaced the constant free carrier lifetime τ fc in Eq.(3) with the expression in Eq.(1). We then numerically integrated Eqs.(2-4) by changing τ 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 dBm and ∆ν = −10 GHz. The best value of τ 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 τ th = 270(10) ns. Remarkably, we had to increase the FCD coefficient up to σ FCD = −9.5 × 10 −27 m 3 in order to observe SP, which is ∼ 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 Fig.4(a,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 maps 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 (σ FCD , τ th ), we generated the same maps in Fig.4(a,b), and we evaluated the overlap with the ones in Fig.2(a,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 ∝ cp τth( dn dT ) √ σ FCD τ fc , where dn dT is the thermo-optic coefficient and c p the specific heat. Since τ 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, ∆ν). This produced the values of c p (P, ∆ν) and τ th (P, ∆ν) 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 ∼ 7 − 10 times higher than its tabulated reference value c p,ref = 700 J/(kgK) . Furthermore, τ th lies between 350 − 650 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 τ 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, ∆ν). Using the average values τ th = 450 and c p = 8c p,ref , we obtain q ∼ 5q 0 , where q 0 is the value of q generating the maps in Fig.4(a,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 ∼ 3.5 times higher than the ones reported in [25,16], 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.

The 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 mW for 20 µ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 ∆T is not a single exponential. At this point, one could fix the discrepancy by directly coupling Eq.(2,3) to the FEM solver, and integrate the temperature profile inside the MR to lump the spatial variables, i.e., ∆T (x, y, z, t) → ∆T (t). This is certainly possible, but on one hand is time consuming, and on the other 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 ∆T = C s exp − t τs +C f exp solution. Here, we call τ s τ 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 . Equation (7) is homogeneous, so we must provide the driving terms related to the absorbed power. Heuristically, we propose the form: which finds justification from the following considerations. At steady state, the solution predicted by Eq. (7) is ∆T = κτsτ f Pabs mcp , which is analogous to the solution of Eq.(4) if we set τ f τ s κ = τ th . As shown in Fig.5(a), d∆T dt must be discontinuous at time t 0 , that is when P abs instantaneously drops from the ON to the OFF state. We can integrate Eq.
where we have used the fact that dPabs dt = ηδ(t − t 0 ). The parameters (κ, η, τ s , τ f ) are found by numerically integrating Eq.(7) and by simultaneously minimizing the least squares with the FEM temperature profile using a PSO. A realistic form of P abs can be obtained by integrating Eqs.(2-4) using any of the values of c p and τ th shown in Fig.4(d). We used as reference P = 8 dBm and ∆ν = 0 GHz. Figure 5(b) shows a very good agreement between the temperature solution of the FEM and the one obtained by integrating Eq. (7). As a comparison, on the same plot we reported two curves, ODE1 (1) and ODE1 (2), which we obtain by integrating Eq.(4) and by setting respectively (c p,ref , τ th = 40 ns) (the choice of τ th = 40 ns minimizes the error with the FEM trace) and (c p = 6.9c p,ref , τ th = 417 ns) (i.e., the pair of values in Fig.4(d)). The single order ODE is subjected to a trade-off between the decay time and the average temperature that is reached at stationary regime. The curve at τ th = 40 ns would benefit from a longer decay time, but this comes at the cost of an higher average temperature T av , since T av ∝ Pabsτth cp , which worsen the similarity with the FEM trace. This can be compensated by simultaneously increasing c p (solid red curve in Fig.5(b)), which is the origin of the anomalously high values reported in Fig.4(d). The tradeoff is removed in Eq.(7) by the inertial term 1 τsτ f ∆T , which decouples the decay rate from the temperature at steady state. Equation (7) also explains why the values of c p and τ th are correlated to the SP period. As shown in Fig.5, the decay of ∆T has two characteristic lifetimes: one τ f governing the initial stage, where the temperature is dropping fast, and another τ s which takes over for t > τ f , creating a smooth decay tail. In the model of Eq.(4), ∆T decays as a single exponential, so ideally the match with Eq. (7) is lost as soon as t > τ f , that is, when the slow exponential takes over. This can be partially recovered by increasing τ th as the overall decay time increases (i.e., the SP period), which compensates the absence of the slower exponential term. Concurrently, c p has to be increased to not alter the average temperature, which is why it gets linearly correlated with τ th , as shown in Fig.4(d).
In a last step, we learned the parameters (κ, η, τ s , τ f ) using the experimental time response of the MR. This is accomplished by running a multi-objective PSO optimizer simultaneously minimizing the least square error between the MR output predicted by Eqs. (2,3,7) and three arbitrarily chosen experimental traces in the plane (P, ∆ν). The result of the optimization yield (κ = 7.45(5) × 10 −4 ns −1 , τ f = 71(2) ns, τ s = 1030(20) ns, η = 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 Fig.6(a,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 p,ref ).
6 Accuracy of the temperature equation Equation (7) fixes the discrepancies between the model and the experiment by introducing an inertial term proportional to ∆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 mW of amplitude, 100% extinction and T ref = 1 µs of period. These conditions well approximate the temporal profile of P abs during the SP regime of a MR. We then changed the period of the square wave to some other value T sp = T ref , and used the FEM solver to evaluate the maximum temperature excursion ∆T max for each SP cycle (initial transient effects are discarded in this analysis). For example, in Fig.5(b), ∆T max ∼ 6 K for T sp = 0.8 µs. Note that ∆T max is a function of T sp . We assume that the temperature variations predicted by FEM represent what occurs in a MR with a SP period of T sp . Using the same driving power P abs in Eq. (7), we obtained the approximated ODE solution, from which we evaluated the time T ODE required to reach the same temperature variation ∆T max . The error on the predicted SP period by Eq. (7) is then determined by T sp − T 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 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 sim predicted by the full model given by Eqs. (2,3,7). To this purpose, we used the set of parameters (κ, η, τ s , τ 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 exp − T sim as a function of T sp are shown in Fig.7 (black diamonds). As a comparison, we also included the errors on the SP period T exp − T 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 µs). If we set the desired accuracy threshold to 10%, the model loses accuracy as long as the period exceeds 2 µs or gets lower than 500 ns. This interval is sufficiently large to accommodate all the SP frequencies reported in Fig.2(b).

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 ∼ 2 − 3 times higher(lower). It is then expected that different MR have to be modeled by different sets (κ, η, τ s , τ 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,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, ∆ν), 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 τ 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 p,ref , but we will have to artificially increase τ fc and/or σ 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 τ 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 τ 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.

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, ∆ν) 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)(3)(4), with the exception of τ fc and τ th , that are defined in the main text. The resonance shift induced by free carriers is modeled as to produce the maps in Fig.4(a,b). The refined expression δ fc = − σFCD,1∆N +σFCD,2∆N 0.8

n0
[35] is used to produce the maps in Fig.6(a,b).

Appendix B: Values of the terms in the Shockley Read Hall expression for τ 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 τ n (τ p ). The terms (a, b, τ 0 ) appearing in Eqs. (2)(3)(4) have the following definitions: a = τ p + τ n τ p (n 1 + n 0 ) + τ n (p 1 + p 0 ) b = 1 n 0 + p 0 τ 0 = τ p n 0 + n 1 n 0 + p 0 + τ n p 0 + p 1 n 0 + p 0 (8) 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).