1. Introduction

Modern optical waveguides offer a growing range of capabilities to generate and manipulate ultrashort pulses. A key enabling material platform is thin-film lithium niobate (TFLN), which has been particularly prominent since the demonstration of ultralow-loss etched waveguides [1] implemented on LiNbO$_3$-on-insulator wafers [2]. Compared to conventional $\chi ^{(2)}$ waveguides such as reverse-proton-exchanged (RPE) waveguides [3], these new TFLN waveguides offer tight mode confinement due to their high index contrast. By combining this property with quasi-phase-matching (QPM), ultra-high normalized efficiencies have been demonstrated [4]. The ability to achieve tight mode confinement enables dispersion engineering, which has been used to match the group velocities of the interacting waves and thereby obtain ultra-wide acceptance bandwidths [5–7]. By simultaneously confining the mode area to the scale of a wavelength (transverse confinement) and utilizing dispersion engineering to maintain temporal overlap between few-optical-cycle pulses over an extended propagation distance, it has become possible to scale down the pulse energy requirements of highly nonlinear optical interactions from the nJ regime to the pJ or even fJ regime.

These advances are of high interest for applications. For example, deploying these new capabilities for high repetition rate optical frequency comb generation [8–10] could benefit applications such as molecular fingerprint spectroscopy [11], precision metrology [12], and optical coherence tomography [13], and supercontinuum generation [14]. Scaling combs to high repetition rates leads to more compact systems, more easily resolvable comb lines, and higher signal-to-noise ratio in dual-comb spectroscopy due to the higher power per comb line [15]. Second-order nonlinear waveguides are particularly appealing for comb generation since they can combine low-energy supercontinuum generation and harmonic generation for f-2f interferometry in the same device [16–22]. Mid-infrared generation can be supported as well by using TFLN devices with sapphire as the cladding material [23]. More generally, the tight bend radii and range of linear, nonlinear, electro-optic functionalities available in these waveguides are extremely promising for novel chip-scale optical systems.

Beyond bandwidth and pulse energy scaling, recent progress in second-order nonlinear photonics has enabled qualitatively new operating regimes for generating, amplifying, or compressing ultrashort pulses. Examples include multi-color solitons [24,25], octave-spanning optical parametric amplifiers [26–28], new approaches to supercontinuum generation [29] (including process that resemble or focus on the generation of high-order harmonics [20,22,30]), and complicated multi-wave oscillators [28]. These advances motivate the development of pulse propagation models to meet two important goals: (i) to capture the details of experimental results, and (ii) as tools to understand and/or design devices. Developing a modeling framework that achieves both of these goals simultaneously is one of the main goals of this work.

Ultra-broadband nonlinear interactions in $\chi ^{(2)}$ waveguides involve several simultaneous processes, including the generation of multiple harmonics of the input pulse. In particular, pulses with octave-spanning spectra can be generated, leading to overlap between the different harmonic components within the waveguide. As a result, the conventional second-harmonic generation (SHG) coupled wave equations are not sufficiently general to capture such $\chi ^{(2)}$ supercontinuum generation experiments. There are at least two possible failure mechanisms of these models: the generation of additional harmonics, and the (possible) breakdown of their validity due to the envelopes having overlapping spectra. Single-spectral-envelope (SSE) models were developed to address such concerns by incorporating all $\chi ^{(2)}$-based driving terms into the propagation equation of a single electric field envelope [31,32]. These models could explain observed SCG results once generalized to include additional effects [33]. However, in these models all nonlinear effects and all harmonics are included at once, making it difficult to isolate individual effects. Therefore, despite the power of these SSE models, we have found it is often preferable to use generalized multi-spectral-envelope (MSE) models instead of, or in parallel to, the SSE models. Such MSE models, when implemented in a sufficiently general way, give greater flexibility in switching on or off different effects to build intuition, iterate on a design, and isolate dominant effects to inform development of accurate reduced models. Yet the question remains of whether these models are sufficiently general, or if important information is lost.

Here we present a common framework incorporating both types of models, and show the utility of both approaches in explaining recent ultra-broadband experiments. We include a discussion of their equivalence, stability, and frequency comb aspects such as the generation and extraction of harmonic combs. A general overarching question which we address is: when is it best to use one or the other numerical approach? We give guidelines for this choice which should be sufficient to cover most cases of practical interest. We will show that in many contexts MSE models can have advantages in terms of efficient device design, understanding essential device physics, computational efficiency, and flexibility in defining the set of nonlinear interactions involved and the associated wave coupling coefficients. Our primary focus is on $\chi ^{(2)}$ effects, but we show the corresponding $\chi ^{(3)}$ terms as well.

In section 2 we summarize the mathematical description of nonlinear mixing of waveguide modes. In section 3 we develop the single-spectral-envelope and multi-spectral-envelope numerical models, denoted SSE and MSE respectively. In sections 4 and 5 we present worked examples based on a periodically poled TFLN-on-sapphire waveguide that was used recently for supercontinuum generation experiments [34]. In section 6 we discuss the results and relevant considerations when choosing an appropriate numerical model in different situations. Our work will help enable the continued development and understanding of broadband $\chi ^{(2)}$ devices and their deployment in chip-scale photonic systems.

2. Theoretical framework

Nonlinear interactions in optical waveguides are generally modeled within the framework of perturbative nonlinear optics [35]. Consequently, the electromagnetic field evolves according to linear-optical propagation effects plus a series of perturbation terms. Although linear optics permits propagation in either direction of the waveguide, it is computationally much simpler to model only one direction. A framework for such unidirectional propagation models was presented in [36], and here we begin with a minimal description of this formalism to introduce relevant notation. Although this section mainly summarizes prior work, we include it to have a consistent and self-contained theoretical framework.

2.1 Nonlinear propagation of waveguide modes

It is convenient to work primarily in the frequency domain since dispersion effects can then be treated exactly. We denote frequency domain quantities by tilde and the Fourier transform by $\mathscr{F}[]$:

Note, we often give quantities in terms of angular optical frequency $\omega =2\pi \nu$ for reasons of compactness and convention, but the simulations themselves generally use non-angular frequencies for the Fourier transform and its inverse. The inverse Fourier transform of a function expressed in terms of angular frequency $\omega$ is given by $g(t) = \int g(\omega ) \exp (i \omega t) d(\omega / (2\pi ))$. In the following, vectors are indicated by bold-face symbols

Nonlinear effects are governed by the nonlinear polarization $\mathbfcal{P}_{NL}$, which contains various nonlinear susceptibility terms. In general each term involves a convolution-type integral over the interacting frequencies, weighted by the frequency-dependent susceptibility. In the simplified case of instantaneous nonlinear responses the integral is a true convolution, in which case it is possible to use a much simpler time-domain expression. We consider dielectric media which are nonlinear, inhomogeneous, anisotropic, and non-gyrotropic. In this case the constitutive relations are given by

Equations (2) and (3) can be inserted into Maxwell’s equations in order to solve for nonlinear wave propagation. However, the problem becomes significantly more tractable if we assume propagation in a small set of bound waveguide modes. In a waveguide whose linear properties vary with transverse coordinates $x$ and $y$ but not with propagation coordinate $z$, the electric field can be expressed as a summation over the waveguide modes:

The mode profiles and corresponding dispersion relations can be found through a variety of numerical techniques [37]. The orthogonality condition for waveguide modes $n$ and $m$ at a particular frequency can be written as [38]

A propagation equation for each of the modal fields $E_n$ can be determined by applying the reciprocity relations for Maxwell’s equations and orthogonality of the waveguide modes. This type of procedure has been applied for both free space and waveguide configurations [36]. By following that procedure, assuming forward-propagating modes, and converting from mode envelopes $E_n$ to mode amplitudes $A_n$ we obtain the following propagation equation for each waveguide mode:

Neglecting counter-propagating waveguide modes when calculating $\tilde {P}_n$ is appropriate for a large portion of practical experiments. The assumption is particularly appropriate for ultra-broadband nonlinear optics since interactions between counter-propagating waves exhibit a very large group velocity mismatch (GVM) that tends to limit the bandwidth of the counter-propagating waves [39].

2.2 Polarization contributions and coupling coefficients

The frequency dependence of the mode profiles and nonlinear susceptibilities are embedded in $\tilde {P}_n$. There are several types of terms, and each one involves an integral over frequency and summation over waveguide modes.

For $\chi ^{(2)}$ there are SFG-type and DFG-type terms, involving $\chi ^{(2)}_{ijk}(\omega _1; \omega _2, \omega _3)$ with $\omega _1=\omega _2+\omega _3$. Here we assume positive frequencies, with the first frequency argument corresponding to the high frequency wave, and the polarization subscripts have the same ordering as the frequency arguments. For polarization subscripts we will use $(i,j,k,l)$; for mode index subscripts we use $(n,p,q,r)$. Therefore, the SFG contributions have the form $\mathcal{P}_i(\omega _1) \sim \chi ^{(2)}_{ijk}(\omega _1; \omega _2, \omega _1-\omega _2)\mathcal{E}_j(\omega _2)\mathcal{E}_k(\omega _1-\omega _2)$ while DFG contributions have the form $\mathcal{P}_j(\omega _2) \sim \chi ^{(2)}_{ijk}(\omega _1; \omega _2, \omega _1-\omega _2)\mathcal{E}_i(\omega _1)\mathcal{E}_k^*(\omega _1-\omega _2)$.

For $\chi ^{(3)}$ there are SFG, DFG, and four-wave-mixing (FWM) terms. In the SFG and DFG cases the interacting frequencies have the structure $\omega _1=\omega _2+\omega _3+\omega _4$, while the FWM frequencies have structure $\omega _1-\omega _2=\omega _3-\omega _4$. On this basis we use the position of the semicolon to distinguish these cases: the notation $\chi _{ijkl}^{(3)}(\omega _1; \omega _2, \omega _3, \omega _4)$ is used for the SFG/DFG, and $\chi ^{(3)}_{ijkl}(\omega _1, \omega _2; \omega _3, \omega _4)$ is used for FWM. The FWM coefficients include self-phase modulation (SPM) via $\chi ^{(3)}_{iiii}(\omega _1,\omega _1;\omega _1,\omega _1)$ and cross-phase modulation (XPM) via $\chi ^{(3)}_{iijj}(\omega _1,\omega _1;\omega _2,\omega _2)$.

All the polarization terms involve convolution-type integrals weighted by effective mixing coefficients $X_{npq}$ and $X_{npqr}$ which are defined below. The $\chi ^{(2)}$ SFG terms are

Polarization terms which can never generate positive Fourier components have been neglected since these are suppressed by the $u(\omega )$ factor in Eq. (7). Basically, such terms are those where all field products in the integrand are conjugated (since terms that are conjugated have their frequencies subtracted from the total).

The $\chi ^{(2)}$ effective nonlinear coefficients $X_{npq}$, which are frequency-dependent scalars, account for the modal overlap between the three frequency components involved:

Here we allow $\chi ^{(2)}$ to vary in the propagation coordinate $z$ to model QPM (discussed further in section 3.5), and also in coordinates $x$ and $y$ in order to account for the transverse properties of the waveguide. It is worthwhile to note that the conjugate relations between the mode profiles in this equation match the photon-energy conservation $\hbar {}\omega =\hbar {}\omega ' + \hbar {}(\omega -\omega ')$. The coefficient $X_{npq}(\omega ; \omega ', \omega -\omega ')$ will apply to SFG into field $\tilde {A}_n(\omega )$ driven by the product of fields $\tilde {A}_p(\omega ')$ and $\tilde {A}_q(\omega -\omega ')$. The corresponding DFG processes into fields $\tilde {A}_p(\omega ')$ and $\tilde {A}_q(\omega - \omega ')$ will be driven by the conjugate of the same coefficient, i.e. $X_{npq}(\omega ;\omega ',\omega -\omega ')^*$.

Analogously to the definition of $X_{npq}$, the $\chi ^{(3)}$ effective nonlinear coefficients for SFG/DFG are

Here $\chi ^{(3)}$ is assumed $z$-independent but the full transverse profile of $\chi ^{(3)}$ can be included (accounting for the presence of different materials in the waveguide, if applicable).

The total polarization can be expressed as a summation of the various terms defined above:

It is worthwhile to note that the dimensions of $X_{npq}$ and $X_{npqr}$ depend on how one chooses to scale the mode profiles $\mathbf {F}_n$, and this choice is arbitrary since a frequency-dependent scale factor applied to $\mathbf {F}_n(x,y,\omega )$ can be absorbed into $\tilde {A}_n(z,\omega )$. Hence, $g_n(\omega )$ can be chosen so as to minimize the variation of the effective nonlinear coefficients with respect to frequency. This is useful because when the dependence of the $X$ coefficients on the integrand frequencies can be neglected, Eqs. (14)–(16) are simply convolutions and can be expressed as the Fourier transform of a time-domain product, e.g. $\tilde {P}_{npq}^{(SFG)}\sim X_{npq}\mathscr{F}[A_p(t)A_q(t)]$.

