Cinco: Método de Monte Carlo en simulación de experimentos

La complejidad de situaciones experimentales hacen necesario un análisis detallado de un posible experimento, especialmente si la geometría y componentes son variados. El método de Monte Carlo nos permite abordar estos casos en forma general.
Ejemplo: consideremos un neutrón de 14 MeV incidente en un material dado $H_2O.$ Nos preguntamos: cuán probable es capturar (a energías térmicas) el neutrón debido a la reacción

$$\begin{eqnarray*} n + p &\rightarrow& d + \gamma \ \ (E_{\gamma} = 2.2 \ \ MeV). \\ \end{eqnarray*}$$

Hagamos una lista ("historia") de cómo esto puede suceder.
1) Obtenga $<L>(E)$, el camino medio libre a energía $E: \ \ \frac{N(x)}{N_o}&=& \exp^{-\frac{x}{\overline L}}$
2) Decida si colisión en $H_2O$ es en $H$ o en $O$ (de acuerdo a respectivas $\sigma_{tot}$ )
3) Obtenga ángulo $\theta$ de la correspondiente distribución angular $d\sigma/d\Omega$
4) Calcule nueva energía del neutrón usando cinemática de dos cuerpos.
5) Vea si $E$ es menor que la energía t'ermica, $\frac{3}{2}kT = .025 \ \ eV$
a)Sí. Continúe en 6).
b) No. Vuelva a 1).

6) Propague el neutrón térmico comparando la sección eficaz de captura con la elástica hasta que se absorba (emitiendo un $\gamma$ de 2.2 MeV) o se escape del volumen.

De este ejemplo se ve que es necesario tomar decisiones. Estas decisiones se hacen al azar (de ahí el nombre "Monte Carlo"), basadas en información previa del fenómeno físico, es decir, en información sobre la distribución involucrada. Pero los computadores pueden generar una variable $r$ distribuída uniformemente en $0<r<1$. Cómo obtenemos de $r$ la correspondiente variable al azar de interés?

TEOREMA / FUNDAMENTAL PARA MUESTREO.-
Se tiene la distribución de probabilidad $\frac{dP(x)}{dx}$, es decir,
$dP(x)$ = probabilidad que $x$ ocurra entre $x$ y $x+dx$.
Supongamos, por generalidad, que $x$ está limitada al intervalo ($x_{min}$, $x_{max}$).
Entonces, la distribución acumulativa normalizada es una distribución uniforme:

$$\begin{eqnarray*} \frac{\int_{x_{min}}^x dP(t)} {\int_{x_{min}}^{x_{max}} dP(t)} &=& r \ \ \ \ \ \ (1) \\ \end{eqnarray*}$$

Véase Fig. del Green que ilustra la equipartición de áreas implícita en eq. (1)

No demostraremos este Teorema, que para algunos puede ser "manifiestamente verdadero". Notemos, por lo menos, que para $x=x_{min}, r=0$, y para $x=x_{max}, r=1$, es decir, la distribución (1) satisface correctamente los límites correspondientes a la variable $r$.

El Teorema nos dice que el valor buscado $x$ se obtiene resolviendo ("invertiendo") eq. (1) para $x$.

Veamos a continuación casos de interés físico en que (1) se puede resolver analíticamente para $x$ en función de $r$.

1) La distribución es uniforme en un intervalo, por ejemplo, el caso del ángulo azimutal $0<\phi<2\pi$ en scattering

$$\begin{eqnarray*} \frac{\int_{0}^x dt} {\int_{0}^{2\pi} dt} &=& r \ \ \ \ \ \ (... ... 2\pi r, \ \ \ \ r \ \ uniforme \ \ entre \ \ 0 \ \ y \ 2\pi \\ \end{eqnarray*}$$

2) Angulo polar $\theta$para scattering de Rutherford para ángulos pequños.

