---
title: "Modulational instability"
sort_title: "Modulational instability"
date: 2021-02-26
categories:
- Physics
- Fiber optics
- Optics
- Perturbation
- Nonlinear optics
layout: "concept"
---

In fiber optics, **modulational instability** (MI)
is a nonlinear effect that leads to the exponential amplification
of background noise in certain frequency regions.
It only occurs in the [anomalous dispersion regime](/know/concept/dispersive-broadening/)
($$\beta_2 < 0$$), which we will prove shortly.
The underlying physical process causing it is *degenerate four-wave mixing*.

Consider the following simple solution to the nonlinear Schrödinger equation:
a time-invariant constant power $$P_0$$ at the carrier frequency $$\omega_0$$,
experiencing [self-phase modulation](/know/concept/self-phase-modulation/):

$$\begin{aligned}
    A(z,t)
    = \sqrt{P_0} \exp( i \gamma P_0 z)
\end{aligned}$$

We add a small perturbation $$\varepsilon(z,t)$$ to this signal,
representing background noise:

$$\begin{aligned}
    A(z,t)
    = \big(\sqrt{P_0} + \varepsilon(z,t)\big) \exp( i \gamma P_0 z)
\end{aligned}$$

We insert this into the nonlinear Schrödinger equation to get a perturbation equation,
which we linearize by assuming that $$|\varepsilon|^2$$ is negligible compared to $$P_0$$,
such that all higher-order terms of $$\varepsilon$$ can be dropped, leaving:

$$\begin{aligned}
    0
    &= - P_0 \sqrt{P_0} \gamma - P_0 \gamma \varepsilon + i \pdv{\varepsilon}{z}
    - \frac{\beta_2}{2} \pdvn{2}{\varepsilon}{t}
    + \gamma \big(\sqrt{P_0} + \varepsilon\big)^2 \big(\sqrt{P_0} + \varepsilon\big)^*
    \\
    &= i \pdv{\varepsilon}{z}
    - \frac{\beta_2}{2} \pdvn{2}{\varepsilon}{t}
    + \gamma \big( P_0 (\varepsilon + \varepsilon^*) + \sqrt{P_0} |\varepsilon|^2
    + \sqrt{P_0} \varepsilon (\varepsilon + \varepsilon^*) + \varepsilon |\varepsilon|^2 \big)
    \\
    &= i \pdv{\varepsilon}{z} - \frac{\beta_2}{2} \pdvn{2}{\varepsilon}{t} + \gamma P_0 (\varepsilon + \varepsilon^*)
\end{aligned}$$

We split the perturbation into real and imaginary parts
$$\varepsilon(z,t) = \varepsilon_r(z,t) + i \varepsilon_i(z,t)$$,
which we put in this equation.
The point is that $$\varepsilon_r$$ and $$\varepsilon_i$$ are real functions:

$$\begin{aligned}
    0
    &= i \pdv{\varepsilon_r}{z} - \pdv{\varepsilon_i}{z}
    - \frac{\beta_2}{2} \pdvn{2}{\varepsilon_r}{t} - i \frac{\beta_2}{2} \pdvn{2}{\varepsilon_i}{t}
    + 2 \gamma P_0 \varepsilon_r
\end{aligned}$$

Splitting this into its real and imaginary parts gives two PDEs
relating $$\varepsilon_r$$ and $$\varepsilon_i$$:

$$\begin{aligned}
    \pdv{\varepsilon_r}{z}
    \frac{\beta_2}{2} \pdvn{2}{\varepsilon_i}{t}
    \qquad \qquad
    \pdv{\varepsilon_i}{z}
    = - \frac{\beta_2}{2} \pdvn{2}{\varepsilon_r}{t} + 2 \gamma P_0 \varepsilon_r
\end{aligned}$$

We [Fourier transform](/know/concept/fourier-transform/)
these in $$t$$ to turn them into ODEs relating
$$\tilde{\varepsilon}_r(z,\omega)$$ and $$\tilde{\varepsilon}_i(z,\omega)$$:

$$\begin{aligned}
    \pdv{\tilde{\varepsilon}_r}{z}
    = - \frac{\beta_2}{2} \omega^2 \tilde{\varepsilon}_i
    \qquad \qquad
    \pdv{\tilde{\varepsilon}_i}{z}
    = \Big(\frac{\beta_2}{2} \omega^2 + 2 \gamma P_0 \Big) \tilde{\varepsilon}_r
\end{aligned}$$

We are interested in exponential growth, so let us make the following ansatz,
where $$k$$ may be a function of $$\omega$$, as long as it is $$z$$-invariant:

$$\begin{aligned}
    \tilde{\varepsilon}_r(z, \omega)
    = \tilde{\varepsilon}_r(0, \omega) \exp(k z)
    \qquad \qquad
    \tilde{\varepsilon}_i(z, \omega)
    = \tilde{\varepsilon}_i(0, \omega) \exp(k z)
\end{aligned}$$

With this, we can write the system of ODEs for
$$\tilde{\varepsilon}_r(z,\omega)$$ and $$\tilde{\varepsilon}_i(z,\omega)$$
in matrix form:

$$\begin{aligned}
    \begin{bmatrix}
        k & \beta_2 \omega^2 / 2 \\
        \beta_2 \omega^2 / 2 \!+\! 2 \gamma P_0 & - k
    \end{bmatrix}
    \cdot
    \begin{bmatrix} \tilde{\varepsilon}_r(0, \omega) \\ \tilde{\varepsilon}_i(0, \omega) \end{bmatrix}
    =
    \begin{bmatrix} 0 \\ 0 \end{bmatrix}
\end{aligned}$$

This has nonzero solutions if the system matrix' determinant is zero,
which is true when:

$$\begin{aligned}
    k
    = \pm \sqrt{ - \frac{\beta_2}{2} \omega^2 \Big( \frac{\beta_2}{2} \omega^2 + 2 \gamma P_0 \Big) }
\end{aligned}$$

To get exponential growth, it is essential that $$\mathrm{Re}\{k\} > 0$$,
so we discard the negative sign,
and get the following condition for MI:

$$\begin{aligned}
    -\frac{\beta_2}{2} \omega^2 \Big( \frac{\beta_2}{2} \omega^2 + 2 \gamma P_0 \Big)
    > 0
    \qquad \implies \qquad
    \boxed{
        \omega^2
        < -\frac{4 \gamma P_0}{\beta_2}
    }
\end{aligned}$$

Since $$\omega^2$$ is positive, MI can only occur when $$\beta_2$$ is negative.
It is worth noting that $$\beta_2 = \beta_2(\omega_0)$$,
meaning there can only be exponential noise growth
when the parent pulse is in the anomalous dispersion regime,
but that growth may appear in areas of normal dispersion,
as long as the above condition is satisfied by the parent.

This result has been derived using perturbation,
so only holds as long as $$|\varepsilon|^2 \ll P_0$$.
Over time, the noise gets amplified so greatly
that this approximation breaks down.

Next, we define the **gain** $$g(\omega)$$,
which expresses how quickly the
perturbation grows as a function of the frequency offset $$\omega$$:

$$\begin{aligned}
    \boxed{
        g(\omega)
        = \mathrm{Re}\{k\}
        = \mathrm{Re} \bigg\{ \sqrt{ - \frac{\beta_2}{2} \omega^2 \Big( \frac{\beta_2}{2} \omega^2 + 2 \gamma P_0 \Big) } \bigg\}
    }
\end{aligned}$$

The frequencies with maximum gain are then found as extrema of $$g(\omega)$$,
which satisfy:

$$\begin{aligned}
    g'(\omega_\mathrm{max})
    = 0
    \qquad \implies \qquad
    \boxed{
        \omega_\mathrm{max}
        = \pm \sqrt{\frac{2 \gamma P_0}{-\beta_2}}
    }
\end{aligned}$$

A simulation of MI is illustrated below.
The pulse considered was a soliton of the following form
with settings $$T_0 = 10\:\mathrm{ps}$$, $$P_0 = 10\:\mathrm{kW}$$,
$$\beta = -10\:\mathrm{ps}^2/\mathrm{m}$$ and $$\gamma = 0.1/\mathrm{W}/\mathrm{m}$$,
whose peak is approximately flat, so our derivation is valid there,
hence it "wrinkles" in the $$t$$-domain:

$$\begin{aligned}
    A(0, t)
    = \sqrt{P_0} \sech\!\Big(\frac{t}{T_0}\Big)
\end{aligned}$$

{% include image.html file="simulation-full.png" width="100%"
    alt="Modulational instability simulation results" %}

Where $$L_\mathrm{NL} = 1/(\gamma P_0)$$ is the characteristic length of nonlinear effects.
Note that no noise was added to the simulation;
you are seeing pure numerical errors getting amplified.

If one of the gain peaks accumulates a lot of energy quickly ($$L_\mathrm{NL}$$ is small),
and that peak is in the anomalous dispersion regime,
then it can in turn also cause MI in its own surroundings,
leading to a cascade of secondary and tertiary gain areas.
This is seen above for $$z > 30 L_\mathrm{NL}$$.

What we described is "pure" MI, but there also exists
a different type caused by Raman scattering.
In that case, amplification occurs at the strongest peak of the Raman gain $$\tilde{g}_R(\omega)$$,
even when the parent pulse has $$\beta_2 > 0$$.
This is an example of stimulated Raman scattering (SRS).



## References
1.  O. Bang,
    *Numerical methods in photonics: lecture notes*, 2019,
    unpublished.
2.  O. Bang,
    *Nonlinear mathematical physics: lecture notes*, 2020,
    unpublished.