---
title: "Self-steepening"
sort_title: "Self-steepening"
date: 2024-09-29 # Originally 2021-02-26, major rewrite
categories:
- Physics
- Optics
- Fiber optics
- Nonlinear optics
layout: "concept"
---

A laser pulse travelling in an optical fiber
causes a nonlinear change of the material's refractive index,
and the resulting dynamics are described by
the [nonlinear Schrödinger (NLS) equation](/know/concept/nonlinear-schrodinger-equation/),
given in its most basic form by:

$$\begin{aligned}
    0
    = i\pdv{A}{z} - \frac{\beta_2}{2} \pdvn{2}{A}{t} + \gamma_0 |A|^2 A
\end{aligned}$$

Where $$A(z, t)$$ is the modulation profile of the carrier wave,
$$\beta_2$$ is the group velocity dispersion
at the carrier frequency $$\omega_0$$,
and $$\gamma_0 \equiv \gamma(\omega_0)$$ is a nonlinear parameter
involving the material's Kerr coefficient $$n_2$$
and the transverse mode's effective area $$A_\mathrm{eff}$$:

$$\begin{aligned}
    \gamma(\omega)
    \equiv \frac{\omega n_2(\omega)}{c A_\mathrm{eff}(\omega)}
\end{aligned}$$

As a consequence of treating $$\gamma_0$$ as frequency-independent,
only the nonlinear *phase* velocity change is represented,
but not the *group* velocity change.
Unfortunately, this form of the NLS equation
does not allow us to include the full $$\gamma(\omega)$$
(this is an advanced topic, see Lægsgaard),
but a decent approximation is to simply Taylor-expand $$\gamma(\omega)$$ around $$\omega_0$$:

$$\begin{aligned}
    \gamma(\omega)
    = \gamma_0 + \gamma_1 \Omega + \frac{\gamma_2}{2} \Omega^2 + \frac{\gamma_3}{6} \Omega^2 + ...
\end{aligned}$$

Where $$\Omega \equiv \omega - \omega_0$$
and $$\gamma_n \equiv \ipdvn{n}{\gamma}{\omega}|_{\omega=\omega_0}$$.
For pulses with a sufficiently narrow spectrum,
we only need the first two terms.
We insert this into the [Fourier transform (FT)](/know/concept/fourier-transform/)
$$\hat{\mathcal{F}}$$ of the equation,
where $$s = \pm 1$$ is the sign of the FT exponent,
which might vary from author to author
($$s = +1$$ corresponds to a forward-propagating carrier wave and vice versa):

$$\begin{aligned}
    0
    = i\pdv{A}{z} - \frac{\beta_2}{2} (-i s \Omega)^2 A + (\gamma_0 + \gamma_1 \Omega) \hat{\mathcal{F}}\big\{ |A|^2 A \big\}
\end{aligned}$$

If we now take the inverse FT,
the factor $$\Omega$$ becomes an operator $$i s \ipdv{}{t}$$:

$$\begin{aligned}
    0
    = i\pdv{A}{z} - \frac{\beta_2}{2} \pdvn{2}{A}{t} + \Big( \gamma_0 + i s \gamma_1 \pdv{}{t} \Big) |A|^2 A
\end{aligned}$$

In theory, this is the desired new NLS equation,
but in fact most authors make a small additional approximation.
Let us write out the derivative of $$\gamma(\omega)$$:

$$\begin{aligned}
    \pdv{\gamma}{\omega}
    = \frac{n_2}{c A_\mathrm{eff}}
    + \frac{\omega}{c A_\mathrm{eff}} \pdv{n_2}{\omega}
    - \frac{\omega n_2}{c A_\mathrm{eff}^2} \pdv{A_\mathrm{eff}}{\omega}
\end{aligned}$$

In practice, the $$\omega$$-dependence of $$n_2$$ and $$A_\mathrm{eff}$$
is relatively weak, so the first term is dominant
and hence sufficient for our purposes.
We therefore have $$\gamma_1 \approx \gamma_0 / \omega_0$$, leading to:

$$\begin{aligned}
    \boxed{
        0
        = i\pdv{A}{z} - \frac{\beta_2}{2} \pdvn{2}{A}{t} + \gamma_0 \Big( 1 + i \frac{s}{\omega_0} \pdv{}{t} \Big) |A|^2 A
    }
\end{aligned}$$

Beware that this NLS equation does not conserve the total energy
$$E \equiv \int_{-\infty}^\infty |A|^2 \dd{t}$$ anymore,
which is often used to quantify simulation errors.
Fortunately, another value can then be used instead:
it can be shown that the "photon number" $$N$$
is still conserved, defined as:

$$\begin{aligned}
    \boxed{
        N(z)
        \equiv \int_{-\infty}^\infty \frac{|\tilde{A}(z, \Omega)|^2}{\Omega} \dd{\Omega}
    }
\end{aligned}$$


A pulse's intensity is highest at its peak,
so the nonlinear index shift is strongest there,
meaning that the peak travels slightly slower than the rest of the pulse,
leading to **self-steepening** of its trailing edge;
an effect exhibited by our modified NLS equation.
Note that $$s$$ controls which edge is regarded as the trailing one.

Let us make the ansatz below,
consisting of an arbitrary power profile $$P$$ with phase $$\phi$$:

$$\begin{aligned}
    A(z,t)
    = \sqrt{P(z,t)} \, \exp\!\big(i \phi(z,t)\big)
\end{aligned}$$

We assume that $$A$$ has a sufficiently narrow spectrum
that we can neglect dispersion $$\beta_2 = 0$$ over a short distance.
Inserting the ansatz into the NLS equation
with $$\varepsilon \equiv \gamma_0 / \omega_0$$ gives:

$$\begin{aligned}
    0
    &= i \frac{1}{2} \frac{P_z}{\sqrt{P}} - \sqrt{P} \phi_z + \gamma_0 P \sqrt{P}
    + i s \varepsilon \frac{3}{2} P_t \sqrt{P} - s \varepsilon P \sqrt{P} \phi_t
\end{aligned}$$

Since $$P$$ is real, this results in two equations,
for the real and imaginary parts:

$$\begin{aligned}
    0
    &= - \phi_z  + \gamma_0 P - s \varepsilon P \phi_t
    \\
    0
    &= P_z + 3 s \varepsilon P_t P
\end{aligned}$$

The phase $$\phi$$ is not so interesting, so we focus on the latter equation for $$P$$.
You can easily show (by insertion) that it has a general solution of the form below,
which says that more intense parts of the pulse
lag behind the rest, as expected:

$$\begin{aligned}
    P(z,t)
    = f(t - 3 s \varepsilon z P)
\end{aligned}$$

Where $$f(t) \equiv P(0,t)$$ is the initial power profile.
The derivatives $$P_t$$ and $$P_z$$ are given by:

$$\begin{aligned}
    P_t
    &= (1 - 3 s \varepsilon z P_t) \: f'
    \qquad\quad\!\! = \frac{f'}{1 + 3 s \varepsilon z f'}
    \\
    P_z
    &= (-3 s \varepsilon P - 3 s \varepsilon z P_z) \: f'
    = \frac{- 3 s \varepsilon P f'}{1 + 3 s \varepsilon z f'}
\end{aligned}$$

Both expressions blow up when their denominator goes to zero,
which, since $$\varepsilon > 0$$, happens earliest at an extremum of $$f'$$;
either its minimum ($$s = +1$$) or maximum ($$s = -1$$).
Let us call this value $$f_\mathrm{extr}'$$,
located on the trailing edge of the pulse.
At the propagation distance $$z$$ where this occurs, $$L_\mathrm{shock}$$,
the pulse "tips over", creating a discontinuous shock:

$$\begin{aligned}
    0
    = 1 + 3 s \varepsilon z f_\mathrm{extr}'
    \qquad \implies \qquad
    z
    = \boxed{
        L_\mathrm{shock}
        \equiv -\frac{\omega_0}{3 s \gamma_0 f_\mathrm{extr}'}
    }
\end{aligned}$$

In practice, however, this never actually happens,
because as the pulse approaches $$L_\mathrm{shock}$$,
its spectrum becomes so broad that dispersion cannot be neglected:
[dispersive broadening](/know/concept/dispersive-broadening/)
pulls the pulse apart before a shock can occur.
The early steepening is observable though.

A simulation of self-steepening without dispersion is illustrated below
for the following initial power distribution,
with $$T_0 = 25\:\mathrm{fs}$$, $$P_0 = 3\:\mathrm{kW}$$,
$$\beta_2 = 0$$, $$\gamma = 0.1/\mathrm{W}/\mathrm{m}$$,
and a vacuum carrier wavelength $$\lambda_0 \approx 73\:\mathrm{nm}$$
(the latter determined by the simulation's resolution settings):

$$\begin{aligned}
    f(t)
    = P(0,t) = P_0 \exp\!\bigg(\!-\!\frac{t^2}{T_0^2} \bigg)
\end{aligned}$$

The first and second derivatives of this Gaussian $$f$$ are as follows:

$$\begin{aligned}
    f'(t)
    &= - \frac{2 P_0}{T_0^2} t \exp\!\bigg(\!-\!\frac{t^2}{T_0^2} \bigg)
    \\
    f''(t)
    &= \frac{2 P_0}{T_0^2} \bigg( \frac{2 t^2}{T_0^2} - 1 \bigg) \exp\!\bigg(\!-\!\frac{t^2}{T_0^2} \bigg)
\end{aligned}$$

The steepest points of $$f'$$ are the roots of $$f''$$,
clearly located at $$2 t^2 = T_0^2$$,
meaning that $$f_\mathrm{extr}'$$ and $$L_\mathrm{shock}$$
are in this case given by:

$$\begin{aligned}
    f_\mathrm{extr}'
    = \mp \sqrt{2} e^{-1/2} \frac{P_0}{T_0}
    \qquad \implies \qquad
    L_\mathrm{shock}
    = \frac{e^{1/2}}{3 \sqrt{2}} \frac{\omega_0 T_0}{\gamma_0 P_0}
\end{aligned}$$

This example Gaussian pulse therefore has a theoretical
$$L_\mathrm{shock} = 0.847\,\mathrm{m}$$,
which seems to be accurate based on these plots,
although the simulation breaks down just before that point due to insufficient resolution:

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

Unfortunately, self-steepening cannot be simulated perfectly: as the
pulse approaches $$L_\mathrm{shock}$$, its spectrum broadens to infinite
frequencies to represent the singularity in its slope.
The simulation thus collapses into chaos when the edge of the frequency window is reached.
Nevertheless, the trend is nicely visible:
the trailing slope becomes extremely steep, and the spectrum
broadens so much that dispersion can no longer be neglected.

{% comment %}
When self-steepening is added to the nonlinear Schrödinger equation,
it no longer conserves the total pulse energy $$\int |A|^2 \dd{t}$$.
Fortunately, the photon number $$N_\mathrm{ph}$$ is still
conserved, which for the physical envelope $$A(z,t)$$ is defined as:

$$\begin{aligned}
    \boxed{
        N_\mathrm{ph}(z) = \int_0^\infty \frac{|\tilde{A}(z,\omega)|^2}{\omega} \dd{\omega}
    }
\end{aligned}$$
{% endcomment %}



## References

1.  B.R. Suydam,
    [Self-steepening of optical pulses](https://doi.org/10.1007/0-387-25097-2_6),
    2006, Springer.
2.  J. Lægsgaard,
    [Mode profile dispersion in the generalized nonlinear Schrödinger equation](https://doi.org/10.1364/OE.15.016110),
    2007, Optica.