path: root/source/know/concept/nonlinear-schrodinger-equation/index.md

---
title: "Nonlinear Schrödinger equation"
sort_title: "Nonlinear Schrodinger equation" # sic
date: 2024-09-15
categories:
- Physics
- Mathematics
- Fiber optics
- Nonlinear optics
layout: "concept"
---

The **nonlinear Schrödinger (NLS) equation**
is a nonlinear 1+1D partial differential equation
that appears in many areas of physics.
It is often given in its dimensionless form,
where it governs the envelope $$u(z, t)$$
of an underlying carrier wave,
with $$t$$ the transverse coordinate,
and $$r = \pm 1$$ a parameter determining
which of two regimes the equation is intended for:

$$\begin{aligned}
    \boxed{
        i \pdv{u}{z} + \pdvn{2}{u}{t} + r |u|^2 u
        = 0
    }
\end{aligned}$$

Many variants exist, depending on the conventions used by authors.
The NLS equation is used to describe pulses in fiber optics (as derived below),
waves over deep water, local opening of DNA chains, and much more.
Very roughly speaking, it is a valid description of
"all" weakly nonlinear, slowly modulated waves in physics.

It exhibits an incredible range of behaviors,
from "simple" effects such as
[dispersive broadening](/know/concept/dispersive-broadening/),
[self-phase modulation](/know/concept/self-phase-modulation/)
and [first-order solitons](/know/concept/optical-soliton/),
to weirder and more complicated phenomena like
[modulational instability](/know/concept/modulational-instability/),
[optical wave breaking](/know/concept/optical-wave-breaking/)
and periodic *higher-order solitons*.
It is also often modified to include additional physics,
further enriching its results with e.g.
[self-steepening](/know/concept/self-steepening/)
and *soliton self-frequency shifting*.

We only consider fiber optics here;
the NLS equation can be derived in many other ways.
We start from the most general form of the
[electromagnetic wave equation](/know/concept/electromagnetic-wave-equation/),
after assuming the medium cannot be magnetized ($$\mu_r = 1$$):

$$\begin{aligned}
    \nabla \cross \big( \nabla \cross \vb{E} \big)
    = - \mu_0 \varepsilon_0 \pdvn{2}{\vb{E}}{t} - \mu_0 \pdvn{2}{\vb{P}}{t}
\end{aligned}$$

Using the vector identity
$$\nabla \cross (\nabla \cross \vb{E}) = \nabla (\nabla \cdot \vb{E}) - \nabla^2 \vb{E}$$
and [Gauss's law](/know/concept/maxwells-equations/) $$\nabla \cdot \vb{E} = 0$$,
and splitting the polarization $$\vb{P}$$
into linear and nonlinear contributions
$$\vb{P}_\mathrm{L}$$ and $$\vb{P}_\mathrm{NL}$$:

$$\begin{aligned}
    \nabla^2 \vb{E} - \mu_0 \varepsilon_0 \pdvn{2}{\vb{E}}{t}
    &= \mu_0 \pdvn{2}{\vb{P}_\mathrm{L}}{t} + \mu_0 \pdvn{2}{\vb{P}_\mathrm{NL}}{t}
\end{aligned}$$

In general, $$\vb{P}_\mathrm{L}$$ is given by the convolution
of $$\vb{E}$$ with a second-rank response tensor $$\chi^{(1)}$$:

$$\begin{aligned}
    \vb{P}_\mathrm{L}(\vb{r}, t)
    = \varepsilon_0 \int_{-\infty}^\infty \chi^{(1)}(t - t') \cdot \vb{E}(\vb{r}, t') \dd{t'}
\end{aligned}$$

In $$\vb{P}_\mathrm{NL}$$ we only include third-order nonlinearities,
since higher orders are usually negligible,
and second-order nonlinear effects only exist in very specific crystals.
So we "only" need to deal with a fourth-rank response tensor $$\chi^{(3)}$$:

$$\begin{aligned}
    \vb{P}_\mathrm{NL}(\vb{r}, t)
    = \varepsilon_0 \iiint_{-\infty}^\infty \chi^{(3)}(t \!-\! t_1, t \!-\! t_2, t \!-\! t_3)
    \:\vdots\: \vb{E}(\vb{r}, t_1) \vb{E}(\vb{r}, t_2) \vb{E}(\vb{r}, t_3) \dd{t_1} \dd{t_2} \dd{t_3}
\end{aligned}$$

In practice, two phenomena contribute to $$\chi^{(3)}$$:
the *Kerr effect* due to electrons' response to $$\vb{E}$$,
and *Raman scattering* due to nuclei's response,
which is slower because of their mass.
But if the light pulses are sufficiently long (>1ps in silica),
both effects can be treated as fast, so:

$$\begin{aligned}
    \chi^{(3)}(t \!-\! t_1, t \!-\! t_2, t \!-\! t_3)
    &= \chi^{(3)} \delta(t - t_1) \delta(t - t_2) \delta(t - t_3)
\end{aligned}$$

Where $$\delta$$ is the [Dirac delta function](/know/concept/dirac-delta-function/).
To keep things simple,
we consider linearly $$x$$-polarized light $$\vb{E} = \vu{x} |\vb{E}|$$,
such that the tensor can be replaced with its scalar element $$\chi^{(3)}_{xxxx}$$.
Then:

$$\begin{aligned}
    \vb{P}_\mathrm{NL}
    = \varepsilon_0 \chi^{(3)}_{xxxx} \big( \vb{E} \cdot \vb{E} \big) \vb{E}
\end{aligned}$$

For the same reasons, the linear polarization is reduced to:

$$\begin{aligned}
    \vb{P}_\mathrm{L}
    &= \varepsilon_0 \chi^{(1)}_{xx} \vb{E}
\end{aligned}$$

Next, we decompose $$\vb{E}$$ as follows,
consisting of a carrier wave $$e^{-i \omega_0 t}$$
at a constant frequency $$\omega_0$$,
modulated by an envelope $$E$$
that is assumed to be slowly-varying compared to the carrier,
plus the complex conjugate $$E^* e^{i \omega_0 t}$$:

$$\begin{aligned}
    \vb{E}(\vb{r}, t)
    &= \vu{x} \frac{1}{2} \Big( E(\vb{r}, t) e^{- i \omega_0 t} + E^*(\vb{r}, t) e^{i \omega_0 t} \Big)
\end{aligned}$$

Note that no generality has been lost in this step.
Inserting it into the polarizations:

$$\begin{aligned}
    \mathrm{P}_\mathrm{L}
    &= \vu{x} \frac{1}{2} \varepsilon_0 \chi^{(1)}_{xx} \Big( E(\vb{r}, t) e^{- i \omega_0 t} + E^*(\vb{r}, t) e^{i \omega_0 t} \Big)
    \\
    \vb{P}_\mathrm{NL}
    &= \vu{x} \frac{1}{8} \varepsilon_0 \chi^{(3)}_{xxxx} \Big( E e^{- i \omega_0 t} + E^* e^{i \omega_0 t} \Big)^{3}
    \\
    &= \vu{x} \frac{1}{8} \varepsilon_0 \chi^{(3)}_{xxxx}
    \Big( E^3 e^{- i 3 \omega_0 t} + 3 E^2 E^* e^{- i \omega_0 t} + 3 E (E^*)^2 e^{i \omega_0 t} + (E^*)^3 e^{i 3 \omega_0 t} \Big)
\end{aligned}$$

The terms with $$3 \omega_0$$ represent *third-harmonic generation*,
and only matter if the carrier is phase-matched
to the tripled wave, which is generally not the case,
so they can be ignored.
Now, if we decompose the polarizations in the same was as $$\vb{E}$$:

$$\begin{aligned}
    \vb{P}_\mathrm{L}(\vb{r}, t)
    &= \vu{x} \frac{1}{2} \Big( P_\mathrm{L}(\vb{r}, t) e^{- i \omega_0 t} + P_\mathrm{L}^*(\vb{r}, t) e^{i \omega_0 t} \Big)
    \\
    \vb{P}_\mathrm{NL}(\vb{r}, t)
    &= \vu{x} \frac{1}{2} \Big( P_\mathrm{NL}(\vb{r}, t) e^{- i \omega_0 t} + P_\mathrm{NL}^*(\vb{r}, t) e^{i \omega_0 t} \Big)
\end{aligned}$$

Then it is straightforward to see that their envelope functions are given by:

$$\begin{aligned}
    P_\mathrm{L}
    &= \varepsilon_0 \chi^{(1)}_{xx} E
    \\
    P_\mathrm{NL}
    &= \frac{3}{4} \varepsilon_0 \chi^{(3)}_{xxxx} |E|^2 E
\end{aligned}$$

The forward carrier $$e^{- i \omega_0 t}$$
and the backward carrier $$e^{i \omega_0 t}$$
can be regarded as separate channels,
which only interact via $$P_\mathrm{NL}$$.
From now on, we only consider the forward-propagating wave,
so all terms containing $$e^{i \omega_0 t}$$ are dropped;
by taking the complex conjugate of the resulting equations,
the backward-propagating counterparts can always be recovered,
so no information is really lost.
Therefore, the main wave equation becomes:

$$\begin{aligned}
    0
    &= \bigg(
        \nabla^2 E - \mu_0 \varepsilon_0 \pdvn{2}{E}{t} - \mu_0 \pdvn{2}{P_\mathrm{L}}{t} - \mu_0 \pdvn{2}{P_\mathrm{NL}}{t}
    \bigg) e^{-i \omega_0 t}
    \\
    &\approx \bigg(
        \nabla^2 E - \Big( 1 + \chi^{(1)}_{xx} + \frac{3}{4} \chi^{(3)}_{xxxx} |E|^2 \Big) \mu_0 \varepsilon_0 \pdvn{2}{E}{t}
    \bigg) e^{-i \omega_0 t}
\end{aligned}$$

Where we have used our assumption that $$E$$ is slowly-varying
to treat $$|E|^2$$ as a constant,
in order to move it outside the $$t$$-derivative.
We thus arrive at:

$$\begin{aligned}
    0
    &= \bigg( \nabla^2 E - \frac{\varepsilon_r}{c^2} \pdvn{2}{E}{t} \bigg) e^{-i \omega_0 t}
\end{aligned}$$

Where $$c = 1 / \sqrt{\mu_0 \varepsilon_0}$$ is the phase velocity of light in a vacuum,
and the relative permittivity $$\varepsilon_r$$ is defined as shown below.
Note that this is a mild abuse of notation,
since the symbol $$\varepsilon_r$$ is usually reserved for linear materials:

$$\begin{aligned}
    \varepsilon_r
    \equiv 1 + \chi^{(1)}_{xx} + \frac{3}{4} \chi^{(3)}_{xxxx} |E|^2
\end{aligned}$$

Next, we take the [Fourier transform](/know/concept/fourier-transform/)
$$t \to \omega$$ of the wave equation,
again treating $$|E|^2$$ (inside $$\varepsilon_r$$) as a constant.
The constant $$s = \pm 1$$ is included here
to deal with the fact that different authors use different sign conventions:

$$\begin{aligned}
    0
    &= \hat{\mathcal{F}}\bigg\{ \bigg( \nabla^2 E - \frac{\varepsilon_r}{c^2} \pdvn{2}{E}{t} \bigg) e^{-i \omega_0 t} \bigg\}
    \\
    &= \int_{-\infty}^\infty
    \bigg( \nabla^2 E - \frac{\varepsilon_r}{c^2} \pdvn{2}{E}{t} \bigg)
    e^{i s (\omega - \omega_0) t} \dd{t}
    \\
    &= \nabla^2 E + s^2 (\omega - \omega_0)^2 \frac{\varepsilon_r}{c^2} E
\end{aligned}$$

We use $$s^2 = 1$$ and define $$\Omega \equiv \omega - \omega_0$$
as the frequency shift relative to the carrier wave:

$$\begin{aligned}
    0
    &= \nabla^2 E + \frac{\Omega^2 \varepsilon_r}{c^2} E
\end{aligned}$$

This is a so-called *Helmholtz equation* in 3D,
which we will solve using separation of variables,
by assuming that its solution can be written as:

$$\begin{aligned}
    E(\vb{r}, \Omega)
    &= F(x, y) \: A(z, \Omega) \: e^{i \beta_0 z}
\end{aligned}$$

Where $$\beta_0$$ is the wavenumber of the carrier,
which will be determined later.
Inserting this ansatz into the Helmholtz equation yields:

$$\begin{aligned}
    0
    &= \Big( \pdvn{2}{F}{x} + \pdvn{2}{F}{y} \Big) A e^{i \beta_0 z}
    + \pdvn{2}{}{z} \Big( A e^{i \beta_0 z} \Big) F
    + \frac{\Omega^2 \varepsilon_r}{c^2} F A e^{i \beta_0 z}
    \\
    &= \bigg( \Big( \pdvn{2}{F}{x} + \pdvn{2}{F}{y} \Big) A
    + \Big( \pdvn{2}{A}{z} + 2 i \beta_0 \pdv{A}{z} - \beta_0^2 A \Big) F
    + \frac{\Omega^2 \varepsilon_r}{c^2} F A \bigg) e^{i \beta_0 z}
\end{aligned}$$

We divide by $$F A \: e^{i \beta_0 z}$$
and rearrange the terms in a specific way:

$$\begin{aligned}
    \Big( \pdvn{2}{F}{x} + \pdvn{2}{F}{y} \Big) \frac{1}{F} + \frac{\Omega^2 \varepsilon_r}{c^2}
    &= - 2 i \beta_0 \pdv{A}{z} \frac{1}{A} + \beta_0^2
\end{aligned}$$

Now all the $$x$$- and $$y$$-dependence is on the left,
and the $$z$$-dependence is on the right.
We have placed the $$\varepsilon_r$$-term on the left too
because it depends relatively strongly on $$(x, y)$$
to describe the fiber's internal structure,
and weakly on $$z$$ due to nonlinear effects.
Meanwhile, $$\beta_0$$ is on the right because that will lead to
a nicer equation for $$A$$ later.

Note that both sides are functions of $$\Omega$$.
Based on the aforementioned dependences,
in order for this equation to have a solution for all $$(x, y, z)$$,
there must exist a quantity $$\beta(\Omega)$$ that is constant in space,
such that we obtain two separated equations for $$F$$ and $$A$$:

$$\begin{aligned}
    \beta(\omega)
    &= \bigg( \pdvn{2}{F}{x} + \pdvn{2}{F}{y} \bigg) \frac{1}{F} + \frac{\omega^2 \varepsilon_r}{c^2}
    \\
    \beta(\Omega)
    &= - 2 i \beta_0 \pdv{A}{z} \frac{1}{A} + \beta_0^2
\end{aligned}$$

Note that we replaced $$\Omega$$ with $$\omega$$ in $$F$$'s equation
(and redefined $$\beta$$ and $$\varepsilon_r$$ accordingly).
This is not an innocent detail:
the idea is that $$\omega \sqrt{\varepsilon_r} / c$$
would be the light's wavenumber if it had not been trapped in a waveguide,
and that $$\beta$$ is the *confined* wavenumber,
also known as the **propagation constant**.
If we had kept $$\Omega$$,
the meaning of $$\beta$$ would not be so straightforward.

The difference between $$\beta(\omega)$$ and $$\beta_0$$
is simply that $$\beta_0 \equiv \beta(\omega_0)$$.
Our ansatz for separating the variables contained $$\beta_0$$,
such that the full carrier wave $$e^{i \beta_0 z - i \omega_0 t}$$ was represented
(with $$e^{- i \omega_0 t}$$ now hidden inside the Fourier transform).
But later, to properly describe how light behaves inside the fiber,
the full dispersion relation $$\beta(\omega)$$ will be needed.

Multiplying by $$F$$ and $$A$$,
we get the following set of equations,
implicitly coupled via $$\beta$$:

$$\begin{aligned}
    \boxed{
        \begin{aligned}
            0
            &= \pdvn{2}{F}{x} + \pdvn{2}{F}{y} + \bigg( \frac{\omega^2 \varepsilon_r}{c^2} - \beta^2 \bigg) F
            \\
            0
            &= 2 i \beta_0 \pdv{A}{z} + \big( \beta^2 - \beta_0^2 \big) A
        \end{aligned}
    }
\end{aligned}$$

The equation for $$F$$ must be solved first.
To do so, we treat the nonlinearity as a perturbation
to be neglected initially.
In other words, we first solve the following eigenvalue problem for $$\beta^2$$,
where $$n(x, y)$$ is the linear refractive index,
with $$n^2 = 1 + \Real\{\chi^{(1)}_{xx}\} \approx \varepsilon_r$$:

$$\begin{aligned}
    \pdvn{2}{F}{x} + \pdvn{2}{F}{y} + \bigg( \frac{\omega^2 n^2}{c^2} - \beta^2 \bigg) F
    = 0
\end{aligned}$$

This gives us the allowed values of $$\beta$$;
see [step-index fiber](/know/concept/step-index-fiber/) for an example solution.
Now we add the small index change $$\Delta{n}(x, y)$$ due to nonlinear effects:

$$\begin{aligned}
    \varepsilon_r
    = (n + \Delta{n})^2
    \approx n^2 + 2 n \: \Delta{n}
\end{aligned}$$

Then it can be shown using first-order
[perturbation theory](/know/concept/time-independent-perturbation-theory/)
that the eigenfunction $$F$$ is not really affected,
and the eigenvalue $$\beta^2$$ is shifted by $$\Delta(\beta^2)$$, given by:

$$\begin{aligned}
    \Delta(\beta^2)
    = \frac{2 \omega^2}{c^2} \frac{\displaystyle \iint_{-\infty}^\infty n \: \Delta{n} \: |F|^2 \dd{x} \dd{y}}
    {\displaystyle \iint_{-\infty}^\infty |F|^2 \dd{x} \dd{y}}
\end{aligned}$$

But we are more interested in the *wavenumber* shift $$\Delta{\beta}$$
than the *eigenvalue* shift $$\Delta(\beta^2)$$.
They are related to one another as follows:

$$\begin{aligned}
    \beta^2 + \Delta(\beta^2)
    = (\beta + \Delta{\beta})^2
    \approx \beta^2 + 2 \beta \Delta{\beta}
\end{aligned}$$

Furthermore, we assume that the fiber only consists of materials
with similar refractive indices, or in other words,
that it confines the light using only a small index difference,
in which case we can treat $$n$$ as a constant and move it outside the integral.
Then $$\Delta{\beta}$$ becomes:

$$\begin{aligned}
    \Delta{\beta}
    = \frac{\omega^2 n}{\beta c^2} \frac{\displaystyle \iint_{-\infty}^\infty \Delta{n} \: |F|^2 \dd{x} \dd{y}}
    {\displaystyle \iint_{-\infty}^\infty |F|^2 \dd{x} \dd{y}}
\end{aligned}$$

Recall that $$\beta$$ is the wavenumber of the confined mode:
by solving the unperturbed $$F$$-equation,
it can be shown that $$\beta$$'s value is somewhere
between the bulk wavenumbers of the fiber materials.
Since we just approximated $$n$$ as a constant,
this means that $$\omega n / c \approx \beta$$, leading us to
the general "final" form of $$\Delta{\beta}$$,
with all the arguments shown for clarity:

$$\begin{aligned}
    \boxed{
        \Delta{\beta}(\omega)
        = \frac{\omega}{c \mathcal{A}_\mathrm{mode}}
        \iint_{-\infty}^\infty \Delta{n}(x, y, \omega) \: |F(x, y)|^2 \dd{x} \dd{y}
    }
\end{aligned}$$

Where we have defined the *mode area* $$\mathcal{A}_\mathrm{mode}$$ as shown below.
In order for $$\mathcal{A}_\mathrm{mode}$$ to be in units of area,
$$F$$ must be dimensionless,
and consequently $$A$$ has (SI) units of an electric field.

$$\begin{aligned}
    \mathcal{A}_\mathrm{mode}
    \equiv \iint_{-\infty}^\infty |F|^2 \dd{x} \dd{y}
\end{aligned}$$

Now we finally turn our attention to the equation for $$A$$.
Before perturbation, it was:

$$\begin{aligned}
    0
    &= 2 i \beta_0 \pdv{A}{z} + \big( \beta^2 - \beta_0^2 \big) A
\end{aligned}$$

Where $$\beta \approx \beta_0$$, so we can replace
$$\beta^2 - \beta_0^2$$ with $$2 \beta_0 (\beta - \beta_0)$$.
Also including $$\Delta{\beta}$$, we get:

$$\begin{aligned}
    0
    &= i \pdv{A}{z} + \big( \beta + \Delta{\beta} - \beta_0 \big) A
\end{aligned}$$

Usually, we do not know a full expression for $$\beta(\omega)$$,
so it makes sense to expand it around the carrier frequency $$\omega_0$$ as follows,
where $$\beta_n = \idvn{n}{\beta}{\omega} |_{\omega = \omega_0}$$:

$$\begin{aligned}
    \beta(\omega)
    &= \beta_0
    + (\omega - \omega_0) \beta_1
    + (\omega - \omega_0)^2 \frac{\beta_2}{2}
    + (\omega - \omega_0)^3 \frac{\beta_3}{6}
    + \: ...
\end{aligned}$$

Spectrally, the broader the light pulse, the more terms must be included.
Recall that earlier, in order to treat $$\chi^{(3)}$$ as instantaneous,
we already assumed a temporally broad
(spectrally narrow) pulse.
Hence, for simplicity, we can cut off this Taylor series at $$\beta_2$$,
which is good enough in many cases.
Inserting the expansion into $$A$$'s equation:

$$\begin{aligned}
    0
    &= i \pdv{A}{z} + i \frac{\beta_1}{s} (-i s \Omega) A - \frac{\beta_2}{2 s^2} (- i s \Omega)^2 A + \Delta{\beta}_0 A
\end{aligned}$$

Which we have rewritten in preparation for taking the inverse Fourier transform,
by introducing $$s$$ and by replacing $$\Delta{\beta}(\omega)$$
with $$\Delta{\beta_0} \equiv \Delta{\beta}(\omega_0)$$
in order to remove all explicit dependence on $$\omega$$,
i.e. we only keep the first term of $$\Delta{\beta}$$'s Taylor expansion.
After transforming and using $$s^2 = 1$$,
we get the following equation for $$A(z, t)$$:

$$\begin{aligned}
    0
    &= i \pdv{A}{z} + i s \beta_1 \pdv{A}{t} - \frac{\beta_2}{2} \pdvn{2}{A}{t} + \Delta{\beta}_0 A
\end{aligned}$$

The next step is to insert our expression for $$\Delta{\beta}_0$$,
for which we must first choose a specific form for $$\Delta{n}$$
according to which effects we want to include.
Earlier, we approximated $$\varepsilon_r \approx n^2$$,
so if we instead say that $$\varepsilon_r = (n \!+\! \Delta{n})^2$$,
then $$\Delta{n}$$ should include absorption and nonlinearity.
The most commonly used form for $$\Delta{n}$$ is therefore:

$$\begin{aligned}
    \Delta{n}(x, y, \omega)
    = n_2(\omega) \: I(x, y, \omega) + i \frac{c \alpha(\omega)}{2 \omega}
\end{aligned}$$

Where $$I$$ is the intensity (i.e. power per unit area) of the light,
$$n_2$$ is the material's *Kerr coefficient* in units of inverse intensity,
and $$\alpha$$ is the attenuation coefficient
consisting of linear and nonlinear contributions
(see [multi-photon absorption](/know/concept/multi-photon-absorption/)).
Specifically, they are given by:

$$\begin{aligned}
    n_2
    = \frac{3 \Real\{\chi^{(3)}_{xxxx}\}}{4 \varepsilon_0 c n^2}
    \qquad
    \alpha
    = \frac{\omega \Imag\{\chi^{(1)}_{xx}\}}{c n}
    + \frac{3 \omega \Imag\{\chi^{(3)}_{xxxx}\}}{2 \varepsilon_0 c^2 n^2} I
    \qquad
    I
    = \frac{\varepsilon_0 c n}{2} |F|^2 |A|^2
\end{aligned}$$

For simplicity we set $$\Imag\{\chi^{(3)}_{xxxx}\} = 0$$,
which is a good approximation for silica fibers.
Inserting this form of $$\Delta{n}$$ into $$\Delta{\beta_0}$$
and neglecting the $$(x, y)$$-dependence of $$\Delta{n}$$ yields:

$$\begin{aligned}
    \Delta{\beta}_0
    &= i \frac{\alpha}{2} \frac{\mathcal{A}_\mathrm{mode}}{\mathcal{A}_\mathrm{mode}}
    + \frac{\omega_0 \varepsilon_0 c n n_2}{2 c \math