I. INTRODUCTION

Section:

ChooseTop of pageABSTRACTI. INTRODUCTION <<II. PROPAGATION MODELIII. COMPUTATIONAL SCHEMEIV. COMPARISON WITH EXPER...V. CONCLUSIONSUPPLEMENTARY MATERIALREFERENCESBroadband optical spectra find a wide range of applications in optical sensing, spectroscopy, tomography, ultrashort pulse generation, and optical precision metrology from terahertz to ultraviolet frequencies.11. J. M. Dudley, G. Genty, and S. Coen, “Supercontinuum generation in photonic crystal fiber,” Rev. Mod. Phys. 78(4), 1135–1184 (2006). https://doi.org/10.1103/revmodphys.78.1135 Often, laser sources with the desired properties cannot operate directly in the spectral region of interest. Instead, spectral content at the desired wavelength may be generated from a pump laser source through nonlinear optical frequency conversion. Sum-(SFG) and difference-frequency generation (DFG) in non-centrosymmetric optical crystals have enabled, for instance, optical parametric amplifiers and oscillators in the mid-infrared and terahertz regimes, as well as high-harmonic and ultraviolet pulse generation. Extreme spectral broadening of optical pulses, in some cases limited only by material transparency, has been achieved through the effect of supercontinuum generation based on the cubic nonlinearities in step-index, photonic crystal, and specialty fibers.22. T. Sylvestre, E. Genier, A. N. Ghosh, P. Bowen, G. Genty, J. Troles, A. Mussot, A. C. Peacock, M. Klimczak, A. M. Heidt, J. C. Travers, O. Bang, J. M. Dudley, O. Bang, O. Bang, and J. M. Dudley, “Recent advances in supercontinuum generation in specialty optical fibers [Invited],” J. Opt. Soc. Am. B 38(12), F90–F103 (2021). https://doi.org/10.1364/josab.439330More recently, chip-based integrated micro- and nano-photonic waveguides have emerged as a powerful addition to established nonlinear media.3–73. M. R. Lamont, B. Luther-Davies, D.-Y. Choi, S. Madden, and B. J. Eggleton, “Supercontinuum generation in dispersion engineered highly nonlinear (γ = 10/W/m) As2S3 chalcogenide planar waveguide,” Opt. Express 16(19), 14938–14944 (2008). https://doi.org/10.1364/oe.16.0149384. K. Prakash, D. R. Solli, and B. Jalali, “Limiting nature of continuum generation in silicon,” Appl. Phys. Lett. 93(9), 091114 (2008). https://doi.org/10.1063/1.29778725. D. Duchesne, M. Peccianti, M. R. E. Lamont, M. Ferrera, L. Razzari, F. Légaré, R. Morandotti, S. Chu, B. E. Little, and D. J. Moss, “Supercontinuum generation in a high index doped silica glass spiral waveguide,” Opt. Express 18(2), 923–930 (2010). https://doi.org/10.1364/oe.18.0009236. B. Kuyken, X. Liu, R. M. Osgood, R. Baets, G. Roelkens, and W. M. J. Green, “Mid-infrared to telecom-band supercontinuum generation in highly nonlinear silicon-on-insulator wire waveguides,” Opt. Express 19(21), 20172–20181 (2011). https://doi.org/10.1364/OE.19.0201727. R. Halir, Y. Okawachi, J. S. Levy, M. A. Foster, M. Lipson, and A. L. Gaeta, “Ultrabroadband supercontinuum generation in a CMOS-compatible platform,” Opt. Lett. 37(10), 1685–1687 (2012). https://doi.org/10.1364/ol.37.001685 Such waveguides can tightly confine light to a sub-square micrometer cross-section and can be made from highly nonlinear materials. Moreover, lithographically defined on the nanoscale, integrated waveguides offer new opportunities for dispersion engineering via the modification of their dimensions or through subwavelength patterning along the propagation direction (e.g., periodic poling or photonic bandgap engineering). Through advances in thin film technology, a wide variety of materials is available, including materials with transparency in the terahertz, mid-infrared, visible, or ultraviolet wavelength domains. Waveguides with cubic nonlinearity have given rise to ultrabroadband supercontinua8–128. A. R. Johnson, A. S. Mayer, A. Klenner, K. Luke, E. S. Lamb, M. R. E. Lamont, C. Joshi, Y. Okawachi, F. W. Wise, M. Lipson, U. Keller, and A. L. Gaeta, “Octave-spanning coherent supercontinuum generation in a silicon nitride waveguide,” Opt. Lett. 40(21), 5117–5120 (2015). https://doi.org/10.1364/ol.40.0051179. H. Guo, C. Herkommer, A. Billat, D. Grassani, C. Zhang, M. H. P. Pfeiffer, W. Weng, C.-S. Brès, and T. J. Kippenberg, “Mid-infrared frequency comb via coherent dispersive wave generation in silicon nitride nanophotonic waveguides,” Nat. Photonics 12(6), 330–335 (2018). https://doi.org/10.1038/s41566-018-0144-110. M. Yu, B. Desiatov, Y. Okawachi, A. L. Gaeta, and M. Lončar, “Coherent two-octave-spanning supercontinuum generation in lithium-niobate waveguides,” Opt. Lett. 44(5), 1222–1225 (2019). https://doi.org/10.1364/ol.44.00122211. Y. Okawachi, M. Yu, B. Desiatov, B. Y. Kim, T. Hansson, M. Lončar, A. L. Gaeta, A. L. Gaeta, and A. L. Gaeta, “Chip-based self-referencing using integrated lithium niobate waveguides,” Optica 7(6), 702–707 (2020). https://doi.org/10.1364/optica.39236312. J. Rutledge, A. Catanese, D. D. Hickstein, S. A. Diddams, T. K. Allison, A. S. Kowligy, S. A. Diddams, T. K. Allison, T. K. Allison, A. S. Kowligy, and A. S. Kowligy, “Broadband ultraviolet-visible frequency combs from cascaded high-harmonic generation in quasi-phase-matched waveguides,” J. Opt. Soc. Am. B 38(8), 2252–2260 (2021). https://doi.org/10.1364/josab.427086 with simultaneous third-harmonic or triple-sum frequency generation (TFG).1313. E. Obrzud, V. Brasch, T. Voumard, A. Stroganov, M. Geiselmann, F. Wildi, F. Pepe, S. Lecomte, and T. Herr, “Visible blue-to-red 10 GHz frequency comb via on-chip triple-sum frequency generation,” Opt. Lett. 44(21), 5290–5293 (2019). https://doi.org/10.1364/ol.44.005290 Waveguides with quadratic and cubic nonlinearities have given rise to high-harmonic generation and simultaneous spectral broadening.14–1614. S. Wabnitz and V. V. Kozlov, “Harmonic and supercontinuum generation in quadratic and cubic nonlinear optical media,” J. Opt. Soc. Am. B 27(9), 1707–1711 (2010). https://doi.org/10.1364/josab.27.00170715. F. Baronio, M. Conforti, C. De Angelis, D. Modotto, S. Wabnitz, M. Andreana, A. Tonello, P. Leproux, and V. Couderc, “Second and third order susceptibilities mixing for supercontinuum generation and shaping,” Opt. Fiber Technol. 18(5), 283–289 (2012). https://doi.org/10.1016/j.yofte.2012.07.00116. D. D. Hickstein, D. R. Carlson, A. Kowligy, M. Kirchner, S. R. Domingue, N. Nader, H. Timmers, A. Lind, G. G. Ycas, M. M. Murnane, H. C. Kapteyn, S. B. Papp, and S. A. Diddams, “High-harmonic generation in periodically poled waveguides,” Optica 4(12), 1538 (2017). https://doi.org/10.1364/optica.4.001538 Quadratic nonlinearities have also been observed in waveguides made from centrosymmetric materials, where they can originate from surface-induced symmetry breaking or photo-induced periodic redistribution of charges.17,1817. X. Lu, G. Moille, A. Rao, D. A. Westly, and K. Srinivasan, “Efficient photoinduced second-harmonic generation in silicon nitride photonics,” Nat. Photonics 15(2), 131–136 (2021). https://doi.org/10.1038/s41566-020-00708-418. E. Nitiss, J. Hu, A. Stroganov, and C.-S. Brès, “Optically reconfigurable quasi-phase-matching in silicon nitride microresonators,” Nat. Photonics 16(2), 134–141 (2022). https://doi.org/10.1038/s41566-021-00925-5 Designing an integrated waveguide usually relies on numeric simulation of the underlying nonlinear optical processes. However, efficiently simulating all possible cascaded nonlinear optical processes in integrated waveguides is non-trivial as it requires a broadband, densely sampled, uninterrupted frequency interval. Additionally, cascaded nonlinearities can lead to the formation of numeric aliases (i.e., incorrect simulation results) due to the cyclic nature of the discrete Fourier transform, thus necessitating the implementation of anti-aliasing strategies for efficient and accurate simulation.Nonlinear propagation phenomena can be described via the generalized nonlinear Schrödinger equation (GNLSE), which remains valid even for extremely short pulses in the single cycle regime.1,19,201. J. M. Dudley, G. Genty, and S. Coen, “Supercontinuum generation in photonic crystal fiber,” Rev. Mod. Phys. 78(4), 1135–1184 (2006). https://doi.org/10.1103/revmodphys.78.113519. P. V. Mamyshev and S. V. Chernikov, “Ultrashort-pulse propagation in optical fibers,” Opt. Lett. 15(19), 1076–1078 (1990). https://doi.org/10.1364/ol.15.00107620. T. Brabec and F. Krausz, “Nonlinear optical pulse propagation in the single-cycle regime,” Phys. Rev. Lett. 78(17), 3282–3285 (1997). https://doi.org/10.1103/physrevlett.78.3282 Besides self-phase modulation (SPM), additional effects such as Raman scattering, self-steepening, frequency dependence of nonlinear parameters and refractive index can readily be included. The GNLSE has also been adapted in recent years to incorporate additional cubic and quadratic nonlinear effects such as SFG or TFG.15,21–3015. F. Baronio, M. Conforti, C. De Angelis, D. Modotto, S. Wabnitz, M. Andreana, A. Tonello, P. Leproux, and V. Couderc, “Second and third order susceptibilities mixing for supercontinuum generation and shaping,” Opt. Fiber Technol. 18(5), 283–289 (2012). https://doi.org/10.1016/j.yofte.2012.07.00121. M. Kolesik, J. V. Moloney, and M. Mlejnek, “Unidirectional optical pulse propagation equation,” Phys. Rev. Lett. 89(28), 283902 (2002). https://doi.org/10.1103/physrevlett.89.28390222. M. Kolesik and J. V. Moloney, “Nonlinear optical pulse propagation simulation: From Maxwell’s to unidirectional equations,” Phys. Rev. E 70(3), 036604 (2004). https://doi.org/10.1103/PhysRevE.70.03660423. G. Genty, P. Kinsler, B. Kibler, and J. M. Dudley, “Nonlinear envelope equation modeling of sub-cycle dynamics and harmonic generation in nonlinear waveguides,” Opt. Express 15(9), 5382–5387 (2007). https://doi.org/10.1364/oe.15.00538224. M. Conforti, F. Baronio, and C. De Angelis, “Nonlinear envelope equation for broadband optical pulses in quadratic media,” Phys. Rev. A 81(5), 053841 (2010). https://doi.org/10.1103/physreva.81.05384125. M. Conforti, F. Baronio, and C. De Angelis, “Ultrabroadband optical phenomena in quadratic nonlinear media,” IEEE Photonics J. 2(4), 600–610 (2010). https://doi.org/10.1109/jphot.2010.205153726. C. R. Phillips, C. Langrock, J. S. Pelc, M. M. Fejer, I. Hartl, and M. E. Fermann, “Supercontinuum generation in quasi-phasematched waveguides,” Opt. Express 19(20), 18754–18773 (2011). https://doi.org/10.1364/oe.19.01875427. M. Conforti, A. Marini, T. X. Tran, D. Faccio, and F. Biancalana, “Interaction between optical fields and their conjugates in nonlinear media,” Opt. Express 21(25), 31239 (2013). https://doi.org/10.1364/oe.21.03123928. M. Bache, “The nonlinear analytical envelope equation in quadratic nonlinear crystals,” arXiv:1603.00188 [physics] (2016).29. O. Melchert and A. Demircan, “py-fmas: A python package for ultrashort optical pulse propagation in terms of forward models for the analytic signal,” Computer Phys. Commun. 273, 108257 (2021). https://doi.org/10.1016/j.cpc.2021.10825730. “pyNLO, Nonlinear optics modeling for Python.” Github. https://github.com/pyNLO/PyNLO

