Self-Induced Transparency: Area Theorem and Solitons

The McCall–Hahn area theorem governs how the pulse area \(\mathcal{A} = \int_{-\infty}^\infty \Omega(t)\,dt\) evolves as a pulse propagates through a resonant absorbing medium:

\[\frac{d\mathcal{A}}{dz} = -\frac{\alpha}{2}\sin\mathcal{A}\]

This has the same form as the equation of motion for a pendulum, and its fixed points tell us everything about what happens to any pulse entering the medium:

Fixed point	Stability	Outcome
0	stable	weak pulses absorbed
π	unstable	pulses tip toward 0 or 2π
2π	stable	SIT soliton propagates without loss
3π	unstable	one 2π soliton emitted
4π	stable	two 2π solitons
…

Even multiples of π are stable: a pulse with area near \(2n\pi\) reshapes toward exactly \(2n\pi\), splitting into \(n\) solitons that propagate without attenuation. This is self-induced transparency (SIT). Odd multiples are unstable saddle points.

The area theorem is also shape-independent: only the integrated area matters, not the pulse envelope. A Gaussian with area \(1.8\pi\) reshapes toward \(2\pi\) just as a sech pulse does.

This notebook demonstrates the full picture, starting from the fundamental SIT soliton and then exploring what the area theorem predicts for other input areas.

import numpy as np
from maxwellbloch import mb_solve, plot

from maxwellbloch.utility import SECH_FWHM_CONV
t_width = 0.001
t_min = -t_width * 10
t_max = t_width * 40
interaction_strength = 5 / t_width
gauss_fwhm = t_width / SECH_FWHM_CONV
print(f"t_width={t_width}, t_min={t_min:.4f}, t_max={t_max:.4f}, interaction_strength={interaction_strength:.1f}")

t_width=0.001, t_min=-0.0100, t_max=0.0400, interaction_strength=5000.0

The 2π SIT soliton

A sech pulse with area exactly \(2\pi\) is the central object of self-induced transparency. It is an exact analytical solution of the Maxwell-Bloch equations: it propagates through the resonant medium without any attenuation, slowed relative to vacuum by a delay proportional to the optical depth. The medium absorbs energy from the leading edge and re-emits it coherently into the trailing edge — the two processes cancel exactly.

This is the stable fixed point at \(\mathcal{A} = 2\pi\). All other inputs are understood in relation to it.

mbs_2pi = mb_solve.MBSolve().from_json_str(f"""
{{
  "atom": {{
    "fields": [{{"coupled_levels": [[0, 1]], "rabi_freq_t_func": "sech",
      "rabi_freq_t_args": {{"n_pi": 2.0, "centre": 0.0, "width": {t_width}}}}}],
    "num_states": 2,
    "decays": [{{"channels": [[0, 1]], "rate": 1.0}}]
  }},
  "t_min": {t_min}, "t_max": {t_max}, "t_steps": 500,
  "z_min": -0.5, "z_max": 1.5, "z_steps": 200,
  "interaction_strengths": [{interaction_strength}], "savefile": "mbs-sit-2pi"
}}
"""
)
mbs_2pi.mbsolve(recalc=False);

100%|██████████| 200/200 [00:07<00:00, 27.74z/s, Ω_max=1.65e+03]

Saving MBSolve to mbs-sit-2pi.qu

fig = plot.field_spacetime(mbs_2pi, show_z_bounds=[0, 1])
fig.update_layout(title="2π sech — SIT soliton propagates without loss")
fig.show(renderer='notebook_connected')

fig = plot.pulse_area(mbs_2pi)
fig.update_layout(title="2π sech — area constant throughout", yaxis_range=[0, 4])
fig.show(renderer='notebook_connected')

1π — absorbed

A sech pulse with area \(\mathcal{A} < \pi\) lies in the basin of attraction of the zero fixed point: it is simply absorbed. There is no coherent re-emission to balance the absorption, so the pulse decays monotonically as it propagates. The pulse area tracks the area theorem prediction, falling smoothly toward zero.

mbs_1pi = mb_solve.MBSolve().from_json_str(f"""
{{
  "atom": {{
    "fields": [{{"coupled_levels": [[0, 1]], "rabi_freq_t_func": "sech",
      "rabi_freq_t_args": {{"n_pi": 1.0, "centre": 0.0, "width": {t_width}}}}}],
    "num_states": 2,
    "decays": [{{"channels": [[0, 1]], "rate": 1.0}}]
  }},
  "t_min": {t_min}, "t_max": {t_max}, "t_steps": 500,
  "z_min": -0.5, "z_max": 1.5, "z_steps": 200,
  "interaction_strengths": [{interaction_strength}], "savefile": "mbs-sit-1pi"
}}
"""
)
mbs_1pi.mbsolve(recalc=False);

100%|██████████| 200/200 [00:07<00:00, 28.18z/s, Ω_max=116]

Saving MBSolve to mbs-sit-1pi.qu

fig = plot.field_spacetime(mbs_1pi, show_z_bounds=[0, 1])
fig.update_layout(title="1π sech — absorbed")
fig.show(renderer='notebook_connected')

fig = plot.pulse_area(mbs_1pi)
fig.update_layout(title="1π sech — area decays to 0", yaxis_range=[0, 4])
fig.show(renderer='notebook_connected')

3π — one 2π soliton emitted

The \(3\pi\) fixed point is unstable: any perturbation drives the area away from it. Propagation tips a \(3\pi\) input downward toward the nearest stable point at \(2\pi\). One \(2\pi\) soliton forms and the remaining area (\(3\pi - 2\pi = \pi\)) falls below the unstable fixed point at \(\pi\) and is absorbed. The area plot shows the transition clearly.

mbs_3pi = mb_solve.MBSolve().from_json_str(f"""
{{
  "atom": {{
    "fields": [{{"coupled_levels": [[0, 1]], "rabi_freq_t_func": "sech",
      "rabi_freq_t_args": {{"n_pi": 3.0, "centre": 0.0, "width": {t_width}}}}}],
    "num_states": 2,
    "decays": [{{"channels": [[0, 1]], "rate": 1.0}}]
  }},
  "t_min": {t_min}, "t_max": {t_max}, "t_steps": 500,
  "z_min": -0.5, "z_max": 1.5, "z_steps": 200,
  "interaction_strengths": [{interaction_strength}], "savefile": "mbs-sit-3pi"
}}
"""
)
mbs_3pi.mbsolve(recalc=False);

100%|██████████| 200/200 [00:07<00:00, 28.22z/s, Ω_max=3.89e+03]

Saving MBSolve to mbs-sit-3pi.qu

fig = plot.field_spacetime(mbs_3pi, show_z_bounds=[0, 1])
fig.update_layout(title="3π sech — reshapes to one 2π soliton")
fig.show(renderer='notebook_connected')

fig = plot.pulse_area(mbs_3pi)
fig.update_layout(title="3π sech — area settles at 2π", yaxis_range=[0, 4])
fig.show(renderer='notebook_connected')

Shape independence — Gaussian pulses

The area theorem contains no information about pulse shape — only \(\mathcal{A}\) appears in \(d\mathcal{A}/dz = -(\alpha/2)\sin\mathcal{A}\). A Gaussian, square, or any other envelope with the same area undergoes the same fate. A \(1.8\pi\) Gaussian reshapes toward the \(2\pi\) sech soliton, just as a \(1.8\pi\) sech pulse would — the attractor depends only on the area, not the shape.

mbs_g18 = mb_solve.MBSolve().from_json_str(f"""
{{
  "atom": {{
    "fields": [{{"coupled_levels": [[0, 1]], "rabi_freq_t_func": "gaussian",
      "rabi_freq_t_args": {{"n_pi": 1.8, "centre": 0.0, "fwhm": {gauss_fwhm}}}}}],
    "num_states": 2,
    "decays": [{{"channels": [[0, 1]], "rate": 1.0}}]
  }},
  "t_min": {t_min}, "t_max": {t_max}, "t_steps": 500,
  "z_min": -0.5, "z_max": 1.5, "z_steps": 200,
  "interaction_strengths": [{interaction_strength}], "savefile": "mbs-sit-gauss-18pi"
}}
"""
)
mbs_g18.mbsolve(recalc=False)
print(f"Input area: {np.trapezoid(mbs_g18.Omegas_zt[0, 0, :].real, mbs_g18.tlist) / np.pi:.3f}π")

100%|██████████| 200/200 [00:07<00:00, 28.24z/s, Ω_max=1.64e+03]

Saving MBSolve to mbs-sit-gauss-18pi.qu
Input area: 1.800π

fig = plot.field_spacetime(mbs_g18, show_z_bounds=[0, 1])
fig.update_layout(title="1.8π Gaussian — reshapes toward 2π soliton")
fig.show(renderer='notebook_connected')

fig = plot.pulse_area(mbs_g18)
fig.update_layout(title="1.8π Gaussian — area asymptotes to 2π", yaxis_range=[0, 3])
fig.show(renderer='notebook_connected')

The \(1.8\pi\) Gaussian reshapes toward \(2\pi\) just as a sech pulse does. The attractor is always the same \(2\pi\) sech soliton regardless of input shape — the pulse distorts in transit, shedding the excess as radiation. Having established the area theorem for single-soliton outcomes, we now look at inputs large enough to produce multiple solitons.

4π — two solitons

A \(4\pi\) input is near a second stable fixed point. It breaks into two \(2\pi\) solitons, each carrying half the original area. The two solitons emerge with different widths and therefore different group velocities — the narrower soliton travels faster — so they separate as they propagate through the medium.

mbs_4pi = mb_solve.MBSolve().from_json_str(f"""
{{
  "atom": {{
    "fields": [{{"coupled_levels": [[0, 1]], "rabi_freq_t_func": "sech",
      "rabi_freq_t_args": {{"n_pi": 4.0, "centre": 0.0, "width": {t_width}}}}}],
    "num_states": 2,
    "decays": [{{"channels": [[0, 1]], "rate": 1.0}}]
  }},
  "t_min": {t_min}, "t_max": {t_max}, "t_steps": 500,
  "z_min": -0.5, "z_max": 1.5, "z_steps": 200,
  "interaction_strengths": [{interaction_strength}], "savefile": "mbs-sit-4pi"
}}
"""
)
mbs_4pi.mbsolve(recalc=False);

100%|██████████| 200/200 [00:07<00:00, 26.96z/s, Ω_max=5.83e+03]

Saving MBSolve to mbs-sit-4pi.qu

fig = plot.field_spacetime(mbs_4pi, show_z_bounds=[0, 1])
fig.update_layout(title="4π sech — breaks into two 2π solitons")
fig.show(renderer='notebook_connected')

fig = plot.pulse_area(mbs_4pi)
fig.update_layout(title="4π sech — total area conserved", yaxis_range=[0, 6])
fig.show(renderer='notebook_connected')

6π — three solitons

A \(6\pi\) input breaks into three \(2\pi\) solitons. The general rule is clear: a \(2n\pi\) input produces exactly \(n\) solitons. The solitons again emerge with different widths and speeds, fanning out as they propagate.

mbs_6pi = mb_solve.MBSolve().from_json_str(f"""
{{
  "atom": {{
    "fields": [{{"coupled_levels": [[0, 1]], "rabi_freq_t_func": "sech",
      "rabi_freq_t_args": {{"n_pi": 6.0, "centre": 0.0, "width": {t_width}}}}}],
    "num_states": 2,
    "decays": [{{"channels": [[0, 1]], "rate": 1.0}}]
  }},
  "t_min": {t_min}, "t_max": {t_max}, "t_steps": 500,
  "z_min": -0.5, "z_max": 1.5, "z_steps": 200, "z_steps_inner": 2,
  "interaction_strengths": [{interaction_strength}], "savefile": "mbs-sit-6pi"
}}
"""
)
mbs_6pi.mbsolve(recalc=False);

100%|██████████| 200/200 [00:07<00:00, 26.01z/s, Ω_max=9.54e+03]

Saving MBSolve to mbs-sit-6pi.qu

fig = plot.field_spacetime(mbs_6pi, show_z_bounds=[0, 1])
fig.update_layout(title="6π sech — breaks into three 2π solitons")
fig.show(renderer='notebook_connected')

fig = plot.pulse_area(mbs_6pi)
fig.update_layout(title="6π sech — total area conserved", yaxis_range=[0, 8])
fig.show(renderer='notebook_connected')

Soliton collision

The Maxwell-Bloch equations in the sharp-line limit are integrable, which means SIT solitons collide elastically: two solitons passing through each other emerge with their original shapes and speeds intact, shifted only in position. A wide (slow) and a narrow (fast) \(2\pi\) soliton are launched together so the faster one overtakes the slower inside the medium.

t_width_wide   = 2.0 * t_width
t_width_narrow = 1.0 * t_width
t_max_col = t_width * 80

mbs_col = mb_solve.MBSolve().from_json_str(f"""
{{
  "atom": {{"fields": [{{"coupled_levels": [[0, 1]]}}], "num_states": 2, "decays": [{{"channels": [[0, 1]], "rate": 1.0}}]}},
  "t_min": {t_min}, "t_max": {t_max_col}, "t_steps": 500,
  "z_min": -0.5, "z_max": 1.5, "z_steps": 200,
  "interaction_strengths": [{interaction_strength}], "savefile": "mbs-sit-collision"
}}
"""
)

from maxwellbloch import t_funcs

probe = mbs_col.atom.fields[0]
probe.rabi_freq_t_func = lambda t, args: t_funcs.sech(1)(t, args) + t_funcs.sech(2)(t, args)
probe.rabi_freq_t_args = {
    "n_pi_1": 2.0, "centre_1": 0.0, "width_1": t_width_wide,
    "n_pi_2": 2.0, "centre_2": t_width * 5, "width_2": t_width_narrow,
}
mbs_col.atom.build_operators()
mbs_col.init_Omegas_zt()
mbs_col.mbsolve(recalc=False);

100%|██████████| 200/200 [00:07<00:00, 28.11z/s, Ω_max=2.37e+03]

Saving MBSolve to mbs-sit-collision.qu

fig = plot.field_spacetime(mbs_col, show_z_bounds=[0, 1])
fig.update_layout(title="2π soliton collision — elastic scattering")
fig.show(renderer='notebook_connected')

The two solitons briefly merge into a high-amplitude transient, then re-emerge unchanged in width and speed. The only lasting effect is a position shift — each soliton arrives slightly earlier or later than it would have without the collision. This is the hallmark of elastic, soliton-preserving scattering.

Summary

Input	Outcome
2π sech	SIT soliton — propagates without loss
1π sech	absorbed
3π sech	one 2π soliton emitted, remainder absorbed
1.8π Gaussian	reshapes toward 2π (shape-independent)
4π sech	two 2π solitons
6π sech	three 2π solitons

The area theorem unifies all of these: \(d\mathcal{A}/dz = -(\alpha/2)\sin\mathcal{A}\) has a pendulum-like fixed-point structure, and the SIT soliton is the \(2\pi\) stable fixed point realised as an exact propagating solution. Soliton collisions are elastic because the underlying equations are integrable.

References

1. 19. 12. McCall, E. L. Hahn, Self-Induced Transparency by Pulsed Coherent Light, PRL 18, 908 (1967).

2. 19. 12. McCall, E. L. Hahn, Self-Induced Transparency, Phys. Rev. 183, 457 (1969).

3. 12. Allen, J. H. Eberly, Optical Resonance and Two-Level Atoms (Dover, 1987), Ch. 4.