Particle transport C: Monte Carlo methods versus moment equations

This is the concluding of three articles introducing the basic concepts of the Monte Carlo method for radiation transport calculations. The demonstration calculations discussed in the previous articles give a basis for comparing the relative benefits and drawbacks of the Monte Carlo and transport equation approaches. We have seen that both strategies are based on averages over random distributions of particles whose statistical properties are known. Both methods give the same result in the limit of a large number of model particles. The main difference is that the averaging process is performed at the beginning of the transport equation calculation but at the end of the Monte Carlo solution.


  • In the transport equation approach to the two-dimensional random walk, the idea is to seek average quantities n or J and to find relationships between them (like Fick’s first and second laws). These relationships are accurate when there are large numbers of particles. To illustrate the meaning of large, note that the number of electrons in one cubic micrometer of aluminum equals 3 × 10^15. When averages are taken over such large numbers, the transport equations are effectively deterministic.
  • In the Monte Carlo method, the idea is to follow individual particles based on a knowledge of their interaction mechanisms. A practical computer simulation may involve millions of model particles, orders of magnitude below the actual particle number. Therefore, each model particle represents the average behavior of a large group of actual particles. In contrast to transport equations, the accuracy of Monte Carlo calculations is dominated by statistic variations.

An additional benefit of transport equations is that they often have closed-form solutions that lead to scaling relationships like Eq. 22 of the previous article. We could extract an approximation to the relationship from Monte Carlo results, although at the expense of some labor.

Despite the apparently favorable features of the transport equations, Monte Carlo is the primary tool for electron/photon transport. Let’s understand why. One advantage is apparent comparing the relative effort in the demonstration solutions — the Monte Carlo calculation is much easier to understand. A clear definition of physical properties of particle collisions was combined with a few simple rules. The only derivation required was that for the mean free path. The entire physical model was contained in a few lines of code. In contrast, the transport model required considerable insight and the derivation of several equations. In addition, it was necessary to introduce additional results like the divergence theorem. Most of us feel more comfortable staying close to the physics with a minimum of intervening mathematical constructions. This attitude represents good strategy, not laziness. Less abstraction means less chance for error. A computer calculation that closely adheres to the physics is called a simulation. Program managers and funding agents have a warm feeling for simulations.

Beyond the emotional appeal, there is an over-riding practical reason to apply Monte Carlo to electron/photon transport in matter. Transport equations become untenable when the interaction physics becomes complex. For example, consider the following scenario for a demonstration calculation:

In 20% of collisions, a particle splits into two particles with velocity 0.5v0 and 0.2v0. The two particles are emitted at a random angles separated by 60°. Each secondary particle has its own cross section for interaction with the background obstacles.

It would be relatively easy to modify the code of the first article to represent this history and even more complex ones. On the other hand, it would be require consider effort and theoretical insight to modify a transport equation. As a second example, suppose the medium were not uniform but had inclusions with different cross sections and with dimensions less than λ. In this case, the derivation of Fick’s first law is invalid. A much more complex relationship would be needed. Again, it would relatively simple to incorporate such a change in a Monte Carlo model. Although these scenarios may sound arbitrary, they are precisely the type of processes that occur in electron/photon showers.

In summary, the goal in collective physics is to describe behavior of huge numbers of particles. We have discussed two approaches:


  • Monte Carlo method. Define a large but reasonable set of model particles, where each model particle represents the behavior of a group of real particles with similar properties. Propagate the model particles as single particles using known physics and probabilities of interactions. Then, take averages to infer the group behavior.
  • Transport equation method. Define macroscopic quantities, averages over particle distributions. Derive and solve differential equations that describe the behavior of the macroscopic quantities.