The third-order nonlinear terms implicitly include important non-instantaneous nonlinear process such as stimulated Raman scattering (SRS). SRS can be described by a mixed time-frequency picture: the optical intensity exerts a force on the crystal (time domain), which excites phonons according to the material’s Raman transfer function $\chi _R(\Omega )$ (in principle a tensor, but we neglect that aspect here). This frequency dependence can be included while still preserving a computationally-convenient Fourier transform approach to calculating the nonlinear polarization [40] by approximating the Raman contribution to the FWM coefficients as $X_{npqr}^{(FWM)}(\omega,\omega ';\omega '',\omega ''+\omega '-\omega )\approx X_R(\omega -\omega ')$:

2.3 Frequency-shifted envelopes

While intuitive, the $A_n$ envelope choice is computationally inconvenient since the field rapidly accumulates phase and delay. To avoid this, we can define frequency-shifted envelopes and introduce a moving coordinate system. With the standard shifted time coordinate $t'=t-z/{v_{\mathrm {ref}}}$ and $z'=z$ the derivatives become $\partial /\partial z = \partial /\partial z' - (1/{v_{\mathrm {ref}}})\partial / \partial t'$, and $\partial ^n/\partial t^n = \partial ^n / \partial t'^n$. In the frequency domain, this implies

Therefore, if we define a frequency-shifted envelope as follows:

The term in square brackets in Eq. (20) is analogous to the carrier envelope phase shift of a pulse, but it is not equivalent to that quantity since both $\omega _0$ and ${v_{\mathrm {ref}}}$ can be chosen arbitrarily. Next, with the basic propagation equations set up, in section 3 we discuss different models and their trade-offs. The procedure of this subsection will be applied as needed to shift to the computationally-convenient $a_n$-type envelopes.

3. Numerical models

3.1 Overview

The standard approach to narrow-band nonlinear frequency conversion is to use one envelope for each wave or spectral region involved in the mixing process, for example a first harmonic (FH) and second harmonic (SH) [41]. The carrier frequencies of these envelopes are chosen as $\omega _2=2\omega _1$, with $\omega _1=\omega _c$ for an input wave of center frequency $\omega _c$. This choice allows resolution of the incident and generated pulse on as coarse a temporal grid as possible, which is computationally efficient (apart from this consideration, any value of $\omega _1$ could be chosen in principle). A generalization of this approach is to split the modal fields in Eq. (4) into a summation over the harmonic orders involved in the $\chi ^{(2)}$ interaction. We refer to this as a multi-spectral-envelope model.

When modeling ultra-broadband interactions, the optimality of using multiple envelopes representing harmonics of the input wave is not obvious. Once the pulse bandwidth becomes broad enough, the separate envelopes overlap in frequency and one might question whether information is lost in the simulation, or whether the approach is still computationally efficient. To address this issue, a single spectral envelope approach to $\chi ^{(2)}$ interactions was introduced and applied to analyze supercontinuum generation [31–33]. These models can be obtained by adapting generalized unidirectional propagation models (e.g. [36]) to the case of a $\chi ^{(2)}$ nonlinearity instead of $\chi ^{(3)}$. The mode amplitudes from section 2 correspond to a single spectral envelope, since the full field in each waveguide mode is described by only one envelope.

In this section we summarize the single-spectral/multi-spatial envelope approaches (denoted SSE) and introduce a generalized multi-spectral/multi-spatial envelope model (denoted MSE). Figure 1 shows a comparison between the two approaches by illustrating some of the interactions that can occur. Considering an input pulse with center frequency $\omega _c$, the interactions lead to new pulse components (harmonics) with center frequencies $N\omega _c$ for all integers $N\ge 0$; each of these harmonics (and the $N=0$ baseband term) can be associated with a different electric field envelope. Similarly, in a pulse train, each harmonic has a carrier envelope offset (CEO) frequency $f^{(N)}_{CEO} = Nf^{(1)}_{CEO}$, where $f^{(N)}_{CEO}$ denotes the CEO frequency of the $N^{th}$ harmonic [8–10]. Eventually, the envelopes become negligible for $N > N_{max}$, where integer $N_{max}$ is determined by the nonlinear dynamics. In an SSE model all harmonics are automatically generated and can mix together as long as the numerical frequency grid is sufficiently broad. These phenomena are naturally captured by MSE models by choosing $\omega _0=\omega _c$.

 

Fig. 1. Illustration of numerical models. (a) Electric field consisting of a pulse at frequency $\omega _0$, its harmonics, and an optical rectification component. Once the underlying fundamental wave’s bandwidth reaches an octave, the different harmonics overlap and the spectrum becomes phase dependent. The resulting changes in the computed single-shot spectra can be used to evaluate $f_{\mathrm {CEO}}$ beats for a pulse train. (b) The individual harmonic components (envelopes) corresponding to (a) that are included explicitly in an MSE model. Such MSE models must capture how each frequency can have source terms arising from different combinations of envelopes (see dashed arrows).

Download Full Size | PDF

When harmonic components exhibit spectral overlap the output spectrum depends on the carrier envelope phase of the input pulse (illustrated Fig. 1(a)). In an MSE model, only the harmonic components that are explicitly included can mix together, and each product between a set of envelopes can be assigned to the nonlinear polarization of a particular harmonic corresponding to the sum of the carrier frequencies involved. Two examples of this assignment are illustrated by the dashed lines in Fig. 1(b): the top dashed lines illustrate sum frequency generation between two components of the fundamental generating a second-harmonic component. The lower dashed lines illustrate sum frequency generation between a component of the fundamental and a component of the second harmonic generating a third-harmonic component.

3.2 Single spectral envelope model

The SSE model requires little extra derivation since section 2 already uses a single mode amplitude for each waveguide mode. Starting from the underlying model of Eq. (7), we first switch to the moving coordinate system and frequency-shifted envelopes $a_n$, but use a mode-number-independent carrier phase shift $\beta _0(\omega )$ instead of $\beta _n(\omega )$ for simplicity. To arrive at a minimal equation, we neglect $\chi ^{(3)}$ terms ($X_{npqr}=0$) and assume the $X_{npq}$ coefficients are non-dispersive. This latter assumption allows the nonlinear polarization contributions to be expressed in terms of Fourier transforms of time-domain products of the envelopes:

3.3 Introducing harmonic envelopes

A benefit of Eq. (22) is its generality: these SSE models can be applied directly to any incident field. However, it is easy to construct interactions where it is orders of magnitude slower to compute compared with traditional models. For example, in a basic SHG interaction between a pair of transform-limited picosecond pulses a conventional model can use a grid step size of order 1 ps, whereas with Eq. (22) a grid step size of order a few femtoseconds is necessary in order to resolve half the optical carrier frequency (assuming $\omega _0$ is selected mid-way between the FH and SH center frequencies). Another issue is that in an SSE model all sum frequency mixing terms are generated, which may require approaches such as an even finer temporal grid to resolve (or suppress) the small amount of light generated at higher harmonics without aliasing [42]. The issue is compounded when considering multiple QPM orders since highly phase-mismatched terms with short coherence lengths will be involved. The MSE model presented in this subsection offers greater flexibility in choosing the set of interactions included in the simulation.

Our approach is to express the electric field as a summation over a set of waveguide modes $n$ and harmonic orders $N$ with carrier frequencies $N\omega _0$, and to choose a subset of QPM grating orders for each mixing process. Neglected QPM orders can also be included as effective SPM and XPM terms by cascaded-$\chi ^{(2)}$ approximations [43] (see Supplement 1). This strategy transitions more gradually from the most basic interactions (one waveguide mode, one QPM order, and two spectral envelopes for QPM SHG) to ultra-broadband interactions where several waveguide modes, harmonic orders, and QPM orders can be involved. Rather than a single nonlinear term containing many individual mixing processes (such as Eq. (22)), here the meaning of each individual contribution to the total nonlinear polarization is preserved since each term is associated with mixing between a specific set of harmonic envelopes.

The mode amplitude for each waveguide mode is expressed as a sum over harmonic orders:

Inserting the $\tilde {A}_{n;N}$ from Eq. (23) into the nonlinear polarization $\tilde {P}_n$ contributions (Eqs. (10)–(13)) yields many terms. In order to arrive at an array of equations of the following form

In analogy to Eq. (17) there are two types of $\chi ^{(2)}$ polarization terms:

The integration limits in these definitions arise from using analytic signals. Analogous terms can be defined for the $\chi ^{(3)}$ mixing processes.

3.4 Multiple spectral envelopes model

Next, we introduce appropriate carrier frequencies and phases into the electric field envelopes, simplify the effective nonlinear susceptibilities, and re-arrange the summations to express the propagation equations in terms of energy-conserving sets of nonlinear mixing terms.

We first define harmonic mode envelopes $a_{n;N}$ using carrier frequencies $N\omega _0$, carrier propagation constants $\beta _0(N\omega _0)$, and assuming a reference group velocity ${v_{\mathrm {ref}}}$. These definitions are analogous to Eq. (20):

The array of propagation equations is given, in analogy to Eq. (21), by

Since the $\tilde {a}_{n;N}$ are defined with a frequency shift corresponding to harmonic order $N$, their frequency argument $\Omega$ corresponds to optical frequencies $\omega =\Omega +N\omega _0$.

What remains is to obtain a convenient and sufficiently accurate expression for the polarization arrays $\underline {p}_n$. As with Eq. (22), assuming that the nonlinear coefficients are non-dispersive enables all terms to be expressed in terms of Fourier transforms of time domain products. However, these coefficients will now take the form of multidimensional arrays since products of all combinations of envelopes give rise to nonlinear polarization terms.

3.5 Second-order nonlinear terms

As a first step towards determining an expression for $\underline {p}_n$, we simplify the effective susceptibilities by approximating them as independent of frequency within the integrands defined by Eqs. (26), (27). For QPM interactions, the grating is expanded in its Fourier series components [44]. We assume that changes in $\chi ^{(2)}$ are separable in transverse coordinate $(x,y)$ and propagation coordinate $z$, and that $\chi ^{(2)}(x,y,z)=\pm |\chi ^{(2)}(x,y,0)|$. With this assumption the $z$-dependence can be written as a Fourier series, i.e. a summation of the form $\sum _m \bar {d}_m \exp (im\phi _G(z))$ [44], and for a constant QPM duty cycle the normalized Fourier coefficients $\bar {d}_m$ are constant. Based on this Fourier series representation we define an array of $\chi ^{(2)}$ nonlinear coefficients according to

The phase mismatch terms are given by

The above definitions are sufficient to express the $\chi ^{(2)}$ contributions to the effective polarization array $\underline {\tilde {p}}_{n}$. However, it is particularly useful to re-arrange the indices so that the SFG-type and DFG-type terms share a single coefficient within the summand. To do so, some additional notational machinery is needed so that different parts of the summand can correspond to different waveguide indices. To capture this, let $\underline {\mathfrak {p}}$ denote a two-dimensional array whose first index corresponds to waveguide mode and second index to harmonic order. Further, multiplying a summand term by $\delta _{n;N}$ indicates that the term will be added to the $(n,N)$ entry of $\underline {\mathfrak {p}}$. With this notation, and using the permutation symmetries of ${\underline {\bar {X}}}{}_{npq}^{(m)}$, we find

For convenience we recap the notation: $a_{n;N}$ refers to waveguide index $n$, harmonic order $N$. In $({\underline {\bar {X}}}{}_{npq}^{(m)})_{P,Q}$ the $\{n,p,q\}$ refer to waveguide indices, $m$ refers to QPM grating order, and $\{P,Q\}$ refer to harmonic orders.

The system of equations can be further refined by keeping only the unique terms from the summation (multiplying them by appropriate scale factors based on whether terms were duplicated or not in the original summation). Note that the form of the phase mismatch in the ${\bar {X}}{}$ terms is $\exp (+i\Delta \beta z)$ while for the ${\bar {X}}{}^*$ terms it is $\exp (-i\Delta \beta z)$, which ensures the same structure as the usual coupled three-wave mixing equations. In general, the terms as written closely resemble the coupled three-wave mixing equations, with the first term representing SFG and the latter two terms representing the corresponding DFG process. For SHG-like cases ($p=q$ and $P=Q$) the second and third terms are the same and need only be calculated once (though of course multiplied by the appropriate factor in evaluating the sum). The formulation in Eq. (32) has several useful features:

3.6 Third-order nonlinear terms

For the $\chi ^{(3)}$ terms, the different types of contributions to the third-order susceptibility must be considered. Here we will assume that $\chi ^{(3)}$ is comprised of a “fast” electronic component $\chi _E$ and a “slow” Raman component $\chi _R$. A simple model that captures the effects of interest is to approximate the third-order susceptibility according to

We also define $\chi ^{(3)}$ carrier phase mismatch terms as

For SPM and XPM effects, these phase mismatch terms are zero. Using this form of $\Delta \beta$, we can define coefficients for the FWM-type terms (which includes SPM and XPM terms) as follows:

One could re-arrange these polarization terms to obtain a similar structure to Eq. (32), but since our focus is on $\chi ^{(2)}$ processes we do not write this out explicitly. The generalized coupled wave equations defined by Eqs. (29), (32), (38) capture all relevant second- and third-order nonlinear effects.

4. Waveguide modes

In this paper we use a recently demonstrated TFLN-on-sapphire structure as a case study [34]. The structure, its dimensions, and its dispersion characteristics are given in Fig. 2. The waveguide is designed for quasi-static nonlinear interactions involving input pulses around 2 $\mathrm {\mu }\textrm {m}$. As such, it exhibits a near-zero group velocity mismatch between the FH and SH at that wavelength (Fig. 2(d)), and has a near-zero group velocity dispersion for wavelengths from 1000 to 2100 nm (Fig. 2(e)). The dispersion characteristics are calculated with a finite difference model [37]. Fabrication processes for these waveguides are discussed in [23,45].

 

Fig. 2. Thin-film LiNbO$_3$ waveguide examined in this paper. (a) Waveguide structure and fundamental mode profiles for the input wavelength and SHG. The wafer is x-cut. The parameters are those of [34]: film thickness = 930 nm; etch depth = 600 nm; top width = 1520 nm; etch angle = 78$^{\circ }$; length = 5 mm. (b) Effective index of the fundamental mode (red) and various higher order modes (grey). (c) Normalized efficiency for SHG as a function of fundamental wavelength. Black: direct calculation using the normal modes; red: fit with an $\omega ^3$ dependence. (d) Group velocity of the fundamental mode (shifted by the group velocity at the reference wavelength). (e) Group velocity dispersion of the fundamental mode.

Download Full Size | PDF

The waveguide supports both TM and TE modes, but since the wafer is x-cut the TE modes, polarized along the crystallographic Z axis, are of primary interest for the predominant QPM interaction. However, the fundamental TE mode exhibits mode crossings with TM modes at certain frequencies. Near these avoided crossings, the normal modes are in a mixed polarization state, which is referred to as mode hybridization [46]. Thin-film LN waveguides are particularly susceptible to mode hybridization due to birefringence and because they usually lack symmetry along the vertical axis [47,48]. The latter condition applies when there are angled sidewalls, or for air top-cladding and dielectric lower cladding configurations. The normal modes are inconvenient to model because they have a different polarization state on either side of the crossing, and hence different $\chi ^{(2)}$ nonlinear mixing properties. This difficulty can be avoided by defining “pure” TE and TM modes that capture the behavior of a hypothetical waveguide with a true crossing instead of an avoided crossing. Such modes can be found by using a fitting procedure to “join up” the dispersion relation on either side of the crossing. Since the pure TE/TM modes are no longer normal modes of the waveguide, the avoided crossing manifests as a linear coupling between them [49].

In PPLN devices the $d_{33}$ nonlinear coefficient usually dominates, so, for x-cut wafers with y-propagating waveguides, nonlinear coupling involving the TE mode is far stronger. Therefore, it is reasonable to include only linear coupling between the pure-TE and pure-TM modes and assume that nonlinear coupling to the pure-TM mode is zero. Here we simply model even further by including only the pure TE fundamental mode and neglecting the other modes. We discuss this assumption briefly in section 6.4; more detailed studies of how mode hybridization influences broadband $\chi ^{(2)}$ nonlinear optics is deferred to future work.

5. Simulation examples

In this section we perform simulations based on the TFLN-on-sapphire waveguide discussed in section 4 and experimentally demonstrated in [34]. In all simulations we assume the following input-pulse parameters. Wavelength: 2100 nm. Shape: sech$^2$ input intensity. Duration: 50 fs FWHM. In addition, we assume that the normalized efficiency is given by the $\omega ^3$ fit in Fig. 2(c), and use Eqs. (14) and (30) to calculate the corresponding ${\bar {X}}{}$ coefficient for fundamental mode interactions. We choose $g_0(\omega )=((\omega /c)^{1/3} n_{\mathrm {eff}})^{-1}$ to scale the mode profiles in order to render $X_{000}$ weakly dispersive so that we can use the same nonlinear coefficient for all mixing processes involving the fundamental waveguide mode.

There are several possible ways to initialize the SSE and MSE simulations. The former requires a suitable choice of grid size, while the latter requires a suitable choice of harmonic orders. The SSE grid size can be inferred by iteratively performing simulations with finer or coarser grids until consistent solutions are obtained, or one can choose a small enough time step so that the Fourier grid’s bandwidth is sufficient to capture all sum and difference frequency generation source terms supported by optical frequencies within the full transparency window of the material. For MSE models, the harmonic orders needed to capture all non-negligible effects can be determined iteratively by adding or removing envelopes until the changes are negligible, or if an SSE simulation has been performed it can be determined from that.

In some early supercontinuum generation experiments based on reverse proton exchanged (RPE) PPLN waveguides (e.g. [16,17,50]) the essential spectral broadening mechanism involved cascaded quadratic nonlinearities (CQN) [43,51,52]. However, the low index contrast in those devices precluded dispersion engineering, so there was non-negligible group velocity mismatch between the fundamental input pulses and their second-harmonic (SH). Consequently, a large phase mismatch $\Delta k$ was needed to obtain a negative effective nonlinear refractive index compatible with few-cycle pulse compression [53]. This limited the available nonlinear coupling coefficient since the effective $n_2$ from cascading scales like $1/\Delta k$. As a result, it was necessary to consider not only the CQN process but also the influence of various other nonlinear contributions to the total effective $\chi ^{(3)}$ susceptibility. Relevant contributions included the intrinsic $\chi ^{(3)}$ nonlinearity, stimulated Raman scattering, as well as $\chi ^{(2)}$ interactions involving higher-order waveguide modes and QPM orders [33]. Another consequence of these competing nonlinear processes together with the lack of tight mode confinement was that nJ-scale pulse energies were needed to drive the SCG process compared to the pJ-scale energies required in dispersion-engineered devices.

In contrast, because of the quasi-static design of the waveguide considered here [Figs. 2(d) and (e)], it is possible to utilize the full strength of the $\chi ^{(2)}$ nonlinearity by operating, with a very small phase mismatch, in the saturated-SHG regime [29]. Therefore, we assume that the $\chi ^{(2)}$ mixing process involving only the fundamental TE mode is dominant, and focus on this interaction for the following simulations.

In subsection 5.1 we simulate supercontinuum generation as a function of pulse energy, obtaining qualitatively similar features to the recent experiments [34]. In subsection 5.2 we show how near-perfect agreement can be obtained between the MSE and SSE models. In subsection 5.3 we examine the required model complexity to capture different aspects of the spectral broadening process. In subsection 5.4 we discuss a subtlety with MSE models that can cause false excess noise amplification, and how this issue can be avoided.

5.1 Simulation versus pulse energy

We first examined the supercontinuum generation process as a function of input pulse energy and phase mismatch. Qualitatively similar behavior can be seen for various phase mismatch values that are slightly negative and of order $\Delta \beta \sim |\pi /L|$. We choose $\Delta \beta = -0.5$ mm$^{-1}$. The output spectrum as a function of input pulse energy is shown in Fig. 3. Comparison with the data shown in [34] shows reasonable agreement with the shape and pulse-energy dependence of the spectra (up to an overall factor of $2$ in the absolute energy, whose quantification is sensitive to device fabrication and experimental calibration details). For example, the spectra exhibit local minima that move in towards the input frequency around 140 THz, as apparent in the experimental results [34]. Given this evidence that our model can capture the dominant nonlinear physics involved, in the following subsections we use this case to show how the predictions of fully configured SSE and MSE are consistent, and then examine in more detail the implications of choosing various implementations of SSE vs MSE models.

 

Fig. 3. Simulated output spectra as a function of input pulse energy assuming a QPM period such that $\Delta \beta =-0.5\,mm^{-1}$. The color scale is in dB. The simulations are based on the SSE model.

Download Full Size | PDF

5.2 Single- versus multi-envelope models

Since SSE methods are generally used for interactions involving octave-spanning spectra, it is useful to compare here the results of MSE vs SSE applied to the same interaction. In this subsection we show that both SSE and MSE models can yield nearly indistinguishable results. To make a side-by-side comparison we need to extract the underlying harmonic mode envelopes from an SSE model. When the spectral content of different envelopes overlaps, this cannot be done in a single simulation. Here we show that it is possible by performing multiple simulations with different input phases $\phi _0$ since each harmonic $N$ picks up a phase shift $N\phi _0$. Consider the set of envelopes $N=\left [ 1,2,3,0\right ]$. The SSE simulation outputs are related to the underlying harmonic mode amplitude components by

The number $M$ of phases simulated should exceed the number of non-negligible harmonics that are generated so that the system of equations can be solved to determine the $\tilde {A}_{0,N}$ components. Using 8 different phases is typically sufficient, for example. Extraction of the harmonic mode amplitudes also yields the frequency comb components as each term has carrier envelope offset frequency $Nf_{\mathrm {CEO}}$. Note, this procedure works for coherent fields but may fail for quantum noise. One can distinguish whether a component has arisen from noise or not by doing a second series of simulations with semi-classical input noise switched off.

To show an application of Eq. (39) we perform SSE and MSE simulations assuming $\Delta \beta =-0.5$ mm$^{-1}$ and 10-pJ input energy. In Fig. 4(a) we show the envelopes as extracted by 8 SSE simulations followed by inversion via Eq. (39). In Fig. 4(b) we show simulations with the same parameters but with an MSE model containing the harmonic mode envelopes with $N \in \left [1,2,3,4,0\right ]$. The two simulations show excellent agreement. The dashed black line in Fig. 4(a) shows the output of the SSE simulation with input phase $\phi =0$.

 

Fig. 4. Comparison of numerical models. (a) Envelopes extracted by applying Eq. (39) to a series of 8 SSE simulations. (b) Direct MSE simulation.

Download Full Size | PDF

Most of the behavior is captured by the first four harmonics. There is a small amount of light coupled into the $N=0$ envelope by optical rectification processes. In principle there are envelopes with $N<0$ as well, but we do not include these as they are negligible. Such envelopes will usually be negligible as they require mixing between the weak spectral wings of the envelopes.

Comparison of Figs. 4(a) and 4(b) also shows a discrepancy between the envelopes’ spectral wings. In the MSE model the envelopes decay away from their nominal center frequency $N\omega _0$ as one might expect. For the SSE model this occurs as well, but the decay is not complete. When pairing an SSE model with Eq. (39), envelopes’ spectra go to zero only via cancellation between the spectra of difference simulation outputs $\tilde {A}_0(\nu ;\phi _j)$. Thus, small errors in cancellation can lead to envelopes having a finite spectral density even though there is no physical mechanism to generate them. The local peaks in the $N=3$ and $N=4$ envelopes around 120 THz in Fig. 4(a) are likely examples of this issue. Imperfect cancellation may be due to harmonic orders not included in Eq. (39) or to numerical errors in solving the matrix equation. We can infer that, due to the difficulty in getting perfect cancellation of terms, the SSE model paired with Eq. (39) is susceptible to small errors in spectral regions where envelopes are supposed to decay to zero. Nonetheless, the scale of these errors is small (e.g. the aforementioned 120-THz features lie 60 dB below the dominant spectral components).

5.3 Examination of reduced models

Having shown the consistency between SSE and MSE models in section 5.2, in this subsection we examine the required model complexity to capture the main effects. Apart from changing the set of harmonic envelopes, the other simulation parameters are identical to section 5.2. A model containing only the fundamental waveguide mode and $N\in [1,2]$ is essentially a standard SHG model (with certain generalizations and adapted to waveguide modes). The output spectrum for this model is shown in Fig. 5(a). Comparing this to the more general models of Fig. 4, one can see that there is somewhat less spectral broadening. Interestingly, after adding a 3$^{\mathrm {rd}}$ harmonic envelope the amount of broadening increases and the agreement is better [Fig. 5(b)]. Adding a 4$^{\mathrm {th}}$ harmonic on top of that [5(c)] yields even more broadening and close agreement with the simulations of 4(b), indicating that the $N=0$ envelope is not intense enough to influence the other envelopes much.

 

Fig. 5. Series of MSE simulations with different sets of envelopes. (a)-(c) show output spectra. (d)-(f) show propagation of the $N=1$ envelope through the waveguide. The three models simulated are $N\in [1,2]$ [shown in (a) and (d)], $N\in [1,2,3]$ [shown in (b) and (e)], and $N\in [1,2,3,4]$ [shown in (c) and (f)].

Download Full Size | PDF

Although the more general models yield more accurate spectral broadening, the qualitative features of how the fundamental envelope propagates are hardly changed. This can be seen by comparing Figs. 5(d)–5(f). The features are very similar in each case. This shows how a minimal SHG model is often sufficient to capture the essential mechanism involved in the spectral broadening process, even if additional envelopes (or an SSE model) are needed to capture the finer details.

In the present example, Figs. 5(d)-(f) are indicative of the quasi-static regime [29]. In the first couple of millimeters each temporal slice undergoes a quasi-cw SHG process, and one can see that lower intensity parts take longer to undergo one cycle of the forward- and back-conversion process (qualitatively similar to the behavior of the Jacobi-elliptic solutions to the basic CW three-wave mixing equations). After about 2 mm additional features start to appear, but this is not surprising since few-optical-cycle temporal features are present by that position in the waveguide.

5.4 Avoiding artificial processes in multi-envelope models

When considering noise, we take a standard approach of representing it semi-classically with a constant power spectral density corresponding to a single photon per mode, and random spectral phase. In an SSE model all noise frequencies are added to the single envelope. If multiple waveguide modes are involved, independent noise fields can be added to those as well. In contrast, if the field in each spatial mode is represented by multiple envelopes having overlapping Fourier spectra, as in an MSE model, there is not a unique way to assign the noise field to the set of envelopes. In other words, there are different ways to partition the noise so that the sum of noise fields equals the assumed total. In the limit of including an infinite series of harmonic envelopes and all interactions between them, the choice of partition should not matter because the MSE and SSE models converge in this limit. But in a realistic MSE model in which the set of harmonic envelopes is truncated, the scheme used to partition noise between the different envelopes can matter. A natural approach is to partition the semi-classical noise around each envelope’s carrier frequency so that the sum over all envelopes yields the required total noise spectrum (assigning each noise spectral component to the harmonic envelope whose carrier frequency is closest). By this centering of the noise fields around each harmonic envelope, they experience the same set of nonlinear mixing terms as the coherent pulses. Modulation instabilities leading to noise amplification within half an octave of a wave’s center frequency would be captured by this approach, making it applicable to a wide variety of situations. However, there are subtle effects that can come into play. In this subsection we examine artificial processes and how they can be avoided.

For pulses of sufficient bandwidth and intensity, intrapulse DFG can be significant, especially if the device is designed with such effects in mind. However, care is needed when adding an $N=0$ envelope because it can have subtle consequences on the semi-classical noise. We illustrate this by the four different examples shown in Fig. 6. Figure 6(a) is the same $N\in [1,2]$ model depicted in Fig. 5(a) except that quantum noise is added at all frequencies of all envelopes (in contrast to the partitioning approach described above). This distribution artificially exaggerates the total noise because the frequency grids of the different envelopes overlap. Nonetheless, no excess noise amplification occurs (see the noise floors at $\sim 10^{-6}$). The situation changes in Fig. 6(b) when we add an $N=0$ envelope: the $a_{0,0}$ envelope has broadband and strong features around 150 THz. However, this peak in the spectrum arises from an artificial optical parametric amplification (OPA) process. The process can be explained as follows. The main SHG interaction leads to a strong SH pulse around 300 THz. This pulse can drive OPA processes involving noise spectral components of the SH and $N=0$ envelopes provided those noise frequencies sum to correspond, approximately, to the SH frequency. Amplification of noise components around 150 THz has the same phase-mismatch as the main SHG process, so strong amplification can occur. The reason why this process is artificial is because $N=2$ noise components around 150 THz should also mix with the strong $N=1$ envelope via SHG, which would suppress the amplification. But this mixing process is not included because we have not included the $N=3$ envelope in the model. By adding it in Fig. 6(c) we see that the envelopes now resemble those of the general models shown in Fig. 4.

 

Fig. 6. Artificial OPA process and its suppression. (a) Shows a normal SHG interaction except that semi-classical input representing quantum noise is added in all envelopes at all frequencies. Panel (b) is similar to (a) but with the $N=0$ envelope included, to capture intrapulse DFG. The artificial OPA process is evident by the large $N=0$ spectrum around 150 THz. (c) By adding an $N=3$ envelope, the artificial OPA process is suppressed. (d) In this simulation, the input noise is partitioned between envelopes, further improving the agreement with SSE simulations. The dashed lines correspond to the assumed input noise spectrum of the correpsonding envelope.

Download Full Size | PDF

In Fig. 6(d) we show a realistic example where the noise has been partitioned around each harmonic’s center frequency (see dashed lines for the input noise seed applied to each envelope). Such partitioning is more realistic than the artificially exaggerated noise seed used in Figs. 6(a)-(c). This approach also strongly mitigates artificial OPA processes because there is no noise seed at the main frequencies susceptible to that process. When using the same noise-seeding scheme in the $N\in [0,1,2]$ model there is also no artificial OPA peak in the $N=0$ envelope. However, in that model (not shown), the $N=0$ spectrum is below the $10^{-6}$ level so it would not be visible with the vertical scale used in Fig. 6.

The preceding discussion highlights that care should be taken when adding envelopes. The $N=0$ envelope is particularly susceptible to these issues because it can permit a strong envelope (such as the SH) to amplify its own noise floor. Two MSE based solutions to the issue are (1) adding an $N=3$ envelope (so that the SH noise floor can get up-converted by the strong FH pulse) or (2) partitioning the input noise spectrum (so that envelopes do not have input noise at frequencies susceptible to such effects). Aside from these solutions, it is also beneficial to have access to an SSE model when considering coherence properties in order to perform a cross-check without the need to choose a partitioning scheme for the noise fields.

6. Discussion

The results of section 5.2 show how the same results can be obtained from SSE and MSE models provided enough envelopes are included (MSE case) or enough phases are simulated (SSE case). In this section we discuss different pros and cons of these approaches, and give an outlook.

6.1 Device design and analysis

When designing optical devices it is often useful to model a minimal set of interactions initially, for example only the coupled wave equations for SHG. A design determined by such a reduced model can then easily be tested and optimized against a more general one including all non-negligible envelopes. The analysis from Fig. 5 is an example of that approach, where successive envelopes are added to a minimal $N \in [1, 2]$ model until all of the necessary physics are included.

In general, the flexibility of MSE models to include different combinations of envelopes and to isolate the interactions between each envelope helps to better identify the key underlying mixing processes. Similarly, weak interactions can be identified and neglected before exploring a large design parameter space. For mixing processes between different envelopes that are relatively simple or involve disparate length scales, analytical or semi-analytical approximations such as cascaded quadratic nonlinearity approximations can be used to remove terms or lump them into generic self-phase modulation terms, thereby yielding a simpler model. Since all the different coupling terms appear directly, it is relatively straightforward to identify the cause of and contribution to different effects by performing simulations with different combinations of interactions included. This type of analysis may be useful both in the context of device design and in understanding experimental results. This level of flexibility is not supported by SSE models, since they treat all the different SFG and DFG processes as essentially the same. For the above reasons, Eqs. (29) and (32) can serve as an efficient tool for developing and analyzing designs.

6.2 Numerical considerations

In this subsection, we discuss some numerical aspects of the different models. Fast simulations are beneficial when exploring a device design space. The spectral bandwidth required in simulations naturally reflects the bandwidth of the different pulse components. SSE models need a very fine temporal grid to resolve all field components, while an MSE model only needs enough resolution to resolve the spectral content of the individual harmonic components. Therefore, SSE models are very inefficient for narrow-bandwidth pulses and are only relevant for octave-spanning spectra. In contrast, MSE models scale gracefully from the CW to the multiple-octave-spanning-spectrum regimes, and remain valid even for bandwidths well over an octave provided a sufficient number of envelopes are included. They are therefore suited to any $\chi ^{(2)}$ or $\chi ^{(3)}$ interaction, at many levels of approximation. However, for cases where the number of harmonics required is unknown and difficult to estimate, it can be useful to use an SSE model as a cross-check. Furthermore, for bandwidths well over an octave, the frequency grids of the different harmonic envelopes must overlap, and so the time required to perform each $z$ step may become larger than in an SSE model.

Nonetheless, there is still a potential advantage of MSE models, even for octave-spanning spectra, in terms of spatial resolution. Various nonlinear interactions can occur simultaneously, and many of these are highly phase mismatched (particularly in QPM media). With a single envelope, these interactions must be resolved in order to avoid aliasing. In contrast, highly phase mismatched terms in Eq. (32) can be solved semi-analytically using the cascaded quadratic nonlinearities approximation (discussed further in the Supplement 1), or neglected entirely where appropriate: only terms relatively close to phase-matching must be resolved numerically.

In the context of frequency comb stabilization via f-2f interferometry, multiple SSE simulations are needed to extract the underlying harmonic components and thereby estimate or optimize the strength of the radio frequency beat note that could be obtained. However, such simulations can be run in parallel before applying Eq. (39) to the output spectra. Therefore, there is not much of a disadvantage if only a small number of device configurations or pulse parameters need to be simulated. The number of different input phases needed is in the range of 4 - 8, depending on how many harmonics are generated (we used 8 here).

Similar kinds of considerations apply to devices designed for purposes other than supercontinuum generation. For example, in an optical parametric amplification configuration one could associate the different input pulses with the relevant harmonic envelopes. For non-degenerate OPA, it would be possible to derive an alternative set of equations coupling envelopes at a pump, signal, and idler frequency and all combinations of mixing frequencies corresponding to sums and differences of those frequencies. However, exploring such models is beyond the scope of this work.

The many considerations involved make it difficult and highly technical to compare simulation times in general. Nonetheless, it is useful to give some example for the case analyzed in section 5. Using a laptop computer and fine grids (0.5 fs in time and 0.5 $\mathrm {\mu }\textrm {m}$ in $z$), the SSE takes about 100 seconds while the minimal MSE model, i.e. one with just N=1 and 2 envelopes and gridded the same way as the SSE model, takes about 150 seconds in our MATLAB implementation (which we note is not fully optimized for speed). The more complex MSE models increase the simulation time further. However, if one takes advantage of the more forgiving resolution requirements of the MSE models then the simulations are much faster. For example, considering the minimal $[1,2]$ MSE model, step sizes of 3 fs (in time) and 10 $\mathrm {\mu }\textrm {m}$ (in $z$) are more than sufficient to capture the features of the fine-resolution simulations, while requiring a simulation time of under 10 seconds.

6.3 Accuracy of the nonlinear coefficients

The second- and third-order effective susceptibilities $X_{npq}$ and $X_{npqr}$ defined in Eqs. (14)–(16) can be strongly dispersive, particularly due to the changes in waveguide mode shape and size with optical frequency. The complete dispersion is included in the multiple frequency arguments to those quantities. However, in the simpler calculation given by Eqs. (22) and (32), the nonlinear polarization is found by calculating various products of the electric field envelopes in the time domain, and then multiplying the Fourier transform of these time-domain products by an appropriate coefficient that depends only on the driven (but not the driving) optical frequency.

In this approach, much of the dispersion of the effective susceptibilities is approximated by evaluating the overlap integral at specific frequencies. By a careful choice of mode normalization functions $g_n(\omega )$ the dispersion of the $X_{npq}$ functions can be made minimized, thereby reducing errors incurred by neglecting that dispersion. Different choices of $g_n(\omega )$ are discussed in in the Supplement 1, but for the TFLN-on-sapphire waveguide studied here we found that an $\sim \omega ^{-1/3}$ type scaling was sufficient to account for the leading $\omega ^3$ behavior of the normalized SHG efficiency depicted in Fig. 2.

For higher harmonics and for higher order waveguide modes, choosing suitable $g_n(\omega )$ functions and using constant nonlinear coefficients may not always capture all of the dispersion. In the MSE approach each ${\bar {X}}{}$ can, in principle, be approximated around a set of harmonic frequencies involved in a particular interaction involving a set of modes. Therefore, it is not necessary to assume that the effective susceptibilities remain constant over multiple octaves of bandwidth, and so more accurate coefficients can be used.

While some freedom in this regard is provided by the MSE model, choices of $({\underline {\bar {X}}}{}_{npq})_{P,Q}$ can still be constrained by other considerations. By choosing ${\underline {\bar {X}}}{}$ equal for all $P$ and $Q$ harmonic indices, the SSE and MSE models are the same in the limit of including all envelopes. Furthermore, the artificial OPA effects discussed in section 5.4 will be optimally suppressed when including the $N=3$ envelope only if that process has the same nonlinear coefficient, i.e. $({\underline {\bar {X}}}{}_{000})_{0,2}=({\underline {\bar {X}}}{}_{000})_{1,2}$. Since non-equal coefficients can cause subtleties like this and since the $g_0(\omega )$ function does quite a good job of removing most of the dispersion of $X_{000}$, we typically assume that all the $({\underline {\bar {X}}}{}_{000})_{P,Q}$ coefficients are the same for all harmonic indices $P$ and $Q$. However, this assumption could be relaxed if needed, and it might be particularly useful to do so for higher order waveguide modes. Another case where assuming equal coefficients would not be appropriate is THz generation, due to the large change in linear and nonlinear material properties for that spectral region. In short, the MSE model gives flexibility to account for dispersive nonlinearities, but this flexibility should be used with some caution.

6.4 Outlook

The preceding discussion helps show how both MSE and SSE models can be useful depending on the context. In a broader context of numerical modeling, much of the machinery needed for the simulations is common to both models: initializing envelopes, calculating linear and nonlinear waveguide properties, setting up Eq. (29), and so forth. Therefore, we argue that the best approach is to construct simulation software such that it supports both SSE and MSE approaches interchangeably. Then, the choice of whether to use one approach or the other (or both) for a given situation is a matter of judgement based on considerations such as those discussed in this section. While the framework we have discussed is quite general, in the future it could be extended to explore effects such as mode crossings and the resultant coupling into TM modes. For example, avoided crossings between TE and TM normal modes manifest as a linear coupling between the the pure-TE and pure-TM modes. Preliminary indications suggest that such coupling can act somewhat as a “barrier” for the spectral broadening process, so analysis of these effects may be important for future device designs. Another aspect which warrants further study is the noise and coherence properties of the spectral broadening process. For example, the sign of the group velocity dispersion coefficient is likely important in this context. The noise properties of $\chi ^{(2)}$ based supercontinuum sources may differ significantly from those based on $\chi ^{(3)}$ nonlinearities, and represents an interesting topic for further study.

By highlighting the utility of both models, and how they can be integrated and used together, our results should help in simulating and understanding the complicated dynamical processes which inevitably arise when pulses of high intensity or broad bandwidth (or both) propagate through nonlinear media in general and modern TFLN waveguides in particular. With this understanding, novel devices with lower energy requirements, improved coherence properties, or new nonlinear-optical functionalities can be considered. Such devices hold promise for advances in diverse application areas including quantum optics, precision metrology, frequency comb generation and stabilization, and absorption spectroscopy.