Classically, one can think of radiation trapping as a multiple-scattering phenomena, where a photon is scattered by multiple atoms in a cloud. This motivates treatment as a diffusion problem. As such, one can primarily consider the mean free path of light, defined as the reciprocal of the density of scatterers and the scattering cross section:

${\displaystyle \ell _{\text{mf}}={\frac {1}{\rho \sigma _{\text{sc}}}}.}$ 

One can assume for simplicity that the scattering diagram is isotropic, which ends up being a good approximation for atoms with equally populated sublevels of total angular momentum. In the classical limit, we can think of the electromagnetic energy density as what is being diffused. So, we consider the diffusion constant in three dimensions,

${\displaystyle D={\frac {\ell _{\text{mf}}^{2}}{3\tau _{r}}},}$ 

where ${\displaystyle \tau _{r}}$  is the transport time.^[3] The transport time accounts for both the group delay between scattering events and Wigner's delay time, which is associated with an elastic scattering process.^[4] It is written as

${\displaystyle \tau _{r}={\frac {\ell _{\text{mf}}}{\nu _{\text{g}}}}+\tau _{\text{W}},}$ 

where ${\displaystyle \nu _{\text{g}}}$  is the group velocity. When the photons are near resonance, the lifetime of an excited state in the atomic vapor is equal to the transport time, ${\displaystyle \tau _{r}=\tau _{at}}$ , independent of the detuning.^[5] This comes in handy, since the average number of scattering events is the ratio of the time spent in the system to the lifetime of the excited state (or equivalently, the scattering time). Since in a 3D diffusion process the electromagnetic energy density spreads as ${\displaystyle \langle r^{2}\rangle =6Dt}$ , we can find the average number of scattering events for a photon before it escapes:

${\displaystyle \langle N_{\text{sc}}^{2}\rangle ={\frac {\langle r^{2}\rangle }{6D\tau _{at}}}.}$ 

Finally, the number of scattering events can be related to the optical depth ${\displaystyle b}$  as follows. Since ${\displaystyle {\sqrt {\langle r^{2}\rangle }}\sim b\ell _{\text{mf}}}$ , the number of scattering events scales with the square of the optical depth.^[6]

In 1947, Theodore Holstein attacked the problem of imprisonment of resonance radiation in a novel way. Foregoing the classical method presented in the prior section, Holstein asserted that there could not exist a mean free path for the photons. His treatment begins with the introduction of a probability function ${\displaystyle G(\mathbf {r} ,\mathbf {r} ')\,d\mathbf {r} }$ , which describes the probability that a photon emitted at ${\displaystyle \mathbf {r} }$  is absorbed within the volume element ${\displaystyle d\mathbf {r} }$  about the point ${\displaystyle \mathbf {r} }$ . Additionally, one can enforce atom number conservation to write

${\displaystyle A-B=dt\,d\mathbf {r} \,{\frac {\partial n(\mathbf {r} )}{\partial t}},}$ 

where ${\displaystyle A,B}$  represent the number increase and decrease in population of excited atoms, and ${\displaystyle n(\mathbf {r} )}$  is the number density of excited atoms. If the reciprocal lifetime of an excited atom is given by ${\displaystyle \Gamma }$ , then ${\displaystyle B}$  is given by

${\displaystyle B=\Gamma n(\mathbf {r} )\,d\mathbf {r} \,dt.}$ 

Then ${\displaystyle A}$  is obtained by considering all other volume elements, which is where the introduction of ${\displaystyle G(\mathbf {r} ,\mathbf {r} ')}$  becomes useful. The contribution of an outside volume ${\displaystyle d\mathbf {r} '}$  to the number of excited atoms is given by the number of photons emitted by that outside volume ${\displaystyle d\mathbf {r} '}$  multiplied by the probability that those photons are absorbed within the volume ${\displaystyle d\mathbf {r} }$ . Integration over all outside volume elements yields

${\displaystyle A=\Gamma \,dt\,d\mathbf {r} \,\int d\mathbf {r} '\,n(\mathbf {r} ')G(\mathbf {r} ,\mathbf {r} ').}$ 

Substituting ${\displaystyle A}$  and ${\displaystyle B}$  into the particle conservation law, we arrive at an integral equation for the density of excited atoms – the Holstein equation^[7]

${\displaystyle {\frac {\partial n(\mathbf {r} )}{\partial t}}=-\Gamma n(\mathbf {r} )+\Gamma \int d\mathbf {r} '\,n(\mathbf {r} ')G(\mathbf {r} ,\mathbf {r} ').}$ 

Finding the escape probability of photons from the Holstein equation edit

Now to find the escape probability of the photons, we consider solutions by ansatz of the form

${\displaystyle n(\mathbf {r} ,t)=n(r)e^{\beta }.}$ 

Observing the Holstein equation, one can note that these solutions are subject to the constraint

${\displaystyle (1-\beta /\Gamma )n(\mathbf {r} )=\int d\mathbf {r} '\,n(\mathbf {r} ')G(\mathbf {r} ,\mathbf {r} ').}$ 

Aided by the exchange symmetry of ${\displaystyle G}$ , namely that ${\displaystyle G(\mathbf {r} ,\mathbf {r} ')=G(\mathbf {r} ',\mathbf {r} )}$ , one can use variational methods to assert that ${\displaystyle \delta (\beta /\Gamma )=0}$  leads to

${\displaystyle {\frac {\beta }{\Gamma }}=1-{\frac {\iint d\mathbf {r} \,d\mathbf {r} '\,G(\mathbf {r} ,\mathbf {r} ')n(\mathbf {r} )n(\mathbf {r} ')}{\int d\mathbf {r} \,n^{2}(\mathbf {r} )}}.}$ 

Completing the square and introducing the escape probability ${\displaystyle E(\mathbf {r} )\equiv 1-\int d\mathbf {r} '\,G(\mathbf {r} ,\mathbf {r} ')}$ , whose definition follows from that all particles must either be absorbed or escape with a summed probability of 1, an equation in terms of the escape probability is derived:

${\displaystyle {\frac {\beta }{\Gamma }}={\frac {\int d\mathbf {r} \,n^{2}(\mathbf {r} )E(\mathbf {r} )+{\frac {1}{2}}\iint d\mathbf {r} \,d\mathbf {r} '\,[n(\mathbf {r} )-n(\mathbf {r} ')]^{2}G(\mathbf {r} ,\mathbf {r} ')}{\int d\mathbf {r} \,n^{2}(\mathbf {r} )}}.}$ 

Many contemporary studies in atomic physics utilize numerical solutions to Holstein's equation to both show the presence of radiation trapping in their experimental system and to discuss its effects on the atomic spectra. Radiation trapping has been observed in a variety of experiments, including in the trapping of cesium atoms in a magneto-optical trap (MOT), in the spectroscopic characterization of dense Rydberg gases of strontium atoms, and in lifetime analyses of doped ytterbium(III) oxide for laser improvement.^[8][9][10]

To solve or simulate the Holstein equation, the Monte Carlo method is commonly employed. An absorption coefficient is calculated for an experiment with a certain opacity, atomic species, Doppler-broadened lineshape, etc., and then a test is made to see whether the photon escapes after ${\displaystyle n}$  flights through the atomic vapor (see Figure 1 in the reference).^[11]

Other methods include transforming the Holstein equation into a linear generalized eigenvalue problem, which is more computationally expensive and requires the usage of several simplifying assumptions, including but not limited to that the lowest eigenmode of the Holstein equation is parabolic in shape, the atomic vapor is spherical, the atomic vapor has reached a steady state after the near-resonant laser has been shut off, etc.^[8]