The choice of method depends on the nature of the particles and their interaction mechanisms. Often, practical calculations usually use a combination of the two approaches. For example, consider the three types of calculations required for the design of X-ray devices (supported in our Xenos package):

  • Radiation transport in matter. Photons may be treated with the Monte Carlo technique, but mixed methods are necessary for electrons and positrons. In addition to discrete events (hard interactions) like Compton scattering, energetic electrons in matter undergo small angle scattering and energy loss with a vast number of background electrons (soft interactions). It would be impossible to model each interaction individually. Instead, averages based on transport calculations are used.
  • Heat transfer. Here, particles are the energy transferred from one atom to an adjacent one. Because the interaction model is simple and the mean-free-path is extremely small, transport equations are clearly the best choice.
  • Electric and magnetic fields. The standard approach is through the Maxwell equations. They are transport type equations, derived by taking averages over a large number of charges. On the other other hand, we employ Monte-Carlo-type methods to treat contributions to fields from high-current electron beams.


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

[2] Contact us :

[3] Field Precision home page:

Particle transport B: Monte Carlo methods versus moment equations

This is the second of three articles introducing the basic concepts of the Monte Carlo method for radiation transport calculations. In this article, we’ll consider an alternative to the Monte Carlo treatment of the two-dimensional random walk: the derivation and solution of a transport equation. Here, we define an appropriate quantity averaged over a random distribution of particles and seek a differential equation that describes how the quantity varies. For this calculation, the quantity is the average density of particles n(x,y,t) in the plane with units of number/m^2. The quantity is plotted in Fig. 1 (this figure also appeared in the preceding article). To make a direct comparison with the Monte Carlo results, we must carefully set model constraints:

  • Although the density may vary in space, the distribution of particle velocities is the same at all points. Particles all have constant speed v0 and there is an isotropic distribution of direction vectors.
  • There is a uniform-random background density of scattering objects.
  • Equation 8 of the previous article gives the probability distribution of a (the distance particles travel between collisions) in terms of the mean-free-path λ.

Figure 1. 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.


We want to find how the density changes as particles perform their random walk. Changes occur if, on the average, there is a flow of particles (a flux) from one region of space to another. If the density n is uniform, the same number of particles flow in one direction as the other, so the average flux is zero. Therefore, we expect that fluxes depend on gradients of the particle density. We can find the dependence using the construction of Fig. 2. Assume that the particle density varies in x near a point x0. Using a coordinate system with origin at x0, the first order density variation is given by Eq. 9. The goal is to find an expression for the number of particles per second passing through the line element Δy. To carry out derivation, we assume the following two conditions:


  • The material is homogeneous. Equivalently, λ has the same value everywhere.
  • Over scale length λ, relative changes in n are small.




Construction to relate flux to density gradient

Figure 2. Construction to relate flux to density gradient.

Using polar coordinates shown centered on the line element, consider an element in the plane of area (r Δθ)(Δr)$. We want to find how many particles per second originating from this region pass through Δy. We can write the quantity as the product Jx Δy, where Jx is the linear flow density in units of particles/m-s. On the average, every particle in the calculation volume has the same average number of collisions per second, given by Eq. 10. The rate of scattering events in the area element equals ν times the number of particles in the area (Eq. 11). The fraction of scattered particles aimed at the segment is given by Eq. 12.

Finally, the probability that a particle scattered out of the area element reaches the line element was given in the previous article as exp(-r/λ). Combining this expression with Eqs. 10 and 11, we can determine the current density from all elements surrounding the line segment. Taking the density variation in the form of Eq. 13 leads to the expression of Eq. 14. The integral of the first term in brackets equals zero, so that only the term proportional to the density gradient contributes. Carrying out the integrals, the linear current density is given by Eq. 15. The planar diffusion coefficient (with units m^2/s) is given by Eq. 16. Generalizing to possible variations in both x and y, can write Eq. 15 in the form of Eq. 17.  This relationship between the vector current density and the gradient of density is called Fick’s first law. Equation 18 lists Fick’s second law, a statement of conservation of particles. In the equation, the quantity ∇•J is the divergence of flux from a point and S is the source of particles at that point (particles/m^2-s). Equation 18 is the diffusion equation for particles in a plane. It states that the density at a point changes in time if there is a divergence of flux or a source or sink.

