The Optical Bloch Equations

Introduction

The Optical Bloch Equations (OBEs) describe the dynamics of a two-level quantum system interacting with a coherent electromagnetic field while subject to spontaneous emission and decoherence.

Whereas the wavefunction formalism captures only coherent evolution, the density-matrix approach naturally incorporates population decay, loss of coherence, and steady-state behavior. For this reason, OBEs form the foundational model for understanding how atoms respond to resonant and near-resonant light.

In the previous chapter, we introduced the density matrix of a two-level atom and distinguished pure, mixed, and partially coherent states. We now build on that framework by including:

• coherent driving by a monochromatic laser field,
• spontaneous decay of the excited state, and
• the resulting relaxation of populations and coherences.

From the Hamiltonian in the rotating frame and the Liouville–von Neumann equation, augmented with phenomenological decay terms, we will derive a closed set of coupled equations for the density-matrix elements. These equations—collectively the Optical Bloch Equations—provide the basis for understanding Rabi oscillations, saturation, optical pumping, and the steady-state and transient response of atoms to laser light.

By the end of this chapter, you should be able to:

• Incorporate spontaneous decay and decoherence into the density-matrix evolution.
  
• Derive the Optical Bloch Equations from the Liouville–von Neumann equation.
  
• Interpret the physical meaning of each term in the OBEs (population decay, coherence damping, driving terms).
  
• Solve the OBEs in steady state and identify regimes of saturation and power broadening.
  
• Describe the behavior of a two-level atom under resonant and off-resonant driving, including transient Rabi oscillations and long-time limits.

In the following section, we derive the Optical Bloch Equations step by step and explore their physical content.

---

1. The Optical Bloch Equations (OBEs)

The time evolution of the density matrix of a two-level atom interacting with a classical laser field is obtained by inserting the Hamiltonian from Equation 5 on the page about the Rotating Wave Approximation into the Liouville–von Neumann equation.

In this section we derive the Optical Bloch Equations (OBEs). The derivation proceeds in two main parts:

1. We first supplement the unitary evolution with phenomenological relaxation terms describing spontaneous emission.
   
2. We then insert the Hamiltonian in the laser-rotating frame and explicitly evaluate the commutator, obtaining the full set of OBEs.

---

1.1 Spontaneous decay

Real atoms decay from the excited state even in the absence of external fields. This effect is not contained in the Hamiltonian, because spontaneous emission originates from vacuum fluctuations of the quantized electromagnetic field. Since we are treating the field classically, we incorporate spontaneous decay phenomenologically by modifying the Liouville–von Neumann equation.

Exponential decay of the excited state

Without any driving field, the excited-state population decays exponentially with rate \(\gamma = 1/\tau\):

\[ \rho_{ee}(t) = \rho_{ee}(0)\, e^{-\gamma t}. \]

Taking the derivative gives the spontaneous-decay contribution to the diagonal elements:

\[ \left(\frac{\partial \rho_{ee}}{\partial t}\right)_{\mathrm{spont.dec.}} = -\left(\frac{\partial \rho_{gg}}{\partial t}\right)_{\mathrm{spont.dec.}} = -\gamma \rho_{ee}. \]

Relaxation of the coherences

The off-diagonal elements represent quantum coherence between the states, and coherence is lost when population decays. A more rigorous derivation requires quantized-field perturbation theory, but the well-established result is:

\[ \left(\frac{\partial \rho_{eg}}{\partial t}\right)_{\mathrm{spont.dec.}} = -\frac{\gamma}{2}\rho_{eg}, \qquad \left(\frac{\partial \rho_{ge}}{\partial t}\right)_{\mathrm{spont.dec.}} = -\frac{\gamma}{2}\rho_{ge}. \]

Modified Liouville–von Neumann equation

Including these relaxation terms yields:

\[ \frac{\partial \hat\rho}{\partial t} = \frac{1}{i\hbar}[\hat H,\hat\rho] + \begin{pmatrix} -\gamma\rho_{ee} & -\tfrac{\gamma}{2}\rho_{eg} \\ -\tfrac{\gamma}{2}\rho_{ge} & \gamma\rho_{ee} \end{pmatrix}. \tag{1}\]

---

1.2 Coherent evolution in the rotating frame

To obtain the full Optical Bloch Equations, we now include the coherent atom–field interaction.
We use the Hamiltonian in the laser-rotating frame derived in the RWA section:

• It removes the fast optical-frequency oscillations,
  
• It introduces the detuning \(\delta = \omega_L - \omega_0\),
  
• It expresses the field coupling via the complex Rabi frequency \(\Omega_0\).

Before adding relaxation, we first compute the unitary part of the evolution: the commutator \([\hat H,\hat\rho]\) in the rotating frame.

Evaluating the commutator

Using the Hamiltonian from the Rotating-Wave Approximation, we obtain

\[ \begin{align} \frac{2}{\hbar}[\hat H,\hat\rho] =& \begin{pmatrix} -\delta & -\Omega_0\\ -\Omega_0^* & \delta \end{pmatrix} \begin{pmatrix} \rho_{ee} & \rho_{eg} \\ \rho_{ge} & \rho_{gg} \end{pmatrix} - \begin{pmatrix} \rho_{ee} & \rho_{eg} \\ \rho_{ge} & \rho_{gg} \end{pmatrix} \begin{pmatrix} -\delta & -\Omega_0\\ -\Omega_0^* & \delta \end{pmatrix} \\ =& \begin{pmatrix} \Omega_0^\ast \rho_{eg} - \Omega_0 \rho_{ge} & \Omega_0(\rho_{ee}-\rho_{gg}) - 2\delta\rho_{eg} \\ -\Omega_0^\ast(\rho_{ee}-\rho_{gg}) + 2\delta\rho_{ge} & \Omega_0 \rho_{ge} - \Omega_0^\ast \rho_{eg} \end{pmatrix}. \end{align} \]

Full evolution: unitary + relaxation

Substituting this into Equation 1 gives the matrix form of the OBEs:

\[ \frac{\partial}{\partial t} \begin{pmatrix} \rho_{ee} & \rho_{eg} \\ \rho_{ge} & \rho_{gg} \end{pmatrix} = \begin{pmatrix} -\tfrac{i}{2}(\Omega_0^\ast\rho_{eg}-\Omega_0\rho_{ge}) - \gamma\rho_{ee} & -\tfrac{i\Omega_0}{2}(\rho_{ee}-\rho_{gg}) + i\delta\rho_{eg} - \tfrac{\gamma}{2}\rho_{eg} \\ \tfrac{i\Omega_0^\ast}{2}(\rho_{ee}-\rho_{gg}) - i\delta\rho_{ge} - \tfrac{\gamma}{2}\rho_{ge} & \tfrac{i}{2}(\Omega_0^\ast\rho_{ge}-\Omega_0\rho_{ge}) + \gamma\rho_{ee} \end{pmatrix}. \tag{2}\]

Component form: the Optical Bloch Equations

Writing Equation 2 element-wise:

\[ \frac{\partial\rho_{ee}}{\partial t} = -\frac{i}{2}(\Omega_0^\ast\rho_{eg} - \Omega_0\rho_{ge}) - \gamma\rho_{ee}, \tag{3}\]

\[ \frac{\partial\rho_{gg}}{\partial t} = \frac{i}{2}(\Omega_0^\ast\rho_{eg} - \Omega_0\rho_{ge}) + \gamma\rho_{ee}, \tag{4}\]

\[ \frac{\partial\rho_{eg}}{\partial t} = -\frac{i\Omega_0}{2}(\rho_{ee}-\rho_{gg}) + i\delta\rho_{eg} - \frac{\gamma}{2}\rho_{eg}, \tag{5}\]

\[ \frac{\partial\rho_{ge}}{\partial t} = \frac{i\Omega_0^\ast}{2}(\rho_{ee}-\rho_{gg}) - i\delta\rho_{ge} - \frac{\gamma}{2}\rho_{ge}. \tag{6}\]

These are the Optical Bloch Equations.

---

1.3 Reducing the system: inversion and coherence

Because the density matrix satisfies

\[ \rho_{ee} + \rho_{gg} = 1, \qquad \rho_{eg} = \rho_{ge}^\ast, \]

the OBEs contain redundant information and can be reduced to two coupled differential equations.

Defining the inversion

Define the inversion:

\[ w = \rho_{ee} - \rho_{gg}. \tag{7}\]

Using \(\rho_{gg} = 1 - \rho_{ee}\), we obtain

\[ w = \rho_{ee}-\rho_{gg} = \rho_{ee}-(1-\rho_{ee}) = 2\rho_{ee}-1 \quad\Rightarrow\quad \rho_{ee} = \tfrac12(w+1). \tag{8}\]

Reduced-form OBEs

Then the OBEs reduce to the pair

\[ \frac{\partial\rho_{eg}}{\partial t} = -\frac{i\Omega_0}{2}\, w + \left(i\delta - \frac{\gamma}{2}\right)\rho_{eg}, \tag{9}\]

\[ \frac{\partial w}{\partial t} = -i(\Omega_0^\ast\rho_{eg} - \Omega_0\rho_{ge}) - \gamma(w+1). \tag{10}\]

Equation Equation 9 is the same as Equation 5 (and the complex conjugate of Equation 6). Equation Equation 10 follows from subtracting Equation 4 from Equation 3 and substituting Equation 8.

These two equations fully describe the dynamics of a driven two-level atom in the RWA including spontaneous decay.

---

2. Solutions of the Optical Bloch Equations

The Optical Bloch Equations (OBEs) describe the driven dynamics of a damped two-level system in the rotating frame. Analytic solutions exist in important special cases (e.g., steady state, resonant drive, weak driving), but in the general case the equations are nonlinear and must be solved numerically. Before turning to analytic limits, we illustrate the full dynamics via direct numerical integration.

2.1 General solution (numerical)

The rotating-frame OBEs for a two-level atom driven with Rabi frequency \(\Omega_0\) and detuning \(\delta\) were derived in Equation 3–Equation 6.

These coupled first-order equations are well suited to numerical integration. Below we use a fourth-order Runge–Kutta scheme to compute the population dynamics for several choices of \(\gamma\).

The code is adapted from the MATLAB implementation of Sathyanarayan Rao (2015).

"""
Damped Rabi Oscillations – Rotating Frame (Python)
Adapted from the MATLAB code by Sathyanarayan Rao (2015).
"""

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['text.usetex'] = True
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.size'] = 16

# --- Optical Bloch equations in rotating frame ---
def drho11(t, r12, r21, r22, gam, Omega, delta):
    return 1j*0.5*np.conj(Omega)*r12 - 1j*0.5*Omega*r21 + gam*r22

def drho12(t, r11, r12, r22, gam, Omega, delta):
    return 1j*0.5*Omega*(r11 - r22) - 1j*delta*r12 - 0.5*gam*r12

def drho21(t, r11, r21, r22, gam, Omega, delta):
    return -1j*0.5*np.conj(Omega)*(r11 - r22) + 1j*delta*r21 - 0.5*gam*r21

def drho22(t, r12, r21, r22, gam, Omega, delta):
    return -1j*0.5*np.conj(Omega)*r12 + 1j*0.5*Omega*r21 - gam*r22

# Parameters
Omega = 2*np.pi*10e6        # Rabi frequency 2π×10 MHz
delta = 0                   # detuning
gamma = Omega               # set a typical scale

N = 1000
t = np.linspace(0, 1e-6, N)
dt = t[1] - t[0]

def density_me(Omega, delta, gamma):
    rho11 = np.ones(N, dtype=complex)
    rho22 = np.zeros(N, dtype=complex)
    rho12 = np.zeros(N, dtype=complex)
    rho21 = np.zeros(N, dtype=complex)

    for i in range(1, N):
        k1 = drho22(t[i-1], rho12[i-1], rho21[i-1], rho22[i-1], gamma, Omega, delta)*dt
        l1 = drho11(t[i-1], rho12[i-1], rho21[i-1], rho22[i-1], gamma, Omega, delta)*dt
        m1 = drho12(t[i-1], rho11[i-1], rho12[i-1], rho22[i-1], gamma, Omega, delta)*dt
        n1 = drho21(t[i-1], rho11[i-1], rho21[i-1], rho22[i-1], gamma, Omega, delta)*dt

        k2 = drho22(t[i-1]+dt/2, rho12[i-1]+m1/2, rho21[i-1]+n1/2,
                    rho22[i-1]+k1/2, gamma, Omega, delta)*dt
        l2 = drho11(t[i-1]+dt/2, rho12[i-1]+m1/2, rho21[i-1]+n1/2,
                    rho22[i-1]+k1/2, gamma, Omega, delta)*dt
        m2 = drho12(t[i-1]+dt/2, rho11[i-1]+l1/2, rho12[i-1]+m1/2,
                    rho22[i-1]+k1/2, gamma, Omega, delta)*dt
        n2 = drho21(t[i-1]+dt/2, rho11[i-1]+l1/2, rho21[i-1]+n1/2,
                    rho22[i-1]+k1/2, gamma, Omega, delta)*dt

        k3 = drho22(t[i-1]+dt/2, rho12[i-1]+m2/2, rho21[i-1]+n2/2,
                    rho22[i-1]+k2/2, gamma, Omega, delta)*dt
        l3 = drho11(t[i-1]+dt/2, rho12[i-1]+m2/2, rho21[i-1]+n2/2,
                    rho22[i-1]+k2/2, gamma, Omega, delta)*dt
        m3 = drho12(t[i-1]+dt/2, rho11[i-1]+l2/2, rho12[i-1]+m2/2,
                    rho22[i-1]+k2/2, gamma, Omega, delta)*dt
        n3 = drho21(t[i-1]+dt/2, rho11[i-1]+l2/2, rho21[i-1]+n2/2,
                    rho22[i-1]+k2/2, gamma, Omega, delta)*dt

        k4 = drho22(t[i-1]+dt, rho12[i-1]+m3, rho21[i-1]+n3,
                    rho22[i-1]+k3, gamma, Omega, delta)*dt
        l4 = drho11(t[i-1]+dt, rho12[i-1]+m3, rho21[i-1]+n3,
                    rho22[i-1]+k3, gamma, Omega, delta)*dt
        m4 = drho12(t[i-1]+dt, rho11[i-1]+l3, rho12[i-1]+m3,
                    rho22[i-1]+k3, gamma, Omega, delta)*dt
        n4 = drho21(t[i-1]+dt, rho11[i-1]+l3, rho21[i-1]+n3,
                    rho22[i-1]+k3, gamma, Omega, delta)*dt

        rho22[i] = rho22[i-1] + (k1 + 2*(k2+k3) + k4)/6
        rho11[i] = rho11[i-1] + (l1 + 2*(l2+l3) + l4)/6
        rho12[i] = rho12[i-1] + (m1 + 2*(m2+m3) + m4)/6
        rho21[i] = rho21[i-1] + (n1 + 2*(n2+n3) + n4)/6

    return rho11, rho12, rho21, rho22

_,_,_, rho22_1 = density_me(Omega, delta, 0.1*gamma)
_,_,_, rho22_2 = density_me(Omega, delta, 0.25*gamma)
_,_,_, rho22_3 = density_me(Omega, delta, gamma)

plt.figure(figsize=(7.5,5))
plt.plot(t*Omega/(2*np.pi), np.real(rho22_1), linewidth=2, label=r'$\gamma=0.1\,\Omega_0$')
plt.plot(t*Omega/(2*np.pi), np.real(rho22_2), 'r--', linewidth=2, label=r'$\gamma=0.25\,\Omega_0$')
plt.plot(t*Omega/(2*np.pi), np.real(rho22_3), 'g:', linewidth=2, label=r'$\gamma=\Omega_0$')

plt.axhline(0.5, color='gray', linestyle=':', linewidth=1)
plt.xlabel(r'$t\,\Omega_0/(2\pi)$')
plt.ylabel(r'Excited state population $\rho_{ee}$')
plt.xlim([0, t[-1]*Omega/(2*np.pi)])
plt.ylim([0,1])
plt.legend(fontsize=12)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

Figure 1: Calculated excited-state population vs. time from numerical integration of the OBEs. The laser is resonant (\(\delta =0\)). Curves show \(\gamma=\Omega_0 /10\) (solid), \(\Omega_0 /4\) (dashed), and \(\gamma =\Omega_0\) (dotted).

---

Interpretation of the numerical solution

As shown in Figure 1, spontaneous emission damps the Rabi oscillations. Two characteristic regimes appear:

• Short times \(t \ll 1/\gamma\):
  Dissipation has little effect; the dynamics closely follow coherent Rabi oscillations if \(\Omega_0 > \gamma\).

• Long times \(t \gg 1/\gamma\):
  The system approaches a steady state with no further population transfer. This regime is dominated by the balance between optical pumping and decay.

These behaviors motivate the definition of experimentally meaningful parameters that organize the regimes of operation.

---

2.2 Saturation parameters

Using the definitions of \(\Omega_0\) and \(\gamma\)

\[ \Omega_0=\frac{d_{eg}E_0}{\hbar},\qquad \gamma=A_{eg}=\frac{d_{eg}^2\omega_0^3}{3\pi\varepsilon_0\hbar c^3}, \]

one introduces:

\[ I_{sat} = \frac{\hbar\omega_0^3\,\gamma}{12\pi c^2} = \frac{\pi h c}{3\lambda^3 \tau}, \]

\[ S_0 = \frac{I}{I_{sat}} = \frac{2|\Omega_0|^2}{\gamma^2}, \]

\[ S = \frac{2|\Omega_0|^2}{\gamma^2 + 4\delta^2} = \frac{S_0}{1 + 4\delta^2/\gamma^2}. \tag{11}\]

Physical meaning

• Saturation intensity \(I_{sat}\):
  A property of the transition. For \(I<I_{sat}\) the atom cannot undergo many Rabi cycles before decaying. Typical values for alkali atoms are \(1\)–\(10\) mW/cm².

• Resonant saturation parameter \(S_0\):
  Controlled experimentally via the laser intensity: \(S_0=1\) corresponds to \(I=I_{sat}\).

• General saturation parameter \(S\) (with detuning):
  Reduces to \(S_0\) at \(\delta=0\). Significant Rabi oscillations require \(S\gtrsim 1\).

These will appear naturally in the analytic steady-state solutions of the OBEs in the next section.

---

2.3 Steady-State Solutions

In many practically relevant situations, atoms interact with a near-resonant, monochromatic laser field for times much longer than the natural lifetime of the excited state. In such cases, all transient dynamics have decayed away, and the system reaches a true steady state, meaning that all time derivatives of the density-matrix elements vanish.

Mathematically, the steady state of the Optical Bloch Equations corresponds to imposing

\[ \frac{\partial w}{\partial t} = 0, \qquad \frac{\partial \rho_{eg}}{\partial t} = 0, \]

where \(w\) is the inversion (Eq. Equation 7). We now solve Eqs. Equation 9 and Equation 10 under this condition.

---

Steady-state solutions

Solving the OBEs for \(\dot w=\dot\rho_{eg}=0\) yields the steady-state inversion and coherence:

\[ w_{\mathrm{ss}} = -\,\frac{1}{\,1 + \dfrac{2|\Omega_0|^2}{\gamma^2 + 4\delta^2}\,} = -\frac{1}{1+S}, \] where \[ S \equiv \frac{2|\Omega_0|^2}{\gamma^2 + 4\delta^2} \] is the saturation parameter.

The steady-state coherence follows from Eq. Equation 9:

\[ \rho_{eg}^{(\mathrm{ss})} = \frac{i\Omega_0}{i2\delta - \gamma}\, w_{\mathrm{ss}} = \frac{i\Omega_0}{(\gamma - i2\delta)(1+S)}. \]

These expressions describe the long-time limit of a driven two-level system including spontaneous emission.

---

NoteDerivation

The expression for the steady-state coherence follows directly by solving Eq. Equation 9 for \(\rho_{eg}\) under the condition \(\partial_t\rho_{eg}=0\):

\[ \rho_{eg} = \frac{i\Omega_0}{i2\delta - \gamma}\, w. \]

Because \(\rho_{ge} = \rho_{eg}^\ast\) and the quantities \(\delta\), \(\gamma\), and \(w\) are real, we obtain

\[ \rho_{ge} = \frac{i\Omega_0^\ast}{i2\delta + \gamma}\, w. \]

We now insert \(\rho_{eg}\) and \(\rho_{ge}\) into the right-hand side of Eq. Equation 10 (the equation for \(\dot w\)). The first term becomes

\[ -i\left(\Omega_0^\ast\rho_{eg} - \Omega_0\rho_{ge}\right) = |\Omega_0|^2 w \left( \frac{1}{i2\delta - \gamma} - \frac{1}{i2\delta + \gamma} \right). \]

Evaluating the bracket:

\[ \frac{1}{i2\delta - \gamma} - \frac{1}{i2\delta + \gamma} = -\,\frac{2\gamma}{\gamma^2 + 4\delta^2}. \]

Thus,

\[ -i(\Omega_0^\ast\rho_{eg}-\Omega_0\rho_{ge}) = -\,|\Omega_0|^2 w\, \frac{2\gamma}{\gamma^2 + 4\delta^2}. \]

Including the spontaneous-decay term from Eq. Equation 10 gives

\[ -\,|\Omega_0|^2 w \frac{2\gamma}{\gamma^2+4\delta^2} - \gamma(w+1). \]

In steady state, the left-hand side of Eq. Equation 10 vanishes, so we set the expression above equal to zero and solve for \(w\). The result is

\[ w = -\frac{1}{1 + \dfrac{2|\Omega_0|^2}{\gamma^2 + 4\delta^2}} = -\frac{1}{1+S}. \]

This completes the derivation.

---

Excited-state steady-state population

Using the identity from Eq. Equation 8,

\[ \rho_{ee} = \tfrac12(w+1), \]

we obtain

\[ \rho_{ee}^{(\mathrm{ss})} = \frac{1}{2}\bigl(w_{\mathrm{ss}} + 1\bigr) = \frac{S}{2(1+S)}. \tag{12}\]

A remarkable and important result is that

the excited-state population never exceeds 50 %,

regardless of how intense the driving field becomes.

This saturation behavior is fundamental to two-level systems with spontaneous decay.

---

2.4 Photon Scattering Rate and Intensity Broadening

When an atom is illuminated by near-resonant laser light, it undergoes a continuous cycle of stimulated absorption, stimulated emission, and—at random times—spontaneous decay. The spontaneously emitted photons are the ones detected in a fluorescence experiment (see Fig. Figure 2). A natural question is:

At what rate does the atom emit fluorescence photons?

This rate is called the photon scattering rate \(\Gamma\). Physically, it is simply the spontaneous decay rate \(\gamma\) multiplied by the steady-state probability of finding the atom in the excited state.

Using the steady-state excited-state population \(\rho_{ee}^{(\mathrm{ss})}\) from Equation 12, and the saturation parameter \(S\) (Equation 11), the scattering rate becomes

\[ \Gamma = \gamma \rho_{ee} = \frac{\gamma}{2}\,\frac{S}{1+S} = \frac{\gamma}{2}\, \frac{S_0}{\,1 + S_0 + \dfrac{4\delta^2}{\gamma^2}\,}. \tag{13}\]

---

Figure 2: Sketch of an atom in a laser field that continuously undergoes absorption, stimulated emission, and spontaneous emission. Spontaneous decay (i.e., fluorescence) occurs at the scattering rate \(\Gamma\).

---

Basic properties of the scattering rate

Several important features follow directly from Eq. Equation 13:

• Saturation limit.
  No matter how large the laser intensity becomes,
  \[ \Gamma_{\max} = \frac{\gamma}{2}, \] because the excited-state population never exceeds 1/2 in steady state.

• Linear regime.
  For low intensity (\(S\ll1\)),
  \[ \Gamma \approx \frac{\gamma}{2} S, \] which is proportional to the laser intensity.

• Frequency dependence.
  The detuning \(\delta = \omega_L - \omega_0\) produces a Lorentzian response centered at resonance.

These features form the basis of many fluorescence-based measurement techniques, including Doppler cooling.

---

Lorentzian form and power broadening

In a typical fluorescence or absorption experiment, one scans the laser detuning
\[\delta = \omega_L - \omega_0\]
while keeping the laser intensity fixed. Thus \(\delta\) is the independent variable, and \(S_0\) (the on-resonance saturation parameter) is a constant during the scan.

Starting from Eq. Equation 13, the steady-state scattering rate can be written:

\[ \Gamma = \frac{\gamma}{2} \frac{S_0}{1+S_0 + \dfrac{4\delta^2}{\gamma^2}}. \]

To make the line shape explicit, factor out \((1+S_0)\) from the denominator:

\[ 1 + S_0 + \frac{4\delta^2}{\gamma^2} = (1+S_0)\left[1 + \frac{4\delta^2}{\gamma^2(1+S_0)}\right]. \]

Define the power-broadened linewidth

\[ \gamma' = \gamma\sqrt{\,1+S_0\,}. \]

Then the scattering rate becomes

\[ \Gamma = \frac{\gamma}{2} \frac{S_0}{1+S_0} \frac{\gamma'^2}{\gamma'^2 + (2\delta)^2}. \tag{14}\]

This has the standard Lorentzian form in the variable \(\delta\) (see Figure 3).

---

Physical meaning: why the line broadens

• At low intensity, \(S_0\ll1\), the width is essentially the natural linewidth \(\gamma\).
  
• As intensity increases, \(\gamma'\) grows because the transition is driven more strongly.

This broadening occurs because:

• At resonance, the scattering rate is already close to its saturation value \(\gamma/2\), so increasing intensity cannot increase it much further.
  
• Off resonance, the scattering rate is still far below saturation, so higher intensity gives a noticeable increase.

The net effect is a wider Lorentzian, known as intensity broadening, saturation broadening, or power broadening. This phenomenon plays a central role in laser cooling, optical pumping, and high-resolution spectroscopy.

import numpy as np
import matplotlib.pyplot as plt

# --- Set up plot styling ---
plt.rcParams['text.usetex'] = True
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.size'] = 14


# Define the scattering rate function
def scattering_rate(delta, S0, gamma=1):
    """
    Calculate the photon scattering rate Γ as a function of detuning.
    
    Parameters:
    -----------
    delta : array-like
        Detuning δ = ω_L - ω_0 (in units of γ)
    S0 : float
        On-resonance saturation parameter
    gamma : float
        Natural linewidth (default=1 for normalized units)
    
    Returns:
    --------
    Gamma : array-like
        Scattering rate (in units of γ)
    """
    return (gamma/2) * S0 / (1 + S0 + 4*delta**2/gamma**2)

# Set up the detuning range (in units of γ)
delta = np.linspace(-10, 10, 1000)

# Saturation parameters to plot
S0_values = [0.5, 1, 10, 100]

# Create the figure
fig, ax = plt.subplots(figsize=(7.5, 5))

# Plot scattering rate for each S0
for S0 in S0_values:
    Gamma = scattering_rate(delta, S0)
    ax.plot(delta, Gamma, linewidth=2, label=f'$S_0 = {S0}$')

# Add horizontal line at Γ_max = γ/2
ax.axhline(y=0.5, color='gray', linestyle='--', linewidth=1, 
           label=r'$\Gamma_{\mathrm{max}} = \gamma/2$')

# Formatting
ax.set_xlabel(r'Detuning $\delta/\gamma$', fontsize=14)
ax.set_ylabel(r'Scattering Rate $\Gamma/\gamma$', fontsize=14)
ax.legend(fontsize=12, loc='upper right')
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 0.55)
ax.set_xlim(-10, 10)


plt.tight_layout()
plt.show()

Figure 3: Scattering rate \(\Gamma/\gamma\) as a function of normalized detuning \(\delta/\gamma\) for several fixed values of the on-resonance saturation parameter \(S_0\).

---

Key Takeaways

• The optical Bloch equations (OBEs) describe the full time evolution of a two-level atom interacting with a coherent laser field, incorporating both coherent dynamics and relaxation processes.

• By introducing the density matrix formalism, the OBEs account for population dynamics, atomic coherences, and decoherence due to spontaneous emission.

• The Rabi frequency \(\Omega_0\), detuning \(\delta\), and decay rate \(\gamma\) are the key parameters determining the system’s behavior, from Rabi oscillations to steady-state saturation.

• The OBEs serve as the foundation for understanding a wide range of light–matter phenomena, including optical pumping, coherent population trapping, line broadening, and laser cooling.

• Solutions to the OBEs—either analytically in limiting cases or numerically in general—provide direct predictions for observable quantities such as fluorescence, absorption, and dispersion.