# Particle transport A: Monte Carlo methods versus moment equations

This is the first of three articles giving a brief introduction to Monte Carlo methods applied to electron/photon transport, the foundation of our GamBet package for X-ray science. In particular, I will use a simple example to show how Monte Carlo compares to transport equations (e.g., the diffusion equation). The fundamental issue is how to deal with extremely large numbers of objects. Calculating the history of every particle is beyond the capabilities of even the most powerful computers. Instead, we seek the average properties of large groups.

In the Monte Carlo method, the full set of particles is represented by a calculable set of model particles. In this case, each model particle represents a group. We follow detailed histories of model particles as they undergo random events like collisions with atoms. Characteristically, we use a random-number generator with a known probability distribution to determine the outcomes of the events. In the end, the core assumption is that averages over model particles represent the average behavior of the entire group. The alternative to this approach is the derivation and solution of moment (or transport) equations. The following article covers this technique.

Instead of an abstract discussion, we’ll address a specific example to illustrate the Monte Carlo method. Consider a random walk in a plane. As shown in Fig. 1, particles emerge from a source at the origin with uniform speed v0. They move freely over the surface unless they strike an obstacle. The figure represents the obstacles as circles of diameter w. The obstacles are distributed randomly and drift about so we can never be sure of their position. The velocity of obstacles is much smaller than v0. If a particle strikes an obstacle, we’ll assume it bounces off at a random direction with no change in speed. The obstacles are unaffected by the collisions.

Figure 1. Random walk in a plane.

In a few sentences, we have set some important constraints on the physical model:

• The nature of the particles (constant speed v0).
• The nature of the obstacles (diameter w, high mass compared to the particles),
• The nature of the interaction (elastic collision with isotropic emission from the collision point)

The same type of considerations apply to calculations of radiation transport. The differences are that 1) the model particles have the properties of photons and electrons, 2) the obstacles are the atoms of materials and 3) there are more complex collision models based on experimental data and theory. To continue, we need to firm up the features of the calculation. Let’s assume that 10^10 particles are released at the origin at time t = 0. Clearly, there are too many particles to handle on a computer. Instead, we start Np = 10,000 model particles and assume that they will give a good idea of the average behavior. In this case, each model particle represent 10^6 real particles. We want to find the approximate distribution of particle positions after they make Nc collisions. The logic of a Monte Carlo calculation for this problem is straightforward. The first model particle starts from the origin moving in a random direction. We follow its history through Nc collisions and record its final position. We continue with the other Np – 1 model particles and then interpret the resulting distribution of final positions.

The source position is x = 0, y = 0. To find the emission direction, we use a random number generator, a component of all programing languages and spreadsheets. Typically, the generator returns a random number ? equally likely to occur anywhere over the interval of Eq. 1. Adjusting the range of values to span the range 0 ? 2?, the initial unit direction vector is given by Eq. 2.

The particle moves a distance a from its initial position and then has it first collision. The question is, how do we determine a? It must be a random quantity because we are uncertain how the obstacles are lined up at any time. In this case, we seek the distribution of expectations that the particle has a collision at distance a, where the distance may range from 0 to ?. To answer the question, we’ll make a brief excursion into probability theory.

Let P(a) equal the probability that the particle moves a distance a without a collision with an object. By convention, a probability value of 0.0 corresponds to an impossible event and 1.0 indicates a certain event. Therefore, P(0) = 1.0 (there is no collision if the particle does not move) and P(?) = 0.0 (a particle traveling an infinite distance must encounter an object). We can calculate P(a) from the construction of Figure 2. The probability that a particle reaches a + ?a equals the probability that the particle reaches a times the probability that it passes through the layer of thickness ?a without a collision. The second quantity equals 1.0 minus the probability of a collision.

Figure 2. Probability of a collision in a differential element.

To find the probability of a collision in the layer, consider a segment of height h. If the average surface density of obstacles is N particles/m^2, then the segment is expected to contain Nh ?a obstacles. Each obstacle is a circle of diameter w. The distance range for an interaction with an obstacle is called the cross-section ?. In this case, we will associate the interaction width with the obstacle diameter, or ? = w. The fraction of the height of the segment obscured by obstacles is given by Eq. 3. The exit probability is given by Eq. 4.

A first-order Taylor expansion (Eq. 5) leads to Eq. 6. Equation 6 defines another useful quantity, the macroscopic cross section ? = n ? with dimensions 1/m. Solving Eq. 6 leads to Eq. 7. The new quantity in Eq. 7 is the mean free path, ?. It equals the average value of a for the exponential probability distribution. The ideas of cross section, macroscopic cross section and mean-free path are central to particle transport.

We can now solidify our procedure for a Monte Carlo calculation. The first step is to emit a particle at the origin in the direction determined by Eq. 2. Then we move the particle forward a distance a consistent with the probability function of Eq. 7. One practical question is, how do we create an exponential distribution with a random number generator that produces only a uniform distribution in the interval of Eq, 1? The plot of the probability distribution of Eq. 7 in Fig. 3 suggests a method. Consider the 10% of particles with collision probabilities between P(0.3) and P(0.4). The corresponding range of paths extends from a(0.6)/? = -ln(0.4) = 0.9163 to – \ln(0.3) = 1.204. If we assign path lengths from the uniform random variable according to Eq, 8, then we can be assured that on the average 10% will lie in the range a/? = 0.9163 to 1.204. By extension, if we apply the transformation of Eq. 8 to a uniform distribution, the resulting distribution will be exponential. To confirm, the lower section of Fig. 3 shows a random distribution calculation with 5000 particles.

Figure 3. Exponential distribution. a) Relationship between intervals of P(a) and a. b) Testing assignment of the collision distance, 5000 particles with mean-free-path equal to 1.0.

To continue the Monte Carlo procedure, we stop the particle at a collision point a distance a from the starting point determined by Eq. 8 and then generate a new random number ? to determine the new direction according to Eq. 2. Another call to the random-number generator gives a new propagation distance a from Eq. 8. The particle is moved to the next collision point. After Nc events, we record the final position and start the next particle. The simple programing task with the choice ? = 1 is performed by the following code:

DO Np=1,NShower ! Start from center XOld = 0.0 YOld = 0.0 ! Loop over steps DO Nc=1,NStep ! Random direction CALL RANDOM_NUMBER(Xi) Angle = DTwoPi*Xi ! Random length CALL RANDOM_NUMBER(Xi) Length = -LOG(Xi) ! Add the vector X = XOld + Length*COS(Angle) Y = YOld + Length*SIN(Angle) XOld = X YOld = Y END DO END DO 

Figure 4 shows the results for ? = 1 (equivalently, the plot is scaled in units of mean-free-paths). The left-hand side shows the trajectories of 10 particles for Nc = 100 steps. With only a few particles, there are large statistical variations, making the distribution in angle skewed. We expect that the distribution will become more uniform as the number of particles increases because there is no preferred emission direction. The right-hand side is a plot of final positions for Np = 10,000 particles. The distribution is relatively symmetric, clustered within roughly 15 mean-free-parts of the origin. In comparison, the average total distance traversed by each particle is 100.

Figure 4. Random walk in a plane with mean-free-path equal to 1.0 and 100 collisions. a) Sample trajectories for 10 model particles. b) Final positions of 10,000 model particles.

Beyond the visual indication of Fig. 4, we want quantitative information about how far particles move from the axis. To determine density as a function of radius, we divide the occupied region into radial shells of thickness ?r and count the number of final particle positions in each shell and divide by the area of the shell. Figure 5 shows the results. The circles indicate the relative density of particles in shells of width 0.8?. Such a plot is called a histogram and the individual shells (containers) are called bins. Histograms are one of the primary methods of displaying Monte Carlo results. Note that the points follow a smooth variation at large radius, but that they have noticeable statistical variations at small radius. The reason is that the shells near the origin have smaller area, and therefore contain fewer particles to contribute to the average. Statistical variations are the prime concern for the accuracy of Monte Carlo calculations.

Figure 5. Particle density as a function of radius (distance from the source), with mean-free-path equal to 1.0 and 100 collisions. The solid line is the solution of the two-dimensional diffusion equation and the points are the results of the Monte Carlo solution.

Footnotes

[1] Use this link for a copy of the full report in PDF format: Monte Carlo method report.