We are now ready to compare the predictions of the model with the Monte Carlo results of the previous section. Equation 19 gives is solution to the diffusion equation for particles emission from the origin of the plane. The quantity r equals √(x^2 + y^2). We can verify Eq. 19 by direct substitution by using the cylindrical form of the divergence and gradient operators and taking D as uniform in space. In order to make a comparison with the Monte Carlo calculation, we pick a time value t0 = Nc λ/v0 and evaluate A based on the condition of Eq. 20. The resulting expression for the density at time t0 is given by Eq. 21. The prediction of Eq. 21 s plotted as the solid line in Fig. 1. The results from the two methods show close absolute agreement.

Finally, we can determine the theoretical 1/e radius of the particle cloud from Eq. 21 to yield Eq. 22. In a random walk, the particle spread increases as the square root of the number of transits between collisions. For Nc = 100, the value is re/λ ≅ 14.1.


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

[2] Contact us :

[3] Field Precision home page:

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.

Random walk in a plane

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.

Probability of a collision in a differential element

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.

Exponential distribution: relation between P(a) and a

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
Angle = DTwoPi*Xi
! Random length
Length = -LOG(Xi)
! Add the vector
X = XOld + Length*COS(Angle)
Y = YOld + Length*SIN(Angle)
XOld = X
YOld = Y

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.

Random walk in a plane

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.

Particle density in a plane.

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.


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

[2] Contact us :

[3] Field Precision home page:


Field Precision educational utilities

We support several utility programs to make life easier for Windows users. They are included in all our packages and are distributed for free to the scientific community. This article reviews what’s available and also introduces a new program, Boilerplate.

The utility collection started with the FP Universal Scale. It grew out of my frustration with conventional screen rulers which were either rigidly referenced to screen pixels or absolute units like inches or centimeters. A more useful approach is to reference the ruler to the units of the graph or photograph to be measured. Accordingly, after several years of thought I set out to create an on-screen version of the much-loved Gerber Variable Scale. The implementation involved intensive interactions with the Windows API, so I decided to use RealBasic with purchased plugins to handle screen overlays. During development, the program expanded from a simple screen ruler to a complete screen digitization system for scientists and engineers.

There were three motivations for the next utility, the FP File Organizer:

  • In comparison to sophisticated two-window file managers like Free Commander, I wanted a simple, clean interface that supported the functions I used every work day.
  • Our technical programs involve extended file organization. In discussing file management in tutorials, I wanted a standard reference environment.
  • I needed a general file-manager unit for my MIDI programs.

FP File Organizer has several nice features like fast file searches, full path copy to the clipboard, definable tools, special folders and desktop shortcut creation. I use the program for all my work except for multi GB file transfers. For these, I use xcopy or robocopy.

The Cecil_B program converts an organized set of BMP files into an AVI movie. I developed it in response to a customer request to make animations of solutions in time-domain programs like TDiff and HeatWave. I created the final two utilities, Computer Task Organizer (CTO) and Boilerplate to reduce frustrations I noticed over the last 30 years using Windows. With regard to CTO, I found that most of my work day involved running the same programs with the same documents or going to the same website repetitively. The program reduces the 100 tasks that I perform every day to single button clicks.

The new utility Boilerplate (Figure 1) expands the functions of the Windows clipboard in two ways:

  • You can build a library of standard text selections (i.e., boilerplate) that can be transferred to the clipboard with a single button press — ready to paste into a document.
  • You can recall items previously on the clipboard.

The second feature deals with an irritating limit of the clipboard — it stores only one item at a time. Boilerplate keeps a running record of the last twenty clipboard texts — they can be recalled to the clipboard with a single button click. I got the idea from the old utility Clipboard Magik. The program had a lot of potential, but was difficult to utilize in practice.

Boilerplate screenshot

Figure 1. Boilerplate screenshot.


[1] Contact us :

[2] Field Precision home page:

Magnum for MRI facility shielding

I recently had an inquiry from the president of a company that creates magnetic shielding for MRI facilities. In the past, the company installed iron sheets in response to customer specifications. The company’s current goal is to assume more design responsibilities by developing an in-house computational capability. Figure 1 shows a typical MRI facility along with calculated contours of |B| supplied by the one of company’s customers. A requirement is that |B| must be less than 5 gauss (0.0005 tesla) everywhere outside the central treatment room.


Example of an MRI facility

Figure 1. Example of an MRI facility.

The complicating factor is that the inner workings of the MRI magnet are proprietary, so my contact would not know the geometry of the drive currents and iron poles that generate the field. He would have to base his shielding calculations entirely on fringing field patterns supplied by his customer (like that of Fig. 1). Everything inside the inner line would be a mystery. My contact was being pursued by a sales rep for a well-known alternative to Magnum. I won’t name names, but for the sake of discussion let’s refer to the program as Lucia di Lammermoor (LDL). The sales rep felt he had the perfect solution. On top of the high price of LDL, my contact could buy a special inverse-solution add-on that would determine the unknown magnet configuration from the fringing-field pattern. There are two drawbacks to this approach:

  • James Clerk Maxwell says it’s impossible.
  • It’s totally unnecessary.

The critical insight is that the fringing fields of any solenoid assembly (no matter how complex) approach those of a simple magnetic dipole in the region outside the assembly. This tutorial reviews the theory:

Magnetic Dipole Moment of a Coil Assembly

To emphasize the point, Figure 2 shows contours of |B| calculated by Magnum for a current loop of radius R = 0.5 m carrying current I  = 1000.0 A (the plot plane includes the magnet z axis). The line shapes clearly resemble those of Fig. 1. Complicating the comparison  is the fact that the lines and intervals of Fig. 1 (supplied by my contact’s customer) are physically impossible. There are three possible explanations:

  • LDL gave the wrong answer.
  • The LDL user at the magnet manufacturer did not pay attention to the potentially large effects of computation boundaries on the weak fringing fields.
  • The magnet manufacturer was reluctant to send actual data to my contact, so they had a draftsman create them.



Fringing fields of a circular coil

Figure 2. Fringing fields of a circular coil calculated by Magnum.

Supposing that some day my contact receives a PDF document with real data, here’s how the analysis would proceed.

1) Use the Universal Scale (shown in Fig. 1) to measure the locations (zi) of contours |Bi| along the z axis and an axis normal to z passing through the magnet center (xi).

2) Calculate estimates of the magnet dipole moment from the equations

mi = Bi*zi^3/(μ0/2*π),  mi = 2*Bi*xi^3/(μ0/2*π).

For a valid field distribution, all of the estimated values should be close, giving an average value m.

3) Set up a Magnum solution volume (large compared to the diagnostic room to minimize boundary effects). For the applied field, define a circular coil normal to z at the origin with a radius R comparable to that of the magnet assembly. Assign the coil current according to I = m/(π*R^2) to replicate the fringing field pattern.

4) Add shielding walls as needed, and run a standard Magnum calculation. Analyze the 5 gauss contour to make sure it is everywhere within the diagnostic room.

In the calculation, scaling relationships may be applied to deal with thin sheets of iron following the discussion in this tutorial:

To facilitate the application, we have added two features to Magnum:

  • We have modified the contour-line plot routines in MagView to enable users to enter a set of specific values (e.g., the 5 gauss limit).
  • We have doubled the number of entries in the library of soft magnetic materials supplied with the code to include common shielding materials like M36.

Magnetic material data are available at


[1] You can get information on Magnum at

[2] Contact us :

[3] Field Precision home page: