Categories: Fiber optics, Nonlinear optics, Optics, Perturbation, Physics.

Modulational instability

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 (β2<0\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 P0P_0 at the carrier frequency ω0\omega_0, experiencing self-phase modulation:

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

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

A(z,t)=(P0+ε(z,t))exp(iγP0z)\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 ε2|\varepsilon|^2 is negligible compared to P0P_0, such that all higher-order terms of ε\varepsilon can be dropped, leaving:

0=P0P0γP0γε+iεzβ222εt2+γ(P0+ε)2(P0+ε)=iεzβ222εt2+γ(P0(ε+ε)+P0ε2+P0ε(ε+ε)+εε2)=iεzβ222εt2+γP0(ε+ε)\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 ε(z,t)=εr(z,t)+iεi(z,t)\varepsilon(z,t) = \varepsilon_r(z,t) + i \varepsilon_i(z,t), which we put in this equation. The point is that εr\varepsilon_r and εi\varepsilon_i are real functions:

0=iεrzεizβ222εrt2iβ222εit2+2γP0εr\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 εr\varepsilon_r and εi\varepsilon_i:

εrzβ222εit2εiz=β222εrt2+2γP0εr\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 these in tt to turn them into ODEs relating ε~r(z,ω)\tilde{\varepsilon}_r(z,\omega) and ε~i(z,ω)\tilde{\varepsilon}_i(z,\omega):

ε~rz=β22ω2ε~iε~iz=(β22ω2+2γP0)ε~r\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 kk may be a function of ω\omega, as long as it is zz-invariant:

ε~r(z,ω)=ε~r(0,ω)exp(kz)ε~i(z,ω)=ε~i(0,ω)exp(kz)\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 ε~r(z,ω)\tilde{\varepsilon}_r(z,\omega) and ε~i(z,ω)\tilde{\varepsilon}_i(z,\omega) in matrix form:

[kβ2ω2/2β2ω2/2 ⁣+ ⁣2γP0k][ε~r(0,ω)ε~i(0,ω)]=[00]\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:

k=±β22ω2(β22ω2+2γP0)\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 Re{k}>0\mathrm{Re}\{k\} > 0, so we discard the negative sign, and get the following condition for MI:

β22ω2(β22ω2+2γP0)>0    ω2<4γP0β2\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 ω2\omega^2 is positive, MI can only occur when β2\beta_2 is negative. It is worth noting that β2=β2(ω0)\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 ε2P0|\varepsilon|^2 \ll P_0. Over time, the noise gets amplified so greatly that this approximation breaks down.

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

g(ω)=Re{k}=Re{β22ω2(β22ω2+2γP0)}\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(ω)g(\omega), which satisfy:

g(ωmax)=0    ωmax=±2γP0β2\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 T0=10psT_0 = 10\:\mathrm{ps}, P0=10kWP_0 = 10\:\mathrm{kW}, β=10ps2/m\beta = -10\:\mathrm{ps}^2/\mathrm{m} and γ=0.1/W/m\gamma = 0.1/\mathrm{W}/\mathrm{m}, whose peak is approximately flat, so our derivation is valid there, hence it “wrinkles” in the tt-domain:

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

Modulational instability simulation results

Where LNL=1/(γP0)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 (LNLL_\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>30LNLz > 30 L_\mathrm{NL}.

Here we described “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 g~R(ω)\tilde{g}_R(\omega), even when the parent pulse has β2>0\beta_2 > 0. This is an example of stimulated Raman scattering (SRS).


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