Nonlinear Optics

Local Modes

Any arbitrary field may be described through the linear combinations of an orthogonal set. The concept of local modes extends the power and validity of this technique to \(z\)-dependent waveguide structures. Each point along a waveguide may be represented by a distinct set that smoothly transforms from one point to the next and that locally satisfies Maxwell’s equations, i.e., local modes are the modes of a waveguide’s instantaneous cross section. Such a generalized representation allows the modal expansion to distinguish the geometric properties of the waveguide from the underlying physics.

Local modes satisfy a modified form of the source-free, linear, and lossless Maxwell’s equations where all \(z\) derivatives have been replaced by the propagation-constant eigenvalues:

\[\begin{split}\begin{align} 0 &= \left(\nabla_t - i \, \beta_n \, \hat{\mathbf{z}}\right) \cdot \left(\epsilon_0 \, \epsilon^{\prime} \cdot \hat{\mathbf{e}}_n\right) \\ 0 &= \left(\nabla_t - i \, \beta_n \, \hat{\mathbf{z}}\right) \cdot \hat{\mathbf{h}}_n \\ \mathbf{0} &= \left(\nabla_t - i \, \beta_n \, \hat{\mathbf{z}}\right) \times \hat{\mathbf{e}}_n + i \, \omega \, \mu_0 \, \hat{\mathbf{h}}_n \\ \mathbf{0} &= \left(\nabla_t - i \, \beta_n \, \hat{\mathbf{z}}\right) \times \hat{\mathbf{h}}_n - i \, \omega \, \epsilon_0 \, \epsilon^{\prime} \cdot {\hat{\mathbf{e}}_n} \end{align}\end{split}\]

where \(\epsilon^{\prime}\) represents the relative permittivity with loss removed.

Multi-Mode Propagation Equations

As shown in the following sections, the evolution of the complex spectral amplitude 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 current and nonlinearity} \end{cases} \end{gather}\end{split}\]

Starting from Maxwell’s equations in the frequency-space domain \(\left(\omega, \mathbf{x}\right)\),

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

generalized unidirectional propagation equations (UPE) may be derived by projecting the curl equations onto the fields of a single mode (see the supplemental for more details):

\[\begin{split}\begin{align} \frac{\partial}{\partial z} a_n\!\left[\nu\right] &= - i \, \beta_n \, a_n \\ & \quad - \sum_{m} a_m \iint_{A_\infty} \frac{1}{2} \left( i \, \omega \, \epsilon_0 \, \hat{\mathbf{e}}_n^* \cdot \left(\epsilon - \epsilon^{\prime}\right) \cdot \hat{\mathbf{e}}_m \right) \, dA \\ & \quad - \sum_{m} a_m \iint_{A_\infty} \frac{1}{2} \, \hat{\mathbf{z}} \cdot \left( \hat{\mathbf{e}}_n^* \times \frac{\partial \hat{\mathbf{h}}_m}{\partial z} + \frac{\partial \hat{\mathbf{e}}_m}{\partial z} \times \hat{\mathbf{h}}^*_n \right) dA \\ & \quad -\iint_{A_\infty} \frac{1}{2} \left(\hat{\mathbf{e}}_n^* \cdot \mathbf{J}_{f} + i \, \omega \, \hat{\mathbf{e}}_n^* \cdot \mathbf{P}_\text{NL}\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 and sets the local evolution of the phase:

\[\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 one of the most important factors 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} a_m \iint_{A_\infty} \frac{1}{2} \left( i \ \omega \ \epsilon_0 \ {\hat{\mathbf{e}}_n^*} \cdot \left(\epsilon - \epsilon^{\prime}\right) \cdot {\hat{\mathbf{e}}_m} \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 \cdot \left(\epsilon - \epsilon^\prime\right) \cdot \hat{\mathbf{e}}_m\right) &= -\frac{1}{2} \left(i \, \omega \, \epsilon_0 \, \hat{\mathbf{e}}^*_n \cdot \left(\left(n + i \frac{c}{\omega} \frac{\alpha}{2} \right)^{2} - \left(n\right)^{2}\right) \cdot \hat{\mathbf{e}}_m\right) \\ &= -\frac{1}{2} \left(i \, \omega \, \epsilon_0 \, \hat{\mathbf{e}}^*_n \cdot \left(i \, n \frac{c}{\omega} \, \alpha - \frac{1}{4} \frac{c^2}{\omega^2} \, \alpha^2\right) \cdot \hat{\mathbf{e}}_m\right) \\ &= \epsilon_0 \, c \, n_\text{eff} \ \hat{\mathbf{e}}^*_n \cdot \left(\frac{n}{n_\text{eff}}\frac{\alpha}{2} + \frac{i}{8} \frac{\alpha^2}{\beta_n}\right) \cdot \hat{\mathbf{e}}_m \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 indicates that changes to the relative permittivity contribute to both 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} a_m \left(\frac{\alpha_\text{eff}}{2} + i \, \beta_\alpha\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 \cdot \left(\frac{n}{n_\text{eff}} \alpha\right) \cdot \hat{\mathbf{e}}_m \, dA \\ {\beta_\alpha}_{nm}\!\left[\nu\right] &= \iint_{A_\infty} \left(\epsilon_0 \, c \, n_\text{eff}\right) \, \hat{\mathbf{e}}^*_n \cdot \left(\frac{1}{8} \frac{\alpha^2}{\beta_n}\right) \cdot \hat{\mathbf{e}}_m \, dA \end{align}\end{split}\]

\(\beta_\alpha\) is of the order \(\beta_n^2/\alpha^2\) times smaller than \(\beta_n\). Since \(\beta_n\) is on the order of \(10^{6}\) per meter for near-infrared light, \(\beta_\alpha\) can usually be ignored. However, it may need to 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, \(\epsilon_0 \, c \, n_\text{eff} \, \hat{\mathbf{e}}_n^* \cdot \hat{\mathbf{e}}_m\) only integrates to \(\delta_{nm}\) if the longitudinal components of the electric fields are negligible. Thus, gain and loss can be expected to couple modes together if \(\alpha\) is inhomogeneous and if the \(\hat{e}_z\) components are large.

Changing Mode Profiles

\[\frac{\partial}{\partial z} a_n\!\left[\nu\right] = \ldots - \sum_{m} a_m \iint_{A_\infty} \frac{1}{2} \, \hat{\mathbf{z}} \cdot \left( \hat{\mathbf{e}}_n^* \times \frac{\partial \hat{\mathbf{h}}_m}{\partial z} + \frac{\partial \hat{\mathbf{e}}_m}{\partial z} \times \hat{\mathbf{h}}^*_n \right) dA\]

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. 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 Current 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^* \cdot \mathbf{J}_{f} \, dA & \quad \text{free current} \\ & - \iint_{A_\infty} \frac{1}{2} i \, \omega \, \hat{\mathbf{e}}_n^* \cdot \mathbf{P}_\text{NL} \, 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 typically negligible if the material is not a plasma.

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] &= \iint_{-\infty}^{\infty} \epsilon_0 \, \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] &= \iint_{-\infty}^{\infty} \epsilon_0 \, \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} -\iint_{A_\infty} \frac{1}{2} i \, \omega \, \hat{\mathbf{e}}^*_n \cdot \mathbf{P}_\text{NL}^{\left(2\right)} \ 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^*}_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] &= \iiint_{-\infty}^{\infty} \epsilon_0 \, \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] &= \iiint_{-\infty}^{\infty} \epsilon_0 \, \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} -\iint_{A_\infty} \frac{1}{2} i \, \omega \, \hat{\mathbf{e}}^*_n \cdot \mathbf{P}_\text{NL}^{\left(3\right)} \, 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^*}_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 just 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}}^* \cdot \left(\epsilon - \epsilon^{\prime} \right) \cdot \hat{\mathbf{e}}\right) dA & \quad \text{gain or loss} \\ & \quad - \iint_{A_\infty} \frac{1}{2} \left( i \, \omega \, \hat{\mathbf{e}}^* \cdot \mathbf{P}_\text{NL} \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}}^* \cdot \left(\epsilon - \epsilon^{\prime} \right) \cdot \hat{\mathbf{e}}\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}}^* \cdot \left(\frac{n}{n_\text{eff}} \alpha\right) \cdot \hat{\mathbf{e}} \ dA \\ \beta_{\alpha}\!\left[\nu\right] &= \iint_{A_\infty} \left( \epsilon_0 \, c \, n_\text{eff} \right) \, \hat{\mathbf{e}}^* \cdot \left(\frac{1}{8} \frac{\alpha^2}{\beta_\text{eff}}\right) \cdot \hat{\mathbf{e}} \ 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 given 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] = \iint_{A_\infty} \frac{1}{2} \epsilon_0 \, \chi^{\left(2\right)}_{i j k}\!\left[\nu_1, \nu_2\right] \hat{\mathbf{e}}^*_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}}^*_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- or cross-phase modulation is always phase matched and corresponds to the nonlinear interaction between positive and negative frequency components from a pair of modes \(\left(\omega = \omega + \omega^\prime - \omega^\prime\right)\). In this case the effective nonlinearity collapses to a function of two variables:

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

If one ignores mode overlap and assumes that a particular pump frequency (\(\nu^\prime\)) dominates, 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}\!\left[\nu, \nu\right] \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}\]