Nonlinear Optics

Local Modes

An arbitrary field may be described through the linear combinations of an orthogonal set, and with local modes each point along the waveguide can be represented by a distinct set. Such a generalized representation allows modal expansion to separate the physical properties of the waveguide from the underlying physics under non-longitudinally invariant conditions. Local modes are defined as the set of \(z\) dependent modes that locally satisfy Maxwell’s equations under idealized waveguide conditions. In other words, local modes are the modes of a waveguide’s instantaneous cross section. Local modes at different points along a waveguide are considered the same mode if they are related to one another through an adiabatic transformation.

The expansions used below are restricted to positive mode designators \(\left(n\right)\) and explicitly include both the forward and backward propagating modes. In the frequency-space domain \(\left(\omega, \mathbf{x}\right)\):

\[\begin{split}\begin{align} \mathbf{E} &= \sum_n \left(a_n\!\left[\omega, z\right] \, \hat{\mathbf{e}}_n\!\left[\omega, \mathbf{x}\right]\right) + \left(a_{-n}\!\left[\omega, z\right] \, \hat{\mathbf{e}}_{-n}\!\left[\omega, \mathbf{x}\right]\right) \tag{1a} \\ \mathbf{H} &= \sum_n \left(a_n\!\left[\omega, z\right] \, \hat{\mathbf{h}}_n\!\left[\omega, \mathbf{x}\right]\right) + \left(a_{-n}\!\left[\omega, z\right] \ \hat{\mathbf{h}}_{-n}\!\left[\omega, \mathbf{x}\right]\right) \tag{1b} \end{align}\end{split}\]

The first terms in the sum \(\left(+n\right)\) correspond to the forward propagating modes, and the second terms \(\left(-n\right)\) correspond to the backward propagating modes. Forward and backward propagating modes are related through the following set of equations:

\[\begin{split}\begin{align} \begin{split} {\hat{\mathbf{e}}_{-n}}_t &= +{\hat{\mathbf{e}}_{n}}_t &&& {\hat{e}_{-n}}_z &= -{\hat{e}_n}_z \\ {\hat{\mathbf{h}}_{-n}}_t &= -{\hat{\mathbf{h}}_n}_t &&& {\hat{h}_{-n}}_z &= +{\hat{h}_n}_z \\ && \beta_{-n} =& -\beta_n && \end{split} \tag{2} \end{align}\end{split}\]

At every point along the waveguide, local modes satisfy Maxwell’s equations where \(\epsilon^{\prime}\) is the idealized, lossless and longitudinally-invariant relative permittivity:

\[\begin{split}\begin{align} 0 &= \nabla_t \cdot \left(\epsilon^{\prime}_t \cdot {\hat{\mathbf{e}}_n}_t\right) - i \, \beta_n \, \epsilon^{\prime}_{zz} \ {\hat{e}_n}_z \tag{3a} \\ 0 &= \nabla_t \cdot {\hat{\mathbf{h}}_n}_t - i \, \beta_n \, {\hat{h}_n}_z \tag{3b} \\ \mathbf{0} &= \hat{\mathbf{z}} \times \left(-i \, \beta_n \, {\hat{\mathbf{e}}_n}_t - \nabla_t {\hat{e}_n}_z\right) + i \, \omega \, \mu_0 \, {\hat{\mathbf{h}}_n}_t \tag{3c} \\ \mathbf{0} &= \hat{\mathbf{z}} \times \left(-i \, \beta_n \, {\hat{\mathbf{h}}_n}_t - \nabla_t {\hat{h}_n}_z\right) - i \, \omega \, \epsilon_0 \left( \epsilon^{\prime}_t \cdot {\hat{\mathbf{e}}_n}_t \right) \tag{3d} \\ \mathbf{0} &= \nabla_t \times {\hat{\mathbf{e}}_n}_t + i \, \omega \, \mu_0 \, {\hat{h}_n}_z \hat{\mathbf{z}} \tag{3e} \\ \mathbf{0} &= \nabla_t \times {\hat{\mathbf{h}}_n}_t - i \, \omega \, \epsilon_0 \, \epsilon^{\prime}_{zz} \, {\hat{e}_n}_z \, \hat{\mathbf{z}} \tag{3f} \end{align}\end{split}\]

Multi-Mode Propagation Equations

Maxwell’s equations in the frequency-space domain \(\left(\omega, \mathbf{x}\right)\) are as follows:

\[\begin{split}\begin{align} \rho_{f} &= \nabla \cdot \left(\epsilon_0 \, \left(\mathbf{\epsilon} \cdot \mathbf{E}\right) + \mathbf{P}_\text{NL}\right) \tag{4a} \\ 0 &= \nabla \cdot \mu_0 \mathbf{H} \tag{4b} \\ \mathbf{0} &= \nabla \times \mathbf{E} + i \, \omega \, \mu_0 \ \mathbf{H} \tag{4c} \\ \mathbf{J}_{f} &= \nabla \times \mathbf{H} - i \, \omega \left(\epsilon_0 \, \left(\mathbf{\epsilon} \cdot \mathbf{E}\right) + \mathbf{P}_\text{NL}\right) \tag{4d} \end{align}\end{split}\]

Generalized unidirectional propagation equations (UPE) may be derived for the forward and backward propagating modes by substituting the modal expansion into the curl equations (see the supplemental for more details). As seen in the following results, the evolution of a pulse is driven by the following effects:

\[\begin{split}\begin{gather} \frac{\partial}{\partial z} a_n = \begin{cases} \text{linear dispersion} \\ \text{gain or loss} \\ \text{changing mode profiles} \\ \text{free sources and nonlinearity} \end{cases} \end{gather}\end{split}\]

Forward Mode Equation

The forward mode propagation equation is found by adding the two curl relations (4c + 4d).

\[\begin{split}\begin{align} \frac{\partial}{\partial z} a_n\!\left[\nu\right] &= - i \, \beta_n \, a_n \\ & \quad - \sum_{m} \left(a_m + a_{-m}\right) \iint_{A_\infty} \frac{1}{2} \left( i \, \omega \, \epsilon_0 \, {\hat{\mathbf{e}}_n^*}_t \cdot \left(\epsilon_t - \epsilon^{\prime}_t\right) \cdot {\hat{\mathbf{e}}_m}_t \right) \, dA \\ & \quad - \sum_{m} a_m \iint_{A_\infty} \frac{1}{2} \left( \frac{\partial {\hat{\mathbf{e}}_m}_t}{\partial z} \times {\hat{\mathbf{h}}^*_n}_t + {\hat{\mathbf{e}}_n^*}_t \times \frac{\partial {\hat{\mathbf{h}}_m}_t}{\partial z} \right) \cdot \hat{\mathbf{z}} \, dA \\ & \quad - \sum_{m} a_{-m} \iint_{A_\infty} \frac{1}{2} \left( \frac{\partial {\hat{\mathbf{e}}_{-m}}_t}{\partial z} \times {\hat{\mathbf{h}}^*_n}_t + {\hat{\mathbf{e}}_n^*}_t \times \frac{\partial {\hat{\mathbf{h}}_{-m}}_t}{\partial z} \right) \cdot \hat{\mathbf{z}} \, dA \\ & \quad -\iint_{A_\infty} \frac{1}{2} \left({\hat{\mathbf{e}}_n^*}_t \cdot {\mathbf{J}_{f}}_t + i \, \omega \, {\hat{\mathbf{e}}_n^*}_t \cdot {\mathbf{P}_\text{NL}}_t\right) \, dA \end{align}\end{split}\]

Backward Mode Equation

The backward mode propagation equation is found by subtracting the two curl relations (4c - 4d).

\[\begin{split}\begin{align} \frac{\partial}{\partial z} a_{-n}\!\left[\nu\right] &= - i \, \beta_{-n} \, a_{-n} \\ & \quad + \sum_m \left(a_m + a_{-m}\right) \iint_{A_\infty} \frac{1}{2} \left( i \, \omega \, \epsilon_0 \, {\hat{\mathbf{e}}_n^*}_t \cdot \left(\epsilon_t - \epsilon^{\prime}_t\right) \cdot {\hat{\mathbf{e}}_m}_t \right) \, dA \\ & \quad + \sum_m a_m \iint_{A_\infty} \frac{1}{2} \left( \frac{\partial {\hat{\mathbf{e}}_m}_t}{\partial z} \times {\hat{\mathbf{h}}^*_{-n}}_t + {\hat{\mathbf{e}}_{-n}^*}_t \times \frac{\partial {\hat{\mathbf{h}}_m}_t}{\partial z} \right) \cdot \hat{\mathbf{z}} \, dA \\ & \quad + \sum_m a_{-m} \iint_{A_\infty} \frac{1}{2} \left( \frac{\partial {\hat{\mathbf{e}}_{-m}}_t}{\partial z} \times {\hat{\mathbf{h}}^*_{-n}}_t + {\hat{\mathbf{e}}_{-n}^*}_t \times \frac{\partial {\hat{\mathbf{h}}_{-m}}_t}{\partial z} \right) \cdot \hat{\mathbf{z}} \, dA \\ & \quad +\iint_{A_\infty} \frac{1}{2} \left({\hat{\mathbf{e}}_n^*}_t \cdot {\mathbf{J}_{f}}_t + i \, \omega \, {\hat{\mathbf{e}}_n^*}_t \cdot {\mathbf{P}_\text{NL}}_t\right) \, dA \end{align}\end{split}\]

Linear Dispersion

\[\begin{gather} \frac{\partial}{\partial z} a_n\!\left[\nu\right] = \ldots - i \, \beta_n \, a_n \end{gather}\]

The first term governs linear dispersion. It leads to oscillations with position that evolve at the local value of the phase coefficient:

\[\begin{gather} a_n\!\left[z\right] \approx a_n\!\left[0\right] \, e^{-i \int^z \beta_n \, dz^\prime} \end{gather}\]

Phase Matching

Phase matching is the most important factor determining the strength of a linear or nonlinear interaction. To see this, we examine an alternative form of the modal amplitude that explicitly separates out contributions due to linear dispersion:

\[\begin{split}\begin{align} a_n &= b_n \, e^{-i \int^z \beta_n \, dz^\prime} \\ \frac{\partial a_n}{\partial z} + i \, \beta_n \, a_n &= \frac{\partial b_n}{\partial z} e^{-i \int^z \beta_n \, dz^\prime} \end{align}\end{split}\]

Replacing \(a_n\) in the propagation equations with the equivalent expressions containing \(b_n\) reveals that the phase coefficients add or subtract exponentially whenever two or more modes interact.

The implications of this can be seen in a toy model that only considers the wavenumber mismatch \(\Delta\beta\):

\[\frac{\partial}{\partial z} b = f \, e^{-i \int^z \Delta\beta \, dz^\prime}\]

In this model the placeholder term \(f\) represents the magnitude of the linear or nonlinear coupling, while \(\Delta\beta\) represents the difference between the phase coefficients of the input and output modes. In general, conversion from one (spatial or frequency) mode to another occurs until the phase difference \(\Delta\beta \, z\) accumulates to half a cycle, and the quadrature component of the complex exponential changes sign. The distance that this occurs is called the coherence length \(L_\text{coh}\), and is the maximum length scale at which non-oscillatory transfer can occur between two modes:

\[\Delta\beta \, L_\text{coh} = \pi\]

At distances shorter than the coherence length, the maximum power transfer scales with \(f^2 / \Delta\beta^2\):

\[\begin{split}\begin{align} \int^L \left(\frac{\partial}{\partial z} b \right) dz &= \int^L \left( f \, e^{-i \int^z \Delta\beta \, dz^\prime} \right) dz \\ b\!\left[L\right] &\approx -i \frac{f}{\Delta\beta} \left(1 - e^{-i \, \Delta\beta \, L}\right) \qquad \text{(holding $f$ and $\Delta\beta$ constant)} \\ \Bigl|b\!\left[L\right]\Bigr|^2 &\approx 4 \frac{f^2}{\Delta\beta^2} \sin\!\!\left[\frac{\Delta\beta \, L}{2}\right]^2 = f^2 L^2 \operatorname{sinc}\!\!\left[\frac{\Delta\beta \, L}{2}\right]^2 \end{align}\end{split}\]

When \(\Delta\beta \neq 0\), mode conversion oscillates with propagation distance. The technique of quasi-phase matching interrupts this oscillation by changing the sign of the coupling term \(f\) at intervals of the coherence length. This requires the generation of a grating with period \(\Lambda = 2\pi/\Delta\beta = 2 \, L_\text{coh}\), and is typically fabricated by inverting, or poling, the crystal axis.

When \(\Delta\beta = 0\), the process is phase matched and the mode transfer is only limited by back action affecting the dynamical quantities we have hidden in the coupling term \(f\).

Gain and Loss

\[\begin{gather} \frac{\partial}{\partial z} a_n\!\left[\nu\right] = \ldots - \sum_{m} \left(a_m + a_{-m}\right) \iint_{A_\infty} \frac{1}{2} \left( i \ \omega \ \epsilon_0 \ {\hat{\mathbf{e}}_n^*}_t \cdot \left(\epsilon_t - \epsilon^{\prime}_t\right) \cdot {\hat{\mathbf{e}}_m}_t \right) dA \end{gather}\]

Gain or loss is encoded in the difference between \(\epsilon\) and the \(\epsilon^\prime\) of the idealized waveguide. Since the lossless condition specifies that \(\epsilon^\prime\) is equal to its conjugate transpose, the difference between \(\epsilon\) and \(\epsilon^\prime\) is the anti-Hermitian part of \(\epsilon\):

\[\begin{gather} \Delta \epsilon = \epsilon - \epsilon^\prime = \frac{1}{2}\left(\epsilon - \epsilon^\dagger\right) \end{gather}\]

The relative permittivity can be decomposed in terms of the refractive indices \(n\) and the gain coefficients \(\alpha\) of the bulk materials:

\[\begin{gather} \epsilon = \mathbb{1} + \chi^{(1)} = \left(n + i \frac{c}{\omega} \frac{\alpha}{2} \right)^{2} \end{gather}\]

\[\begin{split}\begin{align} -\frac{1}{2} \left(i \, \omega \, \epsilon_0 \, {\hat{\mathbf{e}}^*_n}_t \cdot \left(\epsilon - \epsilon^\prime\right) \cdot {\hat{\mathbf{e}}_m}_t\right) &= -\frac{1}{2} \left(i \, \omega \, \epsilon_0 \, {\hat{\mathbf{e}}^*_n}_t \cdot \left(\left(n + i \frac{c}{\omega} \frac{\alpha}{2} \right)^{2} - \left(n\right)^{2}\right) \cdot {\hat{\mathbf{e}}_m}_t\right) \\ &= -\frac{1}{2} \left(i \, \omega \, \epsilon_0 \, {\hat{\mathbf{e}}^*_n}_t \cdot \left(i \, n \frac{c}{\omega} \, \alpha - \frac{1}{4} \frac{c^2}{\omega^2} \, \alpha^2\right) \cdot {\hat{\mathbf{e}}_m}_t\right) \\ &= \epsilon_0 \, c \, n_\text{eff} \ {\hat{\mathbf{e}}^*_n}_t \cdot \left(\frac{n}{n_\text{eff}}\frac{\alpha}{2} + \frac{i}{8} \frac{\alpha^2}{\beta_n}\right) \cdot {\hat{\mathbf{e}}_m}_t \end{align}\end{split}\]

where the effective refractive index \(n_\text{eff}\) is that of the output mode \(n\), with \(\beta_n = \omega \, n_\text{eff} / c\).

This form distributes contributions between an effective gain coefficient \(\alpha_\text{eff}\) and a gain-induced phase coefficient \(\beta_\alpha\).

\[\begin{gather} \frac{\partial}{\partial z} a_n\!\left[\nu\right] = \ldots + \sum_{m} \left(\frac{\alpha_\text{eff}}{2} + i \, \beta_\alpha\right) \left(a_m + a_{-m}\right) \end{gather}\]

\[\begin{split}\begin{align} {\alpha_\text{eff}}_{nm}\!\left[\nu\right] &= \iint_{A_\infty} \left(\epsilon_0 \, c \, n_\text{eff}\right) \, {\hat{\mathbf{e}}^*_n}_t \cdot \left(\frac{n}{n_\text{eff}} \alpha\right) \cdot {\hat{\mathbf{e}}_m}_t \, dA \\ {\beta_\alpha}_{nm}\!\left[\nu\right] &= \iint_{A_\infty} \left(\epsilon_0 \, c \, n_\text{eff}\right) \, {\hat{\mathbf{e}}^*_n}_t \cdot \left(\frac{1}{8} \frac{\alpha^2}{\beta_n}\right) \cdot {\hat{\mathbf{e}}_m}_t \, dA \end{align}\end{split}\]

Since \(\beta_n\) is on the order of \(10^{6}\) per meter for near-infrared light, \(\beta_\alpha\) can usually be ignored. However, it should be included when the magnitude of \(\alpha\) is unusually large such as in an absorptive material or in the gain medium of an amplifier.

From the orthogonality relationship, \({\hat{\mathbf{e}}_n^*}_t \cdot {\hat{\mathbf{e}}_m}_t\) only integrates to a Kronecker delta \(\delta_{m n}\) if the longitudinal components of the electric fields (\(\hat{e}_z\)) are negligible. Thus, besides an inhomogeneous \(\alpha\), gain and loss can also couple modes together because of a large \(\hat{e}_z\) component.

Changing Mode Profiles

\[\begin{split}\begin{align} \frac{\partial}{\partial z} a_n\!\left[\nu\right] = \ldots & - \sum_{m} a_m \iint_{A_\infty} \frac{1}{2} \left( \frac{\partial {\hat{\mathbf{e}}_m}_t}{\partial z} \times {\hat{\mathbf{h}}^*_n}_t + {\hat{\mathbf{e}}_n^*}_t \times \frac{\partial {\hat{\mathbf{h}}_m}_t}{\partial z} \right) \cdot \hat{\mathbf{z}} \, dA & \quad \text{f.p.m.} \\ & - \sum_{m} a_{-m} \iint_{A_\infty} \frac{1}{2} \left( \frac{\partial {\hat{\mathbf{e}}_{-m}}_t}{\partial z} \times {\hat{\mathbf{h}}^*_n}_t + {\hat{\mathbf{e}}_n^*}_t \times \frac{\partial {\hat{\mathbf{h}}_{-m}}_t}{\partial z} \right) \cdot \hat{\mathbf{z}} \, dA & \quad \text{b.p.m.} \end{align}\end{split}\]

These terms couple forward and backward propagating modes together due to changes in the modal profiles themselves. These can occur due to changes in the underlying geometry or material properties that determine the local modes. The terms bear a resemblance to the orthogonality relationship and can be thought of as a differential orthogonality relationship that applies when smoothly changing from one set of local modes to another.

Free Sources and Nonlinearity

\[\begin{split}\begin{align} \frac{\partial}{\partial z} a_n\!\left[\nu\right] = \ldots & -\iint_{A_\infty} \frac{1}{2} {\hat{\mathbf{e}}_n^*}_t \cdot {\mathbf{J}_{f}}_t \, dA & \quad \text{free sources} \\ & -i \, \omega \iint_{A_\infty} \frac{1}{2} \, {\hat{\mathbf{e}}_n^*}_t \cdot {\mathbf{P}_\text{NL}}_t \, dA & \quad \text{nonlinearity} \end{align}\end{split}\]

The free-source term couples modes to an external source of radiation. The extent to which oscillating current couples power into a mode depends solely on the transverse overlap between the current distribution and the mode profile. However, in optics the free current term is negligible, as coupling to an optical mode would require it to be oscillating at optical frequencies.

The nonlinear term is modified by \(i \, \omega\), which corresponds to a derivative in the time domain. In general, the nonlinear polarization is a function of the total electric field and is written as the sum of second and third order effects:

\[\begin{gather} \mathbf{P}_\text{NL} = \mathbf{P}^{\left(2\right)}_\text{NL} + \mathbf{P}^{\left(3\right)}_\text{NL} \end{gather}\]

Individually, these terms are convolutions of the time-domain electric field and the nonlinear optical susceptibilities \((\chi^{\left(2\right)}, \chi^{\left(3\right)})\). See Nonlinear Optics (Boyd 2020) for an exhaustive discussion of the important properties and symmetries associated with the nonlinear susceptibilities. In the frequency domain this corresponds to the sum of all possible interactions that yield a particular output frequency. Given the modal expansion, the nonlinear terms are represented as the sum of all possible terms given by the product of the sum of modes.

Second-Order Nonlinearity

\[\begin{split}\begin{align} \frac{\partial}{\partial z} a_n\!\left[\nu\right] = \ldots -i \, \omega \sum_{r s} \iint_{-\infty}^{\infty} & g_{n r s}^{\left(2\right)}\!\left[\nu_1, \nu_2\right] a_r\!\left[\nu_1\right] a_s\!\left[\nu_2\right] \\ & \qquad \delta\!\left[\nu - \left(\nu_1 + \nu_2\right)\right] d\nu_1 \, d\nu_2 \end{align}\end{split}\]

The second-order nonlinear polarization is defined in the time domain as the convolution of the second-order nonlinear susceptibility \(\chi^{\left(2\right)}_{i j k}\) and two delayed copies of the electric field \(\left(\mathbf{E}_j, \mathbf{E}_k\right)\). The second-order susceptibility is a third-rank tensor, and the last two indices sum over all combinations of the electric field’s Cartesian components. In the frequency domain the convolution may be interpreted as summing the contributions from all possible frequency pathways:

\[\begin{split}\begin{align} {\mathbf{P}_\text{NL}^{\left(2\right)}}_{i}\!\left[t\right] &= \epsilon_0 \iint_{-\infty}^{\infty} \chi^{\left(2\right)}_{i j k}\!\left[\tau_1, \tau_2\right] \mathbf{E}_j\!\left[t - \tau_1\right] \mathbf{E}_k\!\left[t - \tau_2\right] d\tau_1 \, d\tau_2 \\ {\mathbf{P}_\text{NL}^{\left(2\right)}}_{i}\!\left[\nu\right] &= \epsilon_0 \iint_{-\infty}^{\infty} \chi^{\left(2\right)}_{i j k}\!\left[\nu_1, \nu_2\right] \mathbf{E}_j\!\left[\nu_1\right] \mathbf{E}_k\!\left[\nu_2\right] \delta\!\left[\nu - \left(\nu_1 + \nu_2\right)\right] d\nu_1 \, d\nu_2 \end{align}\end{split}\]

Expanding the electric field in terms of the modes and substituting the equation for the frequency-domain nonlinearity into the nonlinear term from the propagation equation yields the following:

\[\begin{split}\begin{align} -i \ \omega \iint_{A_\infty} \frac{1}{2} {\hat{\mathbf{e}}^*_n}_t \cdot {\mathbf{P}_\text{NL}^{\left(2\right)}}_t \ dA = -i \ \omega \sum_{r s} \iint_{-\infty}^{\infty} & g_{n r s}^{\left(2\right)}\!\left[\nu_1, \nu_2\right] a_r\!\left[\nu_1\right] a_s\!\left[\nu_2\right] \\ & \qquad \delta\!\left[\nu - \left(\nu_1 + \nu_2\right)\right] d\nu_1 \, d\nu_2 \end{align}\end{split}\]

\[\begin{gather} g_{n r s}^{\left(2\right)}\!\left[\nu_1, \nu_2\right] = \iint_{A_\infty} \frac{1}{2} \epsilon_0 \, \chi^{\left(2\right)}_{i j k}\!\left[\nu_1, \nu_2\right] {{\hat{\mathbf{e}}_n^*}_t}_i\!\left[\nu_1 + \nu_2\right] {\hat{\mathbf{e}}_r}_j\!\left[\nu_1\right] {\hat{\mathbf{e}}_s}_k\!\left[\nu_2\right] dA \end{gather}\]

where \(g_{n r s}^{\left(2\right)}\) is the effective second-order nonlinear parameter, and the sum of all modes \(r\) and \(s\) represents all terms in the product of the sum over modes \(r\) and the sum over modes \(s\) (both forwards and backwards):

\[\begin{gather} \sum_{r s} \ldots = \left(\sum_r \ldots\right) \left(\sum_s \ldots\right) \end{gather}\]

From unit analysis and inspection of the above equations (see the supplemental for a full derivation), the effective second-order nonlinear parameter \(g^{\left(2\right)}\) can be cast in terms of an effective second-order nonlinear susceptibility \(\chi^{\left(2\right)}_\text{eff}\) that takes into account both the bulk nonlinear susceptibility and the overlap of the interacting modes:

\[\begin{gather} g_{n r s}^{\left(2\right)} = \frac{1}{2} \sqrt{\frac{{A_\text{eff}}_n}{{A_\text{eff}}_r \, {A_\text{eff}}_s}} \frac{{\chi^{\left(2\right)}_\text{eff}}_{n r s}}{\sqrt{\epsilon_0 \, c^3 \, {n_\text{eff}}_n \, {n_\text{eff}}_r \, {n_\text{eff}}_s}} \end{gather}\]

The frequency dependence of each mode in the above equation is implicitly assumed. By definition this is the same quantity as given by the first definition of \(g^{\left(2\right)}\).

Third-Order Nonlinearity

\[\begin{split}\begin{align} \frac{\partial}{\partial z} a_n\!\left[\nu\right] = \ldots -i \, \omega \sum_{q r s} \iiint_{-\infty}^{\infty} & g_{n q r s}^{\left(3\right)}\!\left[\nu_1, \nu_2, \nu_3\right] a_q\!\left[\nu_1\right] a_r\!\left[\nu_2\right] a_s\!\left[\nu_3\right] \\ & \qquad \delta\!\left[\nu - \left(\nu_1 + \nu_2 + \nu_3\right)\right] d\nu_1 \, d\nu_2 \, d\nu_3 \end{align}\end{split}\]

The third-order nonlinear polarization is defined in the time domain as the convolution of the third-order nonlinear susceptibility \(\chi^{\left(3\right)}_{i j k l}\) and three delayed copies of the electric field \(\left(\mathbf{E}_j, \mathbf{E}_k, \mathbf{E}_l\right)\). The third-order susceptibility is a fourth-rank tensor, and the last three indices sum over all combinations of the electric field’s Cartesian components. In the frequency domain the convolution may be interpreted as summing the contributions from all possible frequency pathways:

\[\begin{split}\begin{align} {\mathbf{P}_\text{NL}^{\left(3\right)}}_{i}\!\left[t\right] &= \epsilon_0 \iiint_{-\infty}^{\infty} \chi^{\left(3\right)}_{i j k l}\!\left[\tau_1, \tau_2, \tau_3\right] \mathbf{E}_j\!\left[t - \tau_1\right] \mathbf{E}_k\!\left[t - \tau_2\right] \mathbf{E}_l\!\left[t - \tau_3\right] d\tau_1 \ d\tau_2 \ d\tau_3 \\ {\mathbf{P}_\text{NL}^{\left(3\right)}}_{i}\!\left[\nu\right] &= \epsilon_0 \iiint_{-\infty}^{\infty} \chi^{\left(3\right)}_{i j k l}\!\left[\nu_1, \nu_2, \nu_3\right] \mathbf{E}_j\!\left[\nu_1\right] \mathbf{E}_k\!\left[\nu_2\right] \mathbf{E}_l\!\left[\nu_3\right] \delta\!\left[\nu - \left(\nu_1 + \nu_2 + \nu_3\right)\right] d\nu_1 \, d\nu_2 \, d\nu_3 \end{align}\end{split}\]

Expanding the electric field in terms of the modes and substituting the equation for the frequency-domain nonlinearity into the nonlinear term from the propagation equation yields the following:

\[\begin{split}\begin{align} -i \, \omega \iint_{A_\infty} \frac{1}{2} {\hat{\mathbf{e}}^*_n}_t \cdot {\mathbf{P}_\text{NL}^{\left(3\right)}}_t \, dA = -i \, \omega \sum_{q r s} \iiint_{-\infty}^{\infty} & g_{n q r s}^{\left(3\right)}\!\left[\nu_1, \nu_2, \nu_3\right] a_q\!\left[\nu_1\right] a_r\!\left[\nu_2\right] a_s\!\left[\nu_3\right] \\ & \qquad \delta\!\left[\nu - \left(\nu_1 + \nu_2 + \nu_3\right)\right] d\nu_1 \, d\nu_2 \, d\nu_3 \end{align}\end{split}\]

\[\begin{gather} g_{n q r s}^{\left(3\right)}\!\left[\nu_1, \nu_2, \nu_3\right] = \iint_{A_\infty} \frac{1}{2} \epsilon_0 \, \chi^{\left(3\right)}_{i j k l}\!\left[\nu_1, \nu_2, \nu_3\right] {{\hat{\mathbf{e}}_n^*}_t}_i\!\left[\nu_1 + \nu_2 + \nu_3\right] {\hat{\mathbf{e}}_q}_j\!\left[\nu_1\right] {\hat{\mathbf{e}}_r}_k\!\left[\nu_2\right] {\hat{\mathbf{e}}_s}_l\!\left[\nu_3\right] dA \end{gather}\]

where \(g_{n q r s}^{\left(3\right)}\) is the effective third-order nonlinear parameter, and the sum of all modes \(q\), \(r\), and \(s\) represents all terms in the product of the sum over modes \(q\), \(r\), and \(s\) (both forwards and backwards):

\[\begin{gather} \sum_{q r s} \ldots = \left(\sum_q \ldots\right) \left(\sum_r \ldots\right) \left(\sum_s \ldots\right) \end{gather}\]

Similar to the second-order nonlinearity (see the supplemental for more details), the effective third-order nonlinear parameter \(g^{\left(3\right)}\) can also be cast in terms of an effective nonlinear susceptibility \(\chi^{\left(3\right)}_\text{eff}\) that takes into account both the bulk nonlinear susceptibility and the overlap of the interacting modes:

\[\begin{gather} g_{n q r s}^{\left(3\right)} = \frac{1}{2} \sqrt{\frac{{A_\text{eff}}_n}{{A_\text{eff}}_q \, {A_\text{eff}}_r \, {A_\text{eff}}_s}} \frac{{\chi^{\left(3\right)}_\text{eff}}_{n q r s}}{\epsilon_0 \, c^2 \, \sqrt{{n_\text{eff}}_n \, {n_\text{eff}}_q \, {n_\text{eff}}_r \, {n_\text{eff}}_s}} \end{gather}\]

The frequency dependence of each mode in the above equation is implicitly assumed. By definition this is the same quantity as given by the first definition of \(g^{\left(3\right)}\).

Single-Mode UPE

Considering only a single forward propagating mode simplifies the propagation equation. Due to phase mismatch between different modes this is typically a good approximation:

\[\begin{split}\begin{align} \frac{\partial}{\partial z} a\!\left[\nu\right] &= - i \, \beta \, a & \quad \text{linear dispersion} \\ & \quad - a \iint_{A_\infty} \frac{1}{2} \left( i \, \omega \, \epsilon_0 \, \hat{\mathbf{e}}^*_t \cdot \left(\epsilon_t - \epsilon^{\prime}_t \right) \cdot {\hat{\mathbf{e}}}_t\right) dA & \quad \text{gain or loss} \\ & \quad - \iint_{A_\infty} \frac{1}{2} \left( i \, \omega \, \hat{\mathbf{e}}^*_t \cdot {\mathbf{P}_\text{NL}}_t \right) \ dA & \quad \text{nonlinearity} \end{align}\end{split}\]

Gain and Loss

\[\begin{split}\begin{align} \frac{\partial}{\partial z} a\!\left[\nu\right] &= \ldots - a \iint_{A_\infty} \frac{1}{2} \left( i \, \omega \, \epsilon_0 \, \hat{\mathbf{e}}^*_t \cdot \left(\epsilon_t - \epsilon^{\prime}_t \right) \cdot {\hat{\mathbf{e}}}_t\right) dA \\ &= \ldots + \left(\frac{\alpha_\text{eff}}{2} + i \, \beta_\alpha \right) a \end{align}\end{split}\]

The overlap between the change in relative permittivity and the spatial mode determines the effective propagation gain or loss. Each term is a sum weighted with the orthonormal spatial mode of the waveguide:

\[\begin{split}\begin{align} \alpha_\text{eff}\!\left[\nu\right] &= \iint_{A_\infty} \left( \epsilon_0 \, c \, n_\text{eff}\right) \, \hat{\mathbf{e}}^*_t \cdot \left(\frac{n}{n_\text{eff}} \alpha\right) \cdot {\hat{\mathbf{e}}}_t \ dA \\ \beta_{\alpha}\!\left[\nu\right] &= \iint_{A_\infty} \left( \epsilon_0 \, c \, n_\text{eff} \right) \, \hat{\mathbf{e}}^*_t \cdot \left(\frac{1}{8} \frac{\alpha^2}{\beta_\text{eff}}\right) \cdot {\hat{\mathbf{e}}}_t \ dA \end{align}\end{split}\]

The above relationships take a simpler form in the limit that the longitudinal component of the electric field \(\left(\hat{e}_z\right)\) can be ignored and that \(\alpha\) is isotropic and homogenous across the transverse plane:

\[\begin{split}\begin{align} \alpha_\text{eff} &\approx \alpha \\ \beta_\alpha &\approx \frac{1}{8} \frac{\alpha^2}{\beta_\text{eff}} \end{align}\end{split}\]

If all other effects besides linear dispersion are negligible, then the amplitude has the solution below. Positive \(\alpha\) represents gain and negative \(\alpha\) represents loss:

\[\begin{gather} a\!\left[z\right] \approx a\!\left[0\right] e^{\int^{z} \left(-i \ \left(\beta - \beta_\alpha\right) + \frac{1}{2}\alpha\right) \ dz^\prime} \end{gather}\]

Second Order Nonlinearity

\[\begin{split}\begin{align} \frac{\partial}{\partial z} a\!\left[\nu\right] = \ldots - i \, \omega \iint_{-\infty}^{\infty} & g^{\left(2\right)}\!\left[\nu_1, \nu_2\right] a\!\left[\nu_1\right] a\!\left[\nu_2\right] \\ & \qquad \delta\!\left[\nu - \left(\nu_1 + \nu_2\right)\right] d\nu_1 \, d\nu_2 \end{align}\end{split}\]

where \(g^{\left(2\right)}\) is the effective second-order nonlinear parameter defined as:

\[\begin{gather} g^{\left(2\right)}\!\left[\nu_1, \nu_2\right] = \frac{1}{2} \epsilon_0 \iint_{A_\infty} \chi^{\left(2\right)}_{i j k}\!\left[\nu_1, \nu_2\right] {\hat{\mathbf{e}}^*_t}_i\!\left[\nu_1 + \nu_2\right] \hat{\mathbf{e}}_j\!\left[\nu_1\right] \hat{\mathbf{e}}_k\!\left[\nu_2\right] dA \end{gather}\]

In general, the nonlinear parameter may be refactored in terms of the characteristic waveguide properties and an effective nonlinear susceptibility \(\chi^{\left(2\right)}_\text{eff}\) that considers both the bulk nonlinear susceptibility and the overlap of the interacting modes:

\[\begin{gather} g^{\left(2\right)}\!\left[\nu_1, \nu_2\right] = \frac{1}{2} \sqrt{\frac{{A_\text{eff}}\!\left[\nu\right]}{{A_\text{eff}}\!\left[\nu_1\right] \, {A_\text{eff}}\!\left[\nu_2\right]}} \frac{{\chi^{\left(2\right)}_\text{eff}}\!\left[\nu_1, \nu_2\right]}{\sqrt{\epsilon_0 \, c^3 \, {n_\text{eff}}\!\left[\nu\right] \, {n_\text{eff}}\!\left[\nu_1\right] \, {n_\text{eff}}\!\left[\nu_2\right]}} \end{gather}\]

where \(\nu\) is the sum of \(\nu_1\) and \(\nu_2\).

Cascaded Effects

The second-order nonlinearity is best known for sum \((\omega = \omega_1 + \omega_2)\) and difference \((\omega = \omega_1 - \omega_2)\) frequency generation, but in certain cases cascaded effects can also be important. Cascaded nonlinearity produces a higher-order nonlinearity through repeated action of a lower-order nonlinearity. One of the more interesting phenomena is the effective nonlinear refractive index \(n_2^\text{eff}\) due to phase-mismatched second-harmonic generation. This effect has the same properties as the equivalent third-order nonlinearity discussed below. For sufficiently large values of the phase mismatch \(\Delta\beta\):

\[n_2^\text{eff} \approx - \frac{\nu}{\epsilon_0 \, c^2} \frac{\chi^{(2)}_\text{eff}}{n_{2\omega} \, n_\omega^2} \frac{\pi}{\Delta\beta}\]

Depending on the sign of the phase mismatch, the nonlinear index \(n_2^\text{eff}\) can be either positive or negative. In the small phase mismatch regime, this equation breaks down and the nonlinear index is governed by a Kramers-Kronig relationship with the pump depletion rate (see DeSalvo (1992) for more details).

Third-Order Nonlinearity

\[\begin{split}\begin{align} \frac{\partial}{\partial z} a\!\left[\nu\right] = \ldots - i \, \omega \iiint_{-\infty}^{\infty} & g^{\left(3\right)}\!\left[\nu_1, \nu_2, \nu_3\right] a\!\left[\nu_1\right] a\!\left[\nu_2\right] a\!\left[\nu_3\right] \\ & \qquad \delta\!\left[\nu - \left(\nu_1 + \nu_2 + \nu_3\right)\right] d\nu_1 \, d\nu_2 \, d\nu_3 \end{align}\end{split}\]

where \(g^{\left(3\right)}\) is the effective third-order nonlinear parameter defined as:

\[\begin{gather} g^{\left(3\right)}\!\left[\nu_1, \nu_2, \nu_3\right] = \iint_{A_\infty} \frac{1}{2} \epsilon_0 \, \chi^{\left(3\right)}_{i j k l}\!\left[\nu_1, \nu_2, \nu_3\right] {\hat{\mathbf{e}}^*_t}_i\!\left[\nu_1 + \nu_2 + \nu_3\right] \hat{\mathbf{e}}_j\!\left[\nu_1\right] \hat{\mathbf{e}}_k\!\left[\nu_2\right] \hat{\mathbf{e}}_l\!\left[\nu_3\right] dA \end{gather}\]

In general, the nonlinear parameter may be refactored in terms of the characteristic waveguide properties and an effective nonlinear susceptibility \(\chi^{\left(3\right)}_\text{eff}\) that considers both the bulk nonlinear susceptibility and the overlap of the interacting modes:

\[\begin{gather} g^{\left(3\right)}\!\left[\nu_1, \nu_2, \nu_3\right] = \frac{1}{2} \sqrt{\frac{A_\text{eff}\!\left[\nu\right]}{A_\text{eff}\!\left[\nu_1\right] \, A_\text{eff}\!\left[\nu_2\right] \, A_\text{eff}\!\left[\nu_3\right]}} \frac{\chi^{\left(3\right)}_\text{eff}\!\left[\nu_1, \nu_2, \nu_3\right]}{\epsilon_0 \, c^2 \sqrt{n_\text{eff}\!\left[\nu\right] \, n_\text{eff}\!\left[\nu_1\right] \, n_\text{eff}\!\left[\nu_2\right] \, n_\text{eff}\!\left[\nu_3\right]}} \end{gather}\]

where \(\nu\) is the sum of \(\nu_1\), \(\nu_2\), and \(\nu_3\).

Self-Phase Modulation

Of the \(\chi^{\left(3\right)}\) processes, self-phase modulation is always phase matched and corresponds to the nonlinear interaction between positive and negative frequency components from the same three modes \(\left(\omega = \omega + \omega - \omega\right)\). In this case the effective nonlinearity collapses to a function of a single variable:

\[\begin{gather} g_\text{spm}^{\left(3\right)}\!\left[\nu\right] = \frac{1}{2} \frac{\chi^{\left(3\right)}_\text{eff}\!\left[\nu, \nu, -\nu\right]}{\epsilon_0 \, c^2 \, n_\text{eff}^2 \, A_\text{eff}} \end{gather}\]

The frequency dependence of the effective refractive index and area are implicit. If one assumes that phase mismatch suppresses all other interactions, then \(g^{\left(3\right)}_\text{spm}\) can be completely removed from the frequency integrals and the convolution can be cast as a multiplication in the time domain:

\[\begin{gather} \frac{\partial}{\partial z} a\!\left[\nu\right] \approx \ldots - i \, \omega \, g^{\left(3\right)}_\text{spm} \, \mathscr{F}\!\left[a\left[t\right]^3\right] \end{gather}\]

With a little manipulation, one arrives at the familiar nonlinear parameter \(\gamma\):

\[\begin{gather} \gamma\!\left[\nu\right] = \frac{3}{2} \, \omega \, g^{\left(3\right)}_\text{spm} \end{gather}\]

where the factor of \(3\) is due to the degeneracy of the self-phase modulation term \(\left(\omega = \omega_1 + \omega_2 - \omega_3\right)\) in the triple product of the carrier-resolved amplitude \((a\left[t\right]^3)\). The factor of \(1/2\) is due to the \(\gamma\) parameter being defined against the analytic, root-mean-square amplitude \(a'\), where \(a'\!\left[\nu\right] = \sqrt{2} \, a\!\left[\nu\right]\):

\[\begin{gather} \frac{\partial}{\partial z} a'\!\left[\nu\right] \approx \ldots -i \, \gamma \, \mathscr{F}\!\left[a' \, \left|a'\right|^2\right] \end{gather}\]

The expected nonlinear phase shift is proportional to the instantaneous power and propagation distance:

\[\begin{split}\begin{gather} \phi_{NL} = \gamma \, P \, L \\ P = \left|a'\right|^2 = \text{RMS power} \qquad L = \text{propagation length} \end{gather}\end{split}\]

Nonlinear Refractive Index

The \(\gamma\) factor is often cast in terms of the nonlinear Kerr parameter \(n_2\), which describes the change in refractive index due to the instantaneous, time-averaged optical power. For single frequency, linearly polarized light:

\[\begin{gather} n = n_0 + n_2 \left<I\right> \end{gather}\]

where \(I\) is the intensity of the light defined through the Poynting vector, and \(\left<\ldots\right>\) denotes the instantaneous time average as given by the absolute mean square of the electric field’s complex envelope.

The third order nonlinear susceptibility \(\chi^{\left(3\right)}\) may be cast in terms of the Kerr nonlinear index \(n_2\) by relating the two through the definition of the relative permittivity (see the supplemental for more details):

\[\begin{split}\begin{align} \chi^{\left(3\right)} &= \frac{4}{3} \, \epsilon_0 \, c \, n_0^2 \, n_2 \\ \gamma &= \frac{\omega \, n_2}{c \, A_\text{eff}} \end{align}\end{split}\]

Optical Shock

The shock term \(\gamma_1\) typically seen in the literature is an artifact of approximating the time derivative applied to the nonlinear polarization. In practice this means using a constant nonlinear parameter, i.e. \(\gamma_0 = \gamma\!\left[\omega_0\right]\). If everything is properly accounted for the shock term \(\gamma_1\) and the standard nonlinear term \(\gamma_0\) collapse into the frequency-dependent nonlinear parameter \(\gamma\!\left[\omega\right]\) consistent with the previous sections:

\[\begin{gather} \gamma\!\left[\omega\right] = \frac{\omega \, n_2}{c \, A_\text{eff}} = \frac{3}{2} \, \omega \, g^{\left(3\right)}_\text{spm} \end{gather}\]

The following arguments are made using the Fourier transform of the complex envelope:

\[\begin{split}\begin{align} a\!\left[t\right] &= \int_{-\infty}^{\infty} a\!\left[\omega - \omega_0\right] e^{+i \left(\left(\omega - \omega_0\right) t\right)} \frac{d\omega}{2 \, \pi} \\ a\!\left[\omega - \omega_0\right] &= \int_{-\infty}^{\infty} a\!\left[t\right] e^{-i \left(\left(\omega - \omega_0\right) t\right)} \, dt \end{align}\end{split}\]

Starting in the time domain, we expand \(\gamma[t]\) in terms of \(\gamma_0\) and \(\gamma_1\). Note, the sign before \(\gamma_1\) is due to our definition of the Fourier transform:

\[\gamma\!\left[t\right] = \gamma_0 - i \, \gamma_1 \frac{\partial}{\partial t}\]

With \(\gamma_1 = \gamma_0/\omega_0\), the Fourier transform gives the frequency-dependent \(\gamma[\omega]\):

\[\begin{split}\begin{align} \left(\gamma_0 - i \, \gamma_1 \, \frac{\partial}{\partial t}\right) f\!\left[t\right] &= \gamma_0 \left(1 - \frac{i}{\omega_0} \frac{\partial}{\partial t}\right) f\!\left[t\right] \\ \text{Fourier transform $t \to \omega - \omega_0$} & \Rightarrow \gamma_0 \left(1 + \frac{1}{\omega_0} \left(\omega - \omega_0\right)\right) f\!\left[\omega - \omega_0\right] \\ &= \gamma_0 \frac{\omega}{\omega_0} \, f\!\left[\omega - \omega_0\right] \\ &= \gamma\!\left[\omega\right] \, f\!\left[\omega - \omega_0\right] \end{align}\end{split}\]

Raman Nonlinearity

Raman nonlinearity is a third-order effect that couples energy between light and vibrational modes of the substrate. It is represented in terms of a unitless nonlinear response parameter \(R\) that takes two of the three frequencies as input. The symmetrized form is given below:

\[\begin{gather} \chi^{(3)}\!\left[\nu_1, \nu_2, \nu_3\right] \approx \chi^{(3)}_\text{eff} \left( \frac{1}{3} \left( R\!\left[\nu_1 + \nu_2\right] + R\!\left[\nu_2 + \nu_3\right] + R\!\left[\nu_3 + \nu_1\right]\right)\right) \end{gather}\]

The nonlinear response parameter is itself split between the instantaneous nonlinearity of the Kerr effect and the delayed Raman response function \(h_R\):

\[\begin{split}\begin{gather} R\!\left[t\right] = f_R \, h_R\!\left[t\right] + \left(1 - f_R\right) \, \delta\!\left[t\right] \\ 1 = \int_{0}^{\infty} R\!\left[t\right] \, dt = \int_{0}^{\infty} h_R\!\left[t\right] \, dt \end{gather}\end{split}\]

The Raman fraction \(f_R\) sets the relative strength of the instantaneous and delayed effects. The Raman response \(h_R\) primarily acts on the difference frequencies of the two inputs. In silica it is strongest at offsets between 10 to 15 THz.

In the time domain, the nonlinear response parameter acts through convolution with the square of the time-domain amplitude:

\[\begin{gather} \frac{\partial}{\partial z} a\!\left[\nu\right] = \ldots - i \, \omega \, g^{(3)}_\text{eff}\!\left[\nu\right] \, \mathscr{F}\!\left[a\!\left[t\right] \int^{\infty}_{0} R\!\left[\tau\right] a\left[t - \tau\right]^2 d\tau\right] \end{gather}\]

Numerical Implementation

The goal of a fast and accurate numerical implementation is to cast the equations into a form that minimizes both the numerical complexity and the degree of approximation. Besides just reducing the numerical complexity of individual terms, it is also important to use an accurate and efficient numerical differential equation solver. For this purpose, PyNLO uses the embedded 4th-order Runge-Kutta method of Balac and Mahé (2013).

Nonlinear Terms

The numerical complexity of the convolutions required to calculate the second- and third-order nonlinear terms go as \(N^2\) and \(N^3\) respectively. However, if the nonlinear terms can be put in separable form, i.e. as the product of single variable functions, then the convolution can be recast as a product of Fourier transforms. Using the fast Fourier transform (fft) the complexity of the resulting convolutions scale as \(N \log\left[N\right]\). For the number of points typically used in nonlinear simulations, that difference can mean orders of magnitude improvement to the run time.

Vector Decomposition

The most general procedure involves numerically calculating all terms in the 2nd- and 3rd-order \(g\) parameters and then decomposing the tensors as the sum of the product of separable functions. In general, the \(g\) parameters may be written as the product of a function \(h\), dependent on the output frequency \(\nu\), and the sum of the products of functions \(\eta\), dependent on the input frequencies \(\nu_i\):

\[\begin{split}\begin{align} g^{\left(2\right)}\!\left[\nu_1, \nu_2\right] &= h^{\left(2\right)}\!\left[\nu\right] \sum_n \eta_n^{\left(2\right)}\!\left[\nu_1\right] \eta_n^{\left(2\right)}\!\left[\nu_2\right] \\ g^{\left(3\right)}\!\left[\nu_1, \nu_2, \nu_3\right] &= h^{\left(3\right)}\!\left[\nu\right] \sum_n \eta_n^{\left(3\right)}\!\left[\nu_1\right] \eta_n^{\left(3\right)}\!\left[\nu_2\right] \eta_n^{\left(3\right)}\!\left[\nu_3\right] \end{align}\end{split}\]

Before taking the Fourier transform, each function in the sum \((\eta_n)\) is used to weight the spectral amplitudes \(a\):

\[\begin{gather} b_n\!\left[\nu\right] = a \, \eta_n \end{gather}\]

\[\begin{split}\begin{align} \frac{\partial}{\partial z} a\!\left[\nu\right] = ... & - i \, \omega \, h^{\left(2\right)}\!\left[\nu\right] \, \mathscr{F}\!\left[\sum_n b_n\left[t\right]^2 \right] \\ & - i \, \omega \, h^{\left(3\right)}\!\left[\nu\right] \, \mathscr{F}\!\left[\sum_n b_n\left[t\right]^3 \right] \end{align}\end{split}\]

Unit Decomposition

If effects due to mode overlap can be ignored and the nonlinear susceptibility only depends on the output frequency, separability occurs naturally. In this case, the summation collapses to a single \(h\) and \(\eta\) term dependent on the effective refractive index, area, and nonlinear susceptibilities. The decomposition below is based on unit analysis of the effective nonlinear parameters (see the supplemental for more details) and is equivalent to calculating the nonlinear terms as powers of the electric field:

\[\begin{split}\begin{gather} h\!\left[\nu\right] \approx \frac{1}{2} \sqrt{\frac{\epsilon_0 \, A_\text{eff}}{c \, n_\text{eff}}} \, \chi_\text{eff} \\ \eta\!\left[\nu\right] \approx \frac{1}{\sqrt{\epsilon_0 \, c \, n_\text{eff} \, A_\text{eff}}} \end{gather}\end{split}\]

Phase-Matched Pathways

The most straightforward approximation takes advantage of the selective nature of phase matching. Phase-matched frequency combinations dominate over all other nonlinear interactions. For a particular output frequency, the \(g\) parameters that correspond to these interactions can be pulled out of the integral and the unweighted spectral amplitudes can be Fourier transformed. The unweighted terms still contain a multitude of nonlinear pathways, but as phase mismatch suppresses all other contributions only the phase-matched pathway needs to have an accurate coupling term:

\[\begin{split}\begin{align} g^{\left(2\right)}\!\left[\nu_1, \nu_2\right] &\approx g^{\left(2\right)}_\text{eff}\!\left[\nu\right] \\ g^{\left(3\right)}\!\left[\nu_1, \nu_2, \nu_3\right] &\approx g^{\left(3\right)}_\text{eff}\!\left[\nu\right] \end{align}\end{split}\]

\[\begin{split}\begin{align} \frac{\partial}{\partial z} a\!\left[\nu\right] = ... & - i \, \omega \, g^{\left(2\right)}_\text{eff}\!\left[\nu\right] \, \mathscr{F}\!\left[ a\left[t\right]^2 \right] \\ & - i \, \omega \, g^{\left(3\right)}_\text{eff}\!\left[\nu\right] \, \mathscr{F}\!\left[a\left[t\right]^3 \right] \end{align}\end{split}\]

The Analytic Representation

The analytic representation is an important concept for the theoretical and numerical manipulation of optical pulses. Since the physical \(\mathbf{E}\) and \(\mathbf{H}\) fields are real, as mentioned in Fourier Transforms, it is possible to represent them in analytic form using amplitudes from only the positive-frequency side of the spectrum. This is straightforward when only linear interactions are considered but becomes more difficult when nonlinear terms are involved. Nonlinearity mixes positive and negative frequency components together in an inseparable way.

For the second-order nonlinearity the following interactions are possible:

\(\chi^{\left(2\right)}\) Process	Center	Range
\(+\nu_1\), \(+\nu_2\)	\(+2 \, \nu_c\)	\(0 \le \nu \le 2 \, \nu_\text{max}\)
\(-\nu_1\), \(+\nu_2\)	\(0\)	\(-\nu_\text{max} \le \nu \le \nu_\text{max}\)
\(+\nu_1\), \(-\nu_2\)	\(0\)	\(-\nu_\text{max} \le \nu \le \nu_\text{max}\)
\(-\nu_1\), \(-\nu_2\)	\(-2 \, \nu_c\)	\(-2 \, \nu_\text{max} \le \nu \le 0\)

Similarly, for the third-order nonlinearity:

\(\chi^{\left(3\right)}\) Process	Center	Range
\(+\nu_1\), \(+\nu_2\), \(+\nu_3\)	\(+3 \, \nu_c\)	\(0 \le \nu \le 3 \, \nu_\text{max}\)
\(-\nu_1\), \(+\nu_2\), \(+\nu_3\)	\(+\nu_c\)	\(-\nu_\text{max} \le \nu \le 2 \, \nu_\text{max}\)
\(+\nu_1\), \(-\nu_2\), \(+\nu_3\)	\(+\nu_c\)	\(-\nu_\text{max} \le \nu \le 2 \, \nu_\text{max}\)
\(-\nu_1\), \(-\nu_2\), \(+\nu_3\)	\(-\nu_c\)	\(-2 \, \nu_\text{max} \le \nu \le \nu_\text{max}\)
\(+\nu_1\), \(+\nu_2\), \(-\nu_3\)	\(+\nu_c\)	\(-\nu_\text{max} \le \nu \le 2 \, \nu_\text{max}\)
\(-\nu_1\), \(+\nu_2\), \(-\nu_3\)	\(-\nu_c\)	\(-2 \, \nu_\text{max} \le \nu \le \nu_\text{max}\)
\(+\nu_1\), \(-\nu_2\), \(-\nu_3\)	\(-\nu_c\)	\(-2 \, \nu_\text{max} \le \nu \le \nu_\text{max}\)
\(-\nu_1\), \(-\nu_2\), \(-\nu_3\)	\(-3 \, \nu_c\)	\(-3 \, \nu_\text{max} \le \nu \le 0\)

where \(\nu_c\) is the central or average frequency of the pulse.

Analytic Nonlinear Terms

Negative frequency components at the input to a nonlinear process can always be replaced by conjugated positive frequency components, the difficulty arises when restricting the range of output frequencies to positive values. To reduce the numerical complexity, nonlinear terms are best calculated as products in the time domain. However, the time-domain products are not strictly analytic as in general they each span both positive and negative frequencies. If one were to perform these products in the time domain with analytic grids that only support the positive half of the spectrum, the negative frequency components would alias into positive frequencies. This could lead to erroneous results, especially if an aliased term is inadvertently phase matched.

The simplest way around this is to cheat the analyticity and perform the frequency-to-time-domain Fourier transform with enough extra points to allow the negative frequency components to alias outside of the frequency range of interest. The “analytic” nonlinear terms are thus all entries in the above tables with the potential to generate positive output frequencies. After Fourier transforming back to the frequency domain, the extra points are dropped and only the analytic results remain. If the analytic grid has a support from \(0\) to \(\nu_\text{max}\), then the second-order nonlinearity would need a grid that supports at least another \(\nu_\text{max}\) worth of points, while the third-order nonlinearity would need an extra \(2 \, \nu_\text{max}\), requiring a grid that supports a total of \(2 \, \nu_\text{max}\) and \(3 \, \nu_\text{max}\) points respectively.

Analytic Nonlinear Propagation Equations

In the limit of constant mode area and instantaneous nonlinear susceptibility, the analytic second- and third-order nonlinear propagation equations have the following form:

\[\begin{split}\begin{gather} \frac{\partial}{\partial z} a\!\left[\nu\right] = \ldots - i \, \omega \, g^{\left(2\right)}_\text{eff} \, \mathscr{F}\!\left[ 2 \, \left|a\right|^2 \, e^{-i \, \omega_c \, t} + a^2 \, e^{i \, \omega_c \, t} \right] \\ g^{\left(2\right)}_\text{eff} = \frac{1}{2} \frac{\chi^{\left(2\right)}_\text{eff}}{\sqrt{\epsilon_0 \, c^3 \, n_\text{eff}^3 \, A_\text{eff}}} \end{gather}\end{split}\]

\[\begin{split}\begin{gather} \frac{\partial}{\partial z} a\!\left[\nu\right] = \ldots - i \, \omega \, g^{\left(3\right)}_\text{eff} \, \mathscr{F}\!\left[ 3 \, \left|a\right|^2 a + 3 \, \left|a\right|^2 a^* \, e^{-i \, 2 \, \omega_c \, t} + a^3 \, e^{i \, 2 \, \omega_c \, t} \right] \\ g^{\left(3\right)}_\text{eff} = \frac{1}{2} \frac{\chi^{\left(3\right)}_\text{eff}}{\epsilon_0 \, c^2 \, n_\text{eff}^2 \, A_\text{eff}} \end{gather}\end{split}\]

where the nonlinear terms have been constructed using all possible combinations with positive output frequencies. Conjugate amplitudes have been substituted for amplitudes at negative-frequency, i.e. \(a[-\nu] \to a[\nu]^*\).

Each nonlinear process has been frequency shifted by the respective difference between the input and output center frequency \(\nu_c\). This is required because, mathematically speaking, the analytic spectra are centered at baseband. Since time-domain products of waveforms with a frequency range centered about the origin result in waveforms with frequency ranges centered about the origin, frequency shifts are necessary in order to maintain the correct center frequency of the nonlinear process.

PyNLO does not use the analytic representation when calculating nonlinear interactions and it does so for several reasons. The result of extending the analytic representation to nonlinear terms is a set of equations that are more complicated than their equivalent real-valued implementation. Given that the multiplication of complex numbers takes 3 operations to complete, the analytic representation requires 11 to 19 array operations compared to the 1 to 2 of the real-valued representation. For the same number of points, the rffts used in the real-valued representation are also more performant than the ffts used in the analytic representation. Additionally, if the frequency grids are aligned to the origin (as mentioned in the FFT Cookbook), one representation can be easily transformed to the other using array indexing operations. Ultimately, the analytic representation, while useful when modeling linear effects, is more of a burden than a benefit for nonlinear effects.

Analytic Normalization

When comparing against other formulas in the literature it is important to consider the normalization of the analytic amplitudes. There are two common alternatives to formulas given above, which use the amplitudes as given by the positive half of the double-sided spectrum:

The first rescales the analytic representation such that the integral of the absolute-squared amplitude over positive frequencies is equal to the total energy \(E\):

\[\begin{split}\begin{gather} a'\!\left[\nu\right] = \sqrt{2} \, a \\ E = \int_0^{\infty} \left|a'\right|^2 d\nu = \int_{-\infty}^{\infty} \left|a\right|^2 d\nu \end{gather}\end{split}\]

The second is derived using the time-domain Hilbert transform to zero out the negative frequency components:

\[\begin{split}\begin{gather} a''\!\left[\nu\right] = 2 \, a \\ a''\!\left[t\right] = a + i \, \mathcal{H}\!\left[a\right] \end{gather}\end{split}\]