Here, we present a numeric tool for the simulation of light propagation in nonlinear media and waveguides. It implements an extended GNLSE accounting for cubic and quadratic nonlinear interactions, as well as Raman scattering and self-steepening. To avoid computational artifacts in the multi-octave nonlinear simulation, anti-aliasing strategies are implemented. Modeling of frequency- or spatially-dependent (along the propagation direction) parameters is readily possible, for instance for frequency conversion in poled waveguides or materials with spatially-varying dispersion. A dedicated, fifth-order Runge–Kutta solver with adaptive step-size has been implemented for a high rate of convergence and short computation times even on generic desktop computers. To validate the numeric tool, broadband frequency conversion experiments in photonic integrated waveguides with quadratic and cubic nonlinearities are performed and shown to agree well with the numeric results. The numeric simulation tool is made available as an open-source Python package (pychi).

II. PROPAGATION MODEL

Section:

ChooseTop of pageABSTRACTI. INTRODUCTIONII. PROPAGATION MODEL <<III. COMPUTATIONAL SCHEMEIV. COMPARISON WITH EXPER...V. CONCLUSIONSUPPLEMENTARY MATERIALREFERENCESFollowing the approach in previous work,2727. M. Conforti, A. Marini, T. X. Tran, D. Faccio, and F. Biancalana, “Interaction between optical fields and their conjugates in nonlinear media,” Opt. Express 21(25), 31239 (2013). https://doi.org/10.1364/oe.21.031239 a temporal field envelope is defined as a function of propagation coordinate z and time t,A(z,t)=E(z,t)+iH[E(z,t)]e−iω0t+iβ(ω)z,(1)where E(z, t) is the electric field, H denotes the Hilbert transform, ω0 is the center frequency of the simulation, and β(ω) is the frequency-dependent wave vector in the nonlinear medium. Note that Fourier transforms of quantities are denoted with the corresponding lowercase symbols, e.g., F[A(z,t)](ω)=a(z,ω). The evolution of the spectral field envelope along the propagation direction z is given by∂za(z,ω)+iβ(ω+ω0)a(z,ω)=N[a(z,ω)],(2)where N[a(z,ω)] is the general nonlinear operator describing wave-mixing interactions. Keeping all interactions up to cubic order yieldsN[a(z,ω)]=−i(ω+ω0)4n(ω+ω0)cpNL(2)(z,ω)+pNL(3)(z,ω),wherePNL(2)(z,t)=χ(2)A(z,t)2eiω0t(SFG)+2|A(z,t)|2e−iω0t(DFG),andPNL(3)(z,t)=χ(3)23A(z,t)|A(z,t)|2(SPM)+3A*(z,t)|A(z,t)|2e−2iω0t(CKT)+A3(z,t)e2iω0t(TFG).The terms responsible for sum-frequency (SFG), difference-frequency (DFG), and triple-sum frequency generation (TFG), as well as self-phase modulation (SPM) and the conjugated Kerr term (CKT), have been identified; n(ω) is the effective refractive index (including material and geometric dispersion effects), and χ(2) and χ(3) are the quadratic and cubic nonlinear coefficients, respectively. These equations accounts for the frequency-dependence of the nonlinear operator (which is also responsible for self-steepening) and are here slightly modified to also include Raman scattering via the introduction of a frequency-dependent susceptibility,31,3231. K. J. Blow and D. Wood, “Theoretical description of transient stimulated Raman scattering in optical fibers,” IEEE J. Quantum Electron. 25(12), 2665–2673 (1989). https://doi.org/10.1109/3.4065532. R. H. Stolen, W. J. Tomlinson, H. A. Haus, and J. P. Gordon, “Raman response function of silica-core fibers,” J. Opt. Soc. Am. B 6(6), 1159–1166 (1989). https://doi.org/10.1364/josab.6.001159χ(3)(ω′,ω″,ω′′′)≔χ(3)g(ω′−ω″),(3)where ω′, ω″, and ω‴ are the angular frequencies of the mixing fields driving the polarization, withG(t)=αδ(t)+(1−α)τ12+τ22τ1τ22e−t/τ2⁡sin(t/τ1),(4)where α, τ1, and τ2 are material-dependent quantities (typically, α ≈ 0.7 to 0.82, τ1 ≈ 0.0122 ps, and τ2 ≈ 0.032 ps for silica fibers31,3231. K. J. Blow and D. Wood, “Theoretical description of transient stimulated Raman scattering in optical fibers,” IEEE J. Quantum Electron. 25(12), 2665–2673 (1989). https://doi.org/10.1109/3.4065532. R. H. Stolen, W. J. Tomlinson, H. A. Haus, and J. P. Gordon, “Raman response function of silica-core fibers,” J. Opt. Soc. Am. B 6(6), 1159–1166 (1989). https://doi.org/10.1364/josab.6.001159).

III. COMPUTATIONAL SCHEME

Section:

ChooseTop of pageABSTRACTI. INTRODUCTIONII. PROPAGATION MODELIII. COMPUTATIONAL SCHEME <<IV. COMPARISON WITH EXPER...V. CONCLUSIONSUPPLEMENTARY MATERIALREFERENCESThrough the SFG, DFG, CKT, and TFG nonlinear interaction terms, light with frequency components outside of the simulation range can be generated. These frequency components will get folded back into the simulation domain due to the cyclic nature of the discrete Fourier transform. Untreated, this behavior leads to spurious light aliases in the simulation window. Simply extending the simulation window is not a viable strategy as the computational costs will increase accordingly, and further cascaded nonlinear interactions may again lead to the appearance of aliases. To avoid unphysical simulation results, an anti-aliasing scheme is implemented, as exemplified by the quadratic nonlinearity in Fig. 1. Such techniques have been largely applied in the fields of nonlinear electronics and fluid mechanics and can be equivalently used here. At each step, the field is Fourier transformed and padded with zeros in the frequency domain, effectively augmenting the span of the simulation window. If the original simulation bandwidth is Δω, the padding needs to extend the frequency range by a factor N+12+(N−1)ω0Δω, where N is the nonlinearity’s order (note that in the case of mixed nonlinearities, the nonlinear terms are computed separately). The padded field is transformed back into the time domain, and the nonlinear terms are computed, with the resulting frequencies all comprised within the extended simulation window.33,3433. S. A. Orszag, “On the elimination of aliasing in finite-difference schemes by filtering high-wavenumber components,” J. Atmos. Sci. 28(6), 1074 (1971). https://doi.org/10.1175/1520-0469(1971)028<1074:OTEOAI>2.0.CO;234. S. Derevyanko, “The (n + 1)/2 rule for dealiasing in the split-step Fourier methods for n-wave interactions,” IEEE Photonics Technol. Lett. 20(23), 1929–1931 (2008). https://doi.org/10.1109/lpt.2008.2005420 After Fourier transforming the time-domain nonlinear terms, only the frequency components contributing to the original simulation window are retained, providing an alias-free computation. Notably, this scheme preserves the total energy of the pulse as the discarded spectral components are not allowed to build up and thus do not draw energy from the pulse. Further details about the anti-aliasing technique such as comparison against propagation without anti-aliasing scheme, proof of energy conservation, and the choice of padding sizes are given in the supplementary material.In general, nonlinear interactions are computationally expensive to simulate and thus require some degree of numerical optimization in order to achieve fast computation times, particularly in models where multiple nonlinear interactions are simulated. To this end, the efficient fast Fourier transform library pyFFTW and the optimization library numba have been used. A powerful numeric solver that could be used here is the fourth-order Runge–Kutta RK4IP scheme.3535. H.. Johan, “A fourth-order Runge–Kutta in the interaction picture method for simulating supercontinuum generation in optical fibers,” J. Lightwave Technol. 25(12), 3770–3775 (2007). https://doi.org/10.1109/JLT.2007.909373 It leverages the so-called interaction picture of the propagation (corresponding to a change of reference frame analogous to the quantum mechanical interaction picture) in order to reduce the number of computational operations per step. However, it is not an embedded method, i.e., it does not directly provide an error estimation for an efficient adaptive step-size control. Strategies such as step-doubling have been employed to provide error estimation. Although this significantly augments the computational cost per step, the reduction in the required number of steps results in faster overall computation times.36,3736. O. V. Sinkin, R. Holzlohner, J. Zweck, and C. R. Menyuk, “Optimization of the split-step Fourier method in modeling optical-fiber communications systems,” J. Lightwave Technol. 21(1), 61–68 (2003). https://doi.org/10.1109/jlt.2003.80862837. A. M. Heidt, “Efficient adaptive step size method for the simulation of supercontinuum generation in optical fibers,” J. Lightwave Technol. 27(18), 3984–3991 (2009). https://doi.org/10.1109/JLT.2009.2021538 A later adaptation called the ERK4IP3838. S. Balac and A. Fernandez, “Comparison of adaptive step-size control strategies for solving the Generalised Non-Linear Schrodinger Equation in optics by the Interaction Picture method,” Technical report hal-00740771, 2012, https://hal.science/hal-00740771. provided an embedded version of the RK4IP requiring no extra operations for error estimation and was also tested in this work. In our implementation, using this solver yielded a 3–4 times reduction of computation time over the step-doubling RK4IP for SPM-based propagation. When other nonlinear processes are taken into account, we found that even higher order schemes perform better by providing a higher rate of convergence. Consequently, we implemented a new fifth-order Runge–Kutta solver based on the Dormand-Prince method3939. J. R. Dormand and P. J. Prince, “A family of embedded Runge-Kutta formulae,” J. Comput. Appl. Math. 6(1), 19–26 (1980). https://doi.org/10.1016/0771-050x(80)90013-3 adapted in the interaction picture, requiring fewer dispersion operations.With the field in the frequency domain, the algorithm is the following:k1=e−iΔzβNa0ork1=e−iΔzβk7m−1,(6)k2=e−45iΔzβNe45iΔzβ(aI+Δzα21k1),(7)k3=e−710iΔzβNe710iΔzβaI+Δz⋅∑j=12α3jkj,(8)k4=e−15iΔzβNe15iΔzβaI+Δz⋅∑j=13α4jkj,(9)k5=e−19iΔzβNe19iΔzβaI+Δz⋅∑j=14α5jkj,(10)k6=NaI+Δz⋅∑j=15α6jkj,(11)am+1=aI+Δz⋅∑j=16α7jkj,(12)bm+1=aI+Δz⋅∑j=17α8jkj,(14)ε=‖am+1−bm+1‖‖am+1‖,(15)where Δz is the step-size, am is the electric field envelope at step m, k7m−1 is the previously computed value of k7, and the αij coefficients are given in the Butcher tableau, Table I.

TABLE I. Butcher tableau for the solver implemented in this work.

α2j1/5α3j3/409/40α4j44/45−56/1532/9α5j19  372/6561−25 360/218764  448/6561−212/729α6j9017/3168−355/3346  732/524749/176−5103/18  656α7j35/3840500/1113125/192−2187/678411/84α8j5179/57  60007571/16  695393/640−92 097/3 39  200187/21001/40j1234567

This algorithm is of order five in the step-size and provides a relative error ɛ for step-size control. Additionally, it has the first-same-as-last property, i.e., the evaluation of k1 does not require a nonlinear operation after the first step. Compared to the ERK4IP, the proposed scheme requires six nonlinear operator evaluations against four, but is one order higher. It has been found that the total number of nonlinear operator evaluations required to achieve comparable simulation results between the two solvers does not differ. However, the higher rate of convergence of the order five solver leads to improved error estimation and step-size control, which has empirically proven useful for the dynamics of mixed and cascaded nonlinearities.

IV. COMPARISON WITH EXPERIMENTS

Section:

ChooseTop of pageABSTRACTI. INTRODUCTIONII. PROPAGATION MODELIII. COMPUTATIONAL SCHEMEIV. COMPARISON WITH EXPER... <<V. CONCLUSIONSUPPLEMENTARY MATERIALREFERENCESTo verify the numeric simulation, comparisons with experimental results have been performed. In the first experiment, a lithium niobate waveguide (mode area ∼1 µm2, length 6 mm; see the supplementary material for more details), exhibiting both quadratic and cubic nonlinearities, has been pumped with a 100 MHz Erbium-based mode-locked laser. The pump pulses of 80 fs duration and central wavelength of 1560 nm are coupled to the waveguide via a lensed fiber. The resulting optical spectrum obtained with ∼10 pJ on-chip pulse energy is then recorded via an optical spectrum analyzer (see Fig. 2). To reproduce the experimentally obtained spectrum, material parameters and effective refractive indices (simulated by a separate finite element simulation) are used as input parameters of the simulation and a sech2 pulse intensity profile is assumed. The simulated spectrum is shown in Fig. 2 (computed in ∼300 s), as well as a simulated spectrum based on SPM only. Experimental and simulated results are in excellent agreement; in particular, the width, position, and amplitude of the SFG are accurately reproduced. The finer structure apparent in the simulated SFG spectrum results from the phase mismatch between the pump and the SFG, leading to fringes not resolved by the optical spectrum analyzer. This comparison validates the two different nonlinear operators responsible for quadratic interactions, namely, SFG and DFG generation as well as SPM through cubic nonlinearity.Similarly, a silicon nitride waveguide exhibiting a purely cubic intrinsic material nonlinearity has been tested (mode area ∼1 µm2, length 5 mm; see the supplementary material for more details). The experimentally observed spectrum obtained with ∼100 pJ on-chip pulse energy reveals strong, SPM-dominated spectral broadening and TFG (see Fig. 3). Again, material parameters and effective refractive indices of the waveguide are used as input parameters to the numeric simulation. The simulated spectrum (computed in ∼1100 s), along with an SPM-only simulated spectrum for comparison, is shown in Fig. 3. Both simulated spectra reproduce the SPM-based broadening, but only the complete model captures the TFG part of the spectrum. The overall position and envelope of the TFG are well reproduced, while deviations are visible in the finer structure of the TFG spectrum. We attribute those deviations to increasing uncertainty in the material data toward shorter wavelength (impacting the waveguide dispersion), as well as physical effects not included in the simulation such as surface roughness induced scattering and interaction between higher order transverse waveguide modes (e.g., a sharp peak experimentally observed at around 450 nm). These results validate the remaining nonlinear operators in the propagation model, namely, SPM, TFG, and the CKT. We note that the solver has also been validated against previous numeric results1