The force in the $j$ direction exerted on the atoms due to the pair of beams in that direction is given by

$$
F_j = F_j^+ + F_j^-,
$$

where

$$
F_j^\pm = \pm \hslash k \dfrac{\gamma_{22}}{2} \times
$$
$$
\times \dfrac{ I/I_0 }{ 1 + I/I_0 + 4(\delta \mp kv_j \mp \beta_j r_j)^2/\gamma_{22}^2 }.
$$

$\gamma_{22}$ is the rate of spontaneous decay of atoms from state |2⟩ to state |1⟩,
$k$ is the wave number, $I$ is the intensity of the laser beam,
$I_0$ is the saturation intensity of the transition |1⟩ - |2⟩,
$\delta$ is the laser detuning with respect to the transition resonance,
$v$ is the speed of the atom and $\beta_j$ is defined by

$$
\beta_j = \dfrac{g\mu_B}{\hslash} \dfrac{\partial B_j}{\partial j}.
$$

$g$ is the gyromagnetic factor, $\mu_B$ is the Bohr magneton and $B_j$ the $j$ component
of the magnetic field generated by the pair of Helmholtz coils.

The dynamics of atoms is based on Newton's second law and a random variable $\epsilon = 0.5$
m/s was added in order to simulate the random fluctuations in atomic velocities due to spontaneous emissions:

$$
v(t) \rightarrow v(t) + \epsilon.
$$

Read more: