This simulation models the phenomenon of photon echo in a gas of atoms. Each little ball represents a group of atoms with a certain speed. We know that each group interacts differently with the laser pulses, due to the Doppler shift of the light.

Once excited, each atom behaves like an electric dipole, oscillating at a given frequency. The oscillation cycles are represented in the simulation by the gray scale, where from white to black we have a variation of $2\pi$ in the cycle.

The gas of atoms is excited by a sequence of laser pulses with an area defined by

$$ \theta = \dfrac{\mu_{12}}{\hslash}\int_{-\infty}^{\infty} E_0(t)dt, $$

where $\mu_{12} = \left\langle 1 | e\hat{r} | 2 \right\rangle$ is the electronic transition matrix element and $E_0$ is the envelope of the light pulse.

In the standard photon echo configuration, the first pulse, of area $\pi/2$, leaves the gas with maximum coherence. Over time, the decoherence due to Doppler broadening ends up resetting the polarization in a half-life equal to $1/\Delta_D$, where $\Delta_D$ is the Doppler width.

The second pulse, of area $\pi$, reverses the decoherence, causing the sample to generate a pulse after a time equal to the time interval of the incident pulses.

The polarization calculation is given by

$$ P(t) = \eta\mu_{12}\int_{-\infty}^{\infty} \text{Re} [\rho_{12}(t)]g(\delta)d(\delta), $$

where $\eta$ is the atomic density, $g$ is the Maxwell-Bolztamn velocity distribution, $\delta$ is the Doppler shift and $\rho_{12}$ is one of the elements of the density matrix whose evolution is calculated by solving the Liouville-von Neumann equation:

$$ \dfrac{\partial \hat{\rho}}{\partial t} = -\dfrac{i}{\hslash} \left[ \hat{H}, \hat{\rho} \right]. $$

$\hat{H}$ is the Hamiltonian of a two-level system in the electric dipole approximation. To simplify the calculations, we assume pulses with a rectangular envelope, which greatly simplifies the solution of the Bloch equations practically without altering the polarization evolution.

Read more: