Physical origin of higher-order soliton fission in nanophotonic semiconductor waveguides

Supercontinuum generation in Kerr media has become a staple of nonlinear optics. It has been celebrated for advancing the understanding of soliton propagation as well as its many applications in a broad range of fields. Coherent spectral broadening of laser light is now commonly performed in laboratories and used in commercial “white light” sources. The prospect of miniaturizing the technology is currently driving experiments in different integrated platforms such as semiconductor on insulator waveguides. Central to the spectral broadening is the concept of higher-order soliton fission. While widely accepted in silica fibers, the dynamics of soliton decay in semiconductor waveguides is yet poorly understood. In particular, the role of nonlinear loss and free carriers, absent in silica, remains an open question. Here, through experiments and simulations, we show that nonlinear loss is the dominant perturbation in wire waveguides, while free-carrier dispersion is dominant in photonic crystal waveguides.


INTRODUCTION
Solitons are ubiquitous in science. They correspond to localized packets that propagate unperturbed as a consequence of a balance between nonlinear self-focusing and a diffusionlike process. They have been theoretically and experimentally investigated in hydrodynamics [1], plasma physics [2], biology [3] and optics [4]. In the latter, the sech-shaped solutions of the nonlinear Schrödinger equation (NLSE), known as solitons, are the most notable [5].
Initially considered as potential carriers of information, the focus recently shifted to the central role they play in the process of nonlinear spectral broadening [6]. NLSE solitons can be of different orders, with only the first order (fundamental) soliton strictly maintaining a constant profile during propagation. Higher-order solutions instead display a periodically varying shape [7]. When perturbed, higher-order solitons tend to split into several fundamental solitons, each centered at a different wavelength, potentially spanning several octaves. This process, commonly called higher-order soliton fission, is the basic mechanism underlying supercontiuum generation technology that finds application in a broad range of fields [8]. Initially performed in bulk crystals [9], supercontinuum generation has been promoted by the advent of photonic crystal fibers [10]. The long interaction lengths and small mode areas facilitated the observation of soliton fission and triggered a wealth of experimental investigations [8]. While any perturbation to the NLSE will eventually lead to the decay of a higher-order soliton, it is now well understood that the Raman effect and higher-order dispersion (HOD) govern the dynamics of the fission process in optical fibers.
As there currently is a big interest in the miniaturization of supercontinuum sources, many experimental investigations of spectral broadening through fission of higher-order solitons in nanophotonic waveguides have been reported [11][12][13][14][15][16][17]. With most of the focus on the spectral broadening, discussions on the dynamics of the fission remain scarce. The broadest spectra have been obtained on platforms made of wide band-gap media such as silicon nitride [18] and chalcogenide glasses [19] where the nonlinear loss is low. The dynamics of soliton fission in these waveguides is inferred to be very similar to that of optical fibers.
In contrast, soliton fission in narrow band-gap semiconductors such as silicon and indium gallium phosphide (InGaP) is much less understood. The band structure of the medium complexifies the dynamics as nonlinear absorption and free carriers may impact soliton propagation [20,21]. Supercontinuum generation through fission of higher-order solitons in semiconductor waveguides was first proposed in 2007 [22]. Other theoretical studies followed [15,23,24]. Few however discuss the impact of individual perturbative terms on the soliton dynamics.
In our recent report on supercontinuum generation in silicon wires [25] pumped around the 1550 nm telecommunication wavelength, we argued that the main perturbation to the NLSE describing our measurements is two-photon absorption (2PA). Soliton fission provoked by nonlinear loss was shown numerically by Silberberg [26], and in early experimental investigations of one-dimensional spatial soliton propagation in solid-state waveguides [27,28].
Temporal and spatial solitons being mathematically equivalent, 2PA is likely to play an important role in the temporal domain as well. Simulations describing our experiments indeed showed a fairly symmetrical post-fission dynamics indicative of 2PA, barring the emission of a weak dispersive wave (DW). It is well known that solitons can shed energy when perturbed by HOD [29]. We inferred when analyzing our data that the weak symmetry breaking we observed was solely due to DW emission. This lead us to the conclusion that 2PA is the main cause for soliton decay in a silicon wire pumped around 1550 nm [25,30].
In contrast, another recent report claims to experimentally observe free-carrier induced fission of higher-order solitons in InGaP photonic crystal waveguides (PhCWG) [31]. In the supplementary section, through analytical investigation of the perturbed NLSE, the authors suggested that free-carrier dispersion (FCD) can also trigger soliton decay in 1550 nm pumped silicon wires. However, the authors did not consider the role of the wire waveguide geometry, nor 2PA interplay with the FCD effect, and the analysis appears incomplete. Free carriers are well known to impact pulse propagation in integrated waveguides [32]. The nonlinear loss inherent to semiconductor devices results in the excitation of electron-hole pairs that subsequently impact the dynamics in two different ways. They absorb photons (free carrier absorption, FCA) as well as change their wavenumber (free carrier dispersion, FCD) [33]. The latter is the mechanism underlying carrier depletion semiconductor phase and amplitude modulators [34]. It can also alter short pulse propagation through phase modulation on the leading edge of the pulse and was demonstrated to impact soliton propagation [35]. FCA on the other hand has less impact on soliton propagation as the induced loss is low compared to the instantaneous multi-photon absorption processes (see below).
The conclusions of [25] and [31] are difficult to reconcile, leaving the reader wondering what drives the dynamics of ultrashort pulse propagation in semiconductor nanowaveguides.
Here, we aim to clarify the role that both nonlinear loss and FCD play on higher-order soliton fission. We start by investigating the case of pulse propagation in a silicon wire waveguide pumped around the telecommunication wavelength where nonlinear loss is dominated by 2PA. We perform new experiments of supercontinuum generation and confirm that FCD has very little impact on the dynamics on the wire waveguide geometry. Furthermore we theoretically investigate the case of InGaP waveguides. When pumped at 1550 nm, the lowest order nonlinear loss is three-photon absorption (3PA). We review the impact of 3PA on the propagation of a higher-order soliton. We find that nonlinear loss can induce soliton fission in InGaP wire waveguides and that the impact of carriers strongly depends on the waveguide dispersion and input pulse characteristics. From these broad investigations, we derive a general set of conclusions about the roles of nonlinear loss and FCD in higher-order soliton fission in nanophotonic semiconductor waveguides.

RESULTS
Soliton fission in 2PA-limited silicon nanowires. Experiments and numerical simulations. In previous experimental investigations of soliton fission, the main objective was to demonstrate that broad supercontinuum generation can be achieved despite strong nonlinear loss [25,30]. They were therefore focused on waveguides whose dispersion properties maximized the spectral broadening. Specifically, it was shown that carefully choosing the waveguide dimensions allows to precisely control the emission of DWs. In these measurements, DWs were emitted on the blue side of the soliton only, inducing a spectral asymmetry that may be misinterpreted as FCD induced blue shift. Here, we report new experiments performed in silicon wire waveguides designed to have notably different dispersion profiles.
The objective of these experiment is to further probe the role of carriers on supercontinuum generation in the presence of strong 2PA.
We use 6 mm long, 220 nm high silicon-on-insulator wire waveguides. We consider six different waveguides, with widths ranging from 800 nm to 575 nm. We inject 165 fs pulses at 1556 nm with a peak power of 30 W at the input of the nanowires (See methods).
The measured output spectra are shown in Figure 1 (blue curves). We readily observe (i) evidence of soliton fission and dispersive wave emission, and (ii) that the broadening highly depends on the waveguide width. In order to gain insight into our experimental results we perform numerical simulations of the generalized nonlinear Schrödinger equation (GNLSE) describing ultra-short pulse propagation in silicon nanowires [See Equation (5) in methods]. The dispersion profile is computed by use of a mode solver (see methods).
We stress that calculated dispersion profiles are known to deviate from measured ones. For example, the dispersion of a 220 nm high silicon wire waveguide with a width of 500 nm (as measured in [16]) is best reproduced by mode solver simulations when considering a 420 nm wide waveguide with the same height. We account for these deviations by treating the waveguide width as a free parameter in the simulations. Note that all other parameters are fixed. The computed spectra are shown in Figure 1 and are in excellent agreement with the experiments. Also shown in Figure 1 are the dispersion profiles used in simulations.
Detailed spatial evolution of the simulated (temporal and spectral) pulse profiles can be found in the supplementary materials.
The spectrum at the output of the first waveguide [ Figure 1(a)] is consistent with a propagation in the normal dispersion regime as predicted by the computed dispersion. The spectral broadening is induced by self-phase modulation only and the pulse temporally spreads as it propagates. As the input spectrum moves into the anomalous dispersion regime [ Figure 1(b)], we notice an increased but asymmetric spectral broadening. It stems from higher-order soliton compression and subsequent emission of a blue detuned dispersive wave [25]. The position of the latter is well predicted by the zero crossing of the linear wavenumber as expected from theory [29,36]. As we further reduce the waveguide width, the zero-dispersion wavelength (vertical line) moves away from the pump and the dispersive wave progressively shifts to shorter wavelengths. In the narrowest waveguide [ Figure 1(f)], the phase matching point is located well below the band-gap of silicon (≈ 1120 nm) such that the process is hampered.
The excellent agreement between experimental and numerical results constitute a good starting point for discussing the impact of the different perturbations on soliton propagation. We start by introducing a model solely comprising HOD, linear loss and two-photon absorption as perturbations to the NLSE. This equation, describing the evolution of the temporal envelope E(z, t) of the electric field along z, reads: where β k are the Taylor series expansion coefficients of the dispersion profiles shown in Figure 1, α l and α 2PA account respectively for the linear and nonlinear losses and γ is the nonlinear parameter of the waveguide. We compare (in Figure 1) the spectra computed with equation (1) with those obtained with equation (5) that includes among others the FCD term. We find only minor differences between the two, hinting that only 2PA and HOD impact pulse propagation (linear loss is negligible over such a short distance). We contend that the agreement validates the use of equation (1) for ultra-short pulse propagation in silicon waveguides in the presence of strong 2PA, as predicted in [22]. Whether it is 2PA or HOD that dominates is likely to be sensitive to the experimental conditions. In optical fibers for example, Raman and HOD introduce comparable perturbations and which induces the fission highly depends on input pulse duration [8]. While we found that 2PA and HOD suffice to explain our results, it seems pertinent to investigate the impact of each perturbation independently in an effort to gain further physical insight into the dynamics of soliton fission in semiconductor waveguides. Before considering each perturbation independently, we will consider the full GNLSE equation (5) for the sake of completeness.
For the full GNLSE simulations, we focus on the dispersion profile corresponding to the narrowest waveguide where no DW is emitted to clearly discriminate between the impact of HOD on the fission and DW emission [see Figure 1(f)]. The temporal profile evolution as predicted by equation (5) is shown in Figure 2(a). The propagation is reminiscent of the higher-order soliton fission dynamics observed in optical fibers, with signatures of (i) an early temporal compression stage and (ii) subsequent splitting into several pulses. We note that due to the strong nonlinear loss, it is not clear if these pulses qualify as solitons. It was suggested that they can be described using the concept of path-averaged solitons [20].
In contrast, when nonlinear loss is absent (as discussed below), we have confirmed that the Linear loss is systematically included in the simulations but we note that while it may induce fission on its own (see e.g. [37]), we have not observed such behavior over the length of the sample. In Figure 2  While it is apparent that we can describe this current device with only 2PA, we continue our examination of the different perturbations. If we include FCD only [ Figure 2(d)], the dynamics is very different. It is then governed by soliton acceleration as a consequence of carrier-induced blue shift. This is consistent with the dynamics simulated in [31] and confirms that the small acceleration observed in Figure 2(a) is due to the presence of carriers.
In Figure 2(e) we show the impact of HOD. It is well known that HOD may both provoke soliton decay and lead to DW emission [8]. Remember that we purposely investigate the case where no DW is emitted in the experiment to be able to discriminate between the two.
Without nonlinear loss however, the spectrum of the higher-order soliton is much broader at the point of compression such that it emits DW's in the normal dispersion regime on both sides of the pump. We expect that the dynamics shown in Figure 2(e) captures the main features of higher-order soliton propagation in waveguides with negligible nonlinear loss and Raman effect such as SiN-on-insulator nanowires [14,38]. In silicon however, we see no sign of the impact of HOD on the propagation outside of DW emission.
The fission with Raman only is shown in Figure 2 compete with 2PA induced loss. In brief, it is apparent that many of these perturbations could induce fission independently. In our case, however, the simulations indicate that the only perturbation to have a significant impact on the dynamics is 2PA. While in the above discussion we focused on a single set of parameters, we note that this conclusion holds for the other cases shown in Fig. 1 and in 2PA-limited wire waveguides.
Dimensionless analysis Given the strong 2PA in our system, we now examine whether FCD can be the main perturbation under different conditions. We here focus on the fission dynamics under the combined action of 2PA and FCD. We start by formulating a normalized model that accounts for both effects. All other perturbations are omitted. Moreover we neglect carrier recombination as the corresponding lifetime is much longer than the input pulse duration. The equation describing the normalized field envelope A(ζ, τ ) evolution during propagation (in the anomalous dispersion regime) then reads: where ζ and τ are the normalized distance and time parameters (see methods). The terms on the left hand side correspond to the standard NLSE. On the right hand side are the 2PA and FCD perturbations. The latter is proportional to the normalized parameter η.
It is equal to the ratio of the 2PA and FCD lengths (as defined in 31) and only depends on the free-carrier index change coefficient k c , the effective mode area, and the input pulse energy (see Methods). For a fixed material, the higher the effective fluence [Jm −2 ] in the waveguide, the more impact carriers will have. At this point we note that our normalized equation differs from the one used in [31] where 2PA was neglected in order to focus the analysis solely on FCD.
2PA and FCD effects are difficult to compare as they impact pulse propagation in very different ways. The former results in instantaneous, power-dependent loss while the latter produces a non-instantaneous index variation. We infer however that the normalized amplitude of the corresponding term in equation (2) is likely to be a good indicator of FCD impact. We hence compare the strength of both perturbations as a function of time (τ ) and distance (ζ). We consider the parameters of the experiment shown in Figure 1(f). The input pulse corresponds to a soliton of order N = 6. We integrate equation (2) over the waveguide length [starting from A(0, τ ) = sech(τ )] and plot, in Figure 3, the evolution of the amplitude of the 2PA (a) and FCD (b) perturbations. We readily note that both are strongly dependent on time and propagation distance. While FCD is stronger than 2PA at the beginning of the propagation, it decays more rapidly as the pulse advances. We believe that this is the reason we hardly see any impact of the carriers on the dynamics in our current experiments. Similarly, using characteristic length scales in the presence of nonlinear loss can be misleading since they are calculated based on input parameters only, and not dynamic interaction of the nonlinearities.
Next we investigate the potential impact of carriers on soliton propagation. Specifically, we compute, by use of equation (2), the temporal pattern at the output of the waveguide for different values of the normalized FCD parameter η. Our results are shown in Figure 3(c).
The dotted line indicates the value computed with our experimental parameters (η = 1.78). crystals [24]. This can be more easily understood by rewriting the dimensionless FCD parameter as where the constants common to silicon waveguides pumped around 1550 nm and the soliton number clearly appear, leaving the |β 2 |/T 0 ratio as an important parameter allowing to predict the asymmetry. This is illustrated in Figure 4(a) where we plot η as a function of the |β 2 |/T 0 ratio for a silicon waveguide pumped around 1550 nm with a sixth-order soliton.
As discussed above, our experiments correspond to η ≈ 2 and display a mostly symmetrical output wave. In photonic crystal waveguides however, the group-velocity dispersion is substantally larger because slow propagation arises close to the photonic band edges. Despite the longer pulses, these experiments are characterized by a |β 2 |/T 0 ratio two orders of magnitude larger than for channel waveguides. The lack of demonstration of FCD induced soliton fission in silicon photonic crystals is likely due to the strong linear loss [24]. To further stress how important the dimensionless parameter η is for predicting the dynamics, Our choice to describe the asymmetry in the Fourier domain is motivated by the fact that spectrum does not change noticeably post-fission (See e.g. Figure S1). Our results in Figure 4(b) show that the soliton number has very little impact on the spectral asymmetry of the output wave confirming that the normalized parameter η, proportional to the input energy, is the most relevant parameter to determine the dominant perturbation mechanism.
Soliton fission in 3PA-limited InGaP photonic crystal and wire waveguides. Numerical simulations. We now turn our attention to the dynamics of soliton fission in InGaP waveguides as recently reported in [31]. The experiments are performed by pumping a In Figure 5(a) we show a map of the propagation through the waveguide when simulated by the GNLSE. Note that it is a different GNLSE than the one used for simulating pulse propagation in a silicon waveguide (see [31] and methods for more details). We observe a simultaneous acceleration and compression of the input soliton followed by a decay into two separate pulses. Dispersive waves are emitted at the point of maximum compression.
Next we simulate the dynamics with FCD as the lone perturbation. The temporal profile is shown in Figure 5(b). We see both a split of the input soliton and acceleration of the pulses. However, the lack of nonlinear loss noticeably affects the temporal trajectory of the pulses. The increased carrier-induced blue shift causes the first ejected soliton to propagate much faster. The dynamics when we take only 3PA into account is displayed in Figure 5(c).
Interestingly, we remark that the higher-order soliton still decays. This shows how 3PA, just as 2PA, impacts the dynamics beyond simply transferring energy from the optical field to the carriers. As expected, the temporal profile fully satisfies time reversal symmetry throughout propagation. Finally, in Figure 5(d) we show the dynamics when HOD is the lone perturbation to the NLSE. Again, fission occurs. The case with TOD only was already discussed in [31] and displays no fission at all for the low TOD values in that experiment.
The emission of dispersive waves on both sides of the pump is also an indication that FOD is the main linear perturbation [41]. Importantly, we find that all three perturbations are capable of inducing fission, similar to the 2PA case above. Looking at Figure 5 it seems likely that all play a part, and we hence further investigate the role of 3PA and in particular how it compares with FCD.
Dimensionless analysis Similarly for 2PA, we introduce a normalized model including only nonlinear loss and FCD: The perturbation terms on the right hand side are different than in equation (2)  for different values of κ. We start with the normalized parameters N = 2 and Γ = 2 × 10 −2 computed from the physical constants used in the simulations shown in Figure 5. The results are shown in Figure 6(a). For very low values of the FCD parameter κ, the output profile is the same as in Figure 5(c). Even when the impact of carriers is negligible, the input pulse splits at the first point of compression. The fission is induced by 3PA for these low κ values. As we increase κ, the carrier-induced blue shift breaks the time reversal symmetry of the pattern. Most of the energy moves to the front pulse and the distance between the two pulses increases. The value of κ corresponding to the simulations shown in Figure 5 is highlighted by the dotted line. Such an important difference indicates that it is indeed FCD that drives the fission process in these experiments in photonic crystal waveguides. We stress that the normalized variables are highly dependent on the waveguide and input pulse characteristics (see methods).
We next probe the impact of carriers in a different geometrical configuration involving InGaP wire waveguides. Specifically, we simulate the propagation in waveguides similar to those recently used for octave-spanning supercontinuum generation [17]. The physical parameters can be found in the methods. The corresponding nonlinear absorption parameter is Γ = 1.4 × 10 −3 . We use the same normalized distance as before, corresponding to 4.5 dispersion lengths but we set the input soliton order to N = 4 because soliton fission does not occur for lower orders over this distance. We vary κ keeping Γ fixed and integrate equation (4) over the waveguide length. The results are shown in Figure 6(b). The output profile is still very symmetrical for the physical FCD value (dotted line). A ten fold increase of κ is required for the carriers to noticeably impact the propagation. In contrast to the photonic crystal case above, these simulations suggest that 3PA is likely to be the main driver of soliton fission in InGaP wire waveguides. Similar to the above analysis for the case of 2PA driven silicon waveguides, we now study the impact of the normalized parameters (here κ, Γ and N) on the spectral asymmetry of the output wave. Our results, shown in  Note that a fully symmetric spectrum has a null average frequency.

DISCUSSION
We unveiled the physics driving soliton fission in semiconductor waveguides. This research was motivated by recent results in the literature, each claiming different physical origin: 2PA [25] or FCD [31]. Through a series of experiments and accompanying analysis, we showed that the impact of 2PA and 3PA on ultra-short pulse propagation goes beyond energy loss and can induce fission. We further showed that either 2PA or FCD can be the dominant fission mechanisms depending on the experimental parameters. Our dimensionless analysis suggest that the normalized FCD parameters (η for 2PA driven waveguides, κ for 3PA driven waveguides) may be used to identify the dominant perturbation to the propagation of higher order solitons in semiconductor waveguides. We found that nonlinear loss (2PA/3PA) are the dominant perturbations in wire waveguides, while FCD is dominant in photonic crystal waveguides. We opened this discussion considering the important role of soliton fission in supercontinuum generation. As we often desire supercontinua? with a very broad bandwidth, it is much more likely that wire waveguides will be the primary geometry of interest. A key implication of our work is that nonlinear loss is the main driving mechanism in these systems. The main difference with fission induced by other known perturbations, such as the Raman effect or HOD, is that fission induced by nonlinear loss exhibits symmetric spectral and temporal evolution throughout propagation. Modeling of pulse propagation in silicon waveguides We use the following generalized nonlinear Schrödinger equation to simulate pulse propagation in the waveguides [22].

Experiment in silicon waveguides
where E(z, t) describes the slowly varying envelope of the field as a function of propagation distance z and time t andẼ(z, ω) = F [E(z, t)] is the Fourier transform of the field.
is the frequency dependent wave vector and ω 0 denotes the pump frequency. β(ω) is computed by use of a mode solver (Lumerical). The data is available on request. α l characterizes linear loss. It was evaluated at 2 dB/cm by cutback measurements on similar waveguides. σN c correspond to the free carrier induced loss where σ = 1.45 × 10 −21 m 2 [32] and N c is the carrier density, computed by solving: where A 3eff is the effective area related to a third order process, computed with the mode solver. The nonlinear parameter γ = n 2 ω 0 /(cA 3eff ) and 2PA coefficient are extracted from measurements on similar waveguides [42] and scaled through the effective mode area. As the latter is different for each waveguide, so are the nonlinear parameter and is the delayed nonlinear function where f R is the fractional Raman contribution and h R (t) is the Raman response function. They are deduced from the spectral Lorentzian shape of the Raman response of silicon [32]. We use the relation f R = g R (ω 0 )Γ R /[Ω R A eff γ] with g R (ω 0 ) = 3.7 × 10 −10 m/W [43], Ω R /(2π) = 15.6 THz and Γ R /π = 105 GHz [32]. The normalized equation (2) is linked to equation (5) through the following relations: where L D = T 2 0 /|β 2 | and L NL = 1/(γ R P 0 ) are the dispersion and nonlinear lengths as defined at the input of the waveguide, considering an input pulse √ P 0 sech(t/T 0 ). Note that for the normalized model, HOD, FCA, self-steepening as well as the Raman effect are neglected. Moreover we consider that the FCD parameter k c remains constant throughout the waveguide. We use k c = 3.8 × 10 −27 m 3 , computed with the input pulse parameters. Modeling of pulse propagation in InGaP photonic crystal and wire waveguides.
We use the following equation to simulate pulse propagation in InGaP waveguides [31]: where E(z, t) describes the slowly varying envelope of the field as a function of propagation distance z and time t. The equation differs from equation (5). The 2PA coefficient α 2PA is absent as there is no two-photon absorption when pumping around the 1550 nm telecommunication wavelength. Instead the 3PA coefficient α 3P A is included. Also, the FCD parameter k c is considered independent of N c and we neglect the Raman effect as well as self-steepening.
The carrier density is computed by solving where A 5eff is the effective area related to a fifth order process. The values used in our simulations of photonic crystals are listed in table II. The normalized equation (4) is linked to equation (7) through the relations: where L D = T 2 0 /|β 2 | and L NL = 1/(γP 0 ) are the dispersion and nonlinear lengths as defined at the entrance of the waveguide. NA: not applicable. * francleo@ulb.ac.be