$$\begin{eqnarray*} \frac{dP(\theta)} {d\theta} &=& K \frac{1}{\theta^4}, \ \ o, ... ...r(x_{max}^{\alpha+1} - x_{min}^{\alpha+1})}]^\frac{1}{\alpha+1} \end{eqnarray*}$$

Nótese que
a) El factor multiplicativo de la distribución no es relevante en la simulación
b) Para el caso $\alpha=0$, se reduce a la distribución uniforme, como debe ser.

3) Camino medio libre $\Lambda$ (o también decaimiento radioactivo). Estos casos corresponden a una distribución exponencial $\frac{dP(x)}{dx} = K \exp{\frac{-x}{\Lambda}}$:

$$\begin{eqnarray*} \frac{\int_{x_{min}}^x K e^-(t/\Lambda) dt} {\int_{x_{min}}... ...^{-\frac{x_{max}}{\Lambda}} -e^{-\frac{x_{min}}{\Lambda}})] \\ \end{eqnarray*}$$

4) Distribución angular isotrópica.- $\frac{d\sigma}{d\Omega}=const$, es decir, uno desea enviar el mismo número de partículas (en promedio) en cada elemento de ángulo sólido.

$$\begin{eqnarray*} d\Omega = sin \theta d\theta d\phi = d(cos \theta) d\phi \end{eqnarray*}$$

Hemos expresado el elemento de ángulo sólido $d\Omega$ como el producto de dos elementos diferenciales, cada uno independiente y con una distribución uniforme. Para $0<\phi<2\pi$ y $\theta_{min}<\phi<\theta_{max}$, obtenemos, por 1):

$$\begin{eqnarray*} cos\theta &=& cos \theta_{min} + r_1 (cos \theta_{max}-cos \theta_{min}) \\ \phi &=& 2 d\pi r_2, \\ \end{eqnarray*}$$

en que $r_1$ y $r_2$ son dos números al azar y uniformes en el intervalo (0,1).

Si la distribución a ser simulada no tiene inversa, es decir, si no se puede resolver $x$ en función de $r$, es posible siempre usar un método de "trial and error" que eventualmente elije un valor de $x$ con la distribución deseada. Este método se lo conoce con el nombre de "acceptance/rejection method". Véase el Apéndice de Green o el libro de propiedades de las partículas.

El problema de "random walk" (en 3 dimensiones).-
Como se sabe, para el caso de un paso constante $d$, la solución analítica para el promedio del cuadrado de la distancia final desde el origen es:

$$\begin{eqnarray*} \overline {r^2} &=& N d \\ d &=& longitud \ \ del \ \ paso \ ... ...derar \\ r &=& \'ultimo \ \ radio \ \ desde \ \ el \ \ origen. \end{eqnarray*}$$

Procedimiento utilizando Monte Carlo.-
Se parte del origen.
Se elije una dirección $\theta_1,\phi_1$ isotrópica ($cos \theta$ y $\phi$ uniformes en [-1,1] y [0,2$\pi$], respectivamente.
Se avanza una distancia $d$ en la dirección $\overrightarrow{01}$
Se elije una nueva dirección $\theta_2,\phi_2$ con respecto a $\overrightarrow{01}$. Se recomienda usar cosenos directores en esta parte. Las fórmulas para llevar a cabo esto son generales y envuelven sólo consideraciones geométricas.
Se avanza una distancia $d$ en la dirección $\overrightarrow{12}$
Etc., hasta que se llega al valor N, por ejemplo N=10. Se calcula la distancia al origen.

Se repite este procedimiento muchas veces y se obtiene así el promedio $\overline {r^2}$. La idea es comparar el valor obtenido con el dado por solución analítca, $\overline {r^2} &=& N d$.

Vovamos al caso de neutrones de 14 MeV. El procedimiento es similar al de random walk, excepto que el paso $d=\Lambda$, el camino medio libre entre colisiones, varía, ya que la energía del neutrón se degrada en cada colisión lo que implica que la sección eficaz total cambia a su vez.