The previous page showed how to determine photon path lengths and scattering angles using the beam attenuation , the scattering phase function , and a uniform random number generator. This page shows how to combine those two processes to create random photon paths through an absorbing and scattering medium.
There is more than one way to simulate photon paths, and each will give the same answer. However, some techniques can be numerically much more eﬃcient than others. Indeed, a reasonable approach to developing a Monte Carlo algorithm for a particular problem is to
- ﬁrst ﬁgure out how to numerically simulate a process as it occurs in nature, and
- then ﬁgure out how to simulate another, perhaps artiﬁcial, process that will give the same answer as the ”natural” process, but with less computational time.
This page illustrates this two-step development process.
Consider ﬁrst how photons propagate through a medium. Loosely speaking, a photon travels until it interacts with a particle, e.g. a molecule of water or chlorophyll. It is then either absorbed by the particle and disappears, or it is scattered into a new direction and continues on its way until in interacts with another particle.
Recall the albedo of single scattering, . If there is no absorption, and . If there is no scattering, and . thus can be interpreted as the probability of photon survival in any particular interaction. When a photon encounters a particle, we can randomly decide if the photon is to be absorbed or scattered as follows:
- Draw a random number from a distribution.
- If , then the photon is scattered.
- If , then the photon is absorbed.
If the photon is absorbed, tracing stops and a new photon is emitted from the source and tracing begins anew. If the photon is scattered, two new random numbers are drawn and used to determine new polar and azimuthal scattering directions and as shown on the previous page. Another random number is used along with to determine the distance traveled before another interaction.
Figure 1 illustrates this process for two photons, which also introduces the geometry to be used in numerical simulations below. A source emits a collimated beam of photons, which are then recorded by an annular, ring, or ”bullseye” detector some distance away. The red photon is emitted by the source, undergoes one scattering, and is then absorbed by a particle. The green photon is emitted, undergoes two scatterings, and is recorded by a detector.
This process mimics what happens in nature. We will call this ”type 1” photon tracing (there are no standard names for ways of tracing photons). Note that all of the computations used to trace the absorbed photon are wasted because the photon never reached the detector. Nature can aﬀord to trace innumerable photons and waste some by absorption, but that is not advisable for most numerical simulations. We therefore seek other ways to trace photons.
The previous page showed that the mean free path or average distance between interactions with the medium is . These interactions can result in either absorption or scattering of the photon, as just described. Rather than tracing one photon at a time as nature does, consider a source emitting ”packets” of many photons. Then view each interaction as having a fraction of the photons in the packet being absorbed, and the remaining fraction being scattered, all in the same direction. Let the packet be emitted with an initial weight of , which can represent one unit of energy, power, or some number of photons. At each interaction, the current weight is multiplied by to account for the loss of energy or number of photons by absorption (that is, a fraction continues onward). The scattered packet then carries a reduced weight. If the photon packet reaches the target, the current weight is tallied. Another packet is then emitted from the source and traced. This tracing process, which we’ll call ”type 2,” is illustrated by the green photon track in Fig. 2. After two scatterings, as illustrated, the detected photon packet has weight .
A third ray-tracing process can be envisioned. The mean distance traveled between scattering events is . We can thus use and a random number to determine the distance between scattering interactions, and the initial weight of is not changed at each interaction because all photons in the packet are viewed as being scattered. Then, if the photon packet reaches the target, absorption is treated as a continuous process occurring along the entire photon path. Assuming homogeneous water, the ﬁnal weight tallied is then , where is the total path length in meters and is the absorption coeﬃcient. The red track in Fig. 2 illustrates this ”type 3” ray tracing. The red track shows a total path length of , so the ﬁnal weight is .
The two tracks in Fig. 2 are drawn as though each track were generated by exactly the same sequence of random numbers. Because , the individual type 3 photon paths will be greater than the type 2 paths. The scattering angles are the same. Thus these two tracing types clearly lead to diﬀerent results, photon packet by photon packet. However, numerical simulation of many photon packets shows all three of these photon tracing types yields the same distribution of energy at the detector.
To summarize, the three types of ray tracing considered here are
- Type 1:
- Individual photons are tracked, and photons can be absorbed.
- Type 2:
- Photon packets are tracked, with a packet weight being multipled by at each interaction.
- Type 3:
- Photon packets are tracked, with track lengths determined by the mean free path for scattering, no weighting at scattering events, and absorption treated as a continuous process based on total path length.
To illustrate the results obtained for diﬀerent ways of tracking photons, a Monte Carlo code was written to simulate the energy received by an annular target as illustrated in Figs. 1 and 2. For the simulations shown here, the IOPs were deﬁned by , , and a Fournier-Forand phase function with parameter values = (index of refraction, slope of Junge distribution) chosen to give a good ﬁt to the Petzold average particle phase function. Thus , and optical distance is numerically the same as geometrical distance in meters. A run was made with photons being sent from the source and using type 1 ray tracing. Photons were traced until they were absorbed. Figure 3 shows some of the resulting statistics.
The red histogram shows the percent of photons that traveled an optical distance between interactions, for a bin size of . The theoretically expected fraction of photons traveling an optical distance between and between interactions is
The red dots in the ﬁgure are the expected values given by this formula. The shortest photon path length between interactions was and the longest was 16.69. The blue histogram shows the distribution of total distances traveled until the photons were absorbed. Thus the value for the ﬁrst bar shows that about 18% of the photons were absorbed after going a total optical distance between 0 and 1. Note that this distance can represent more than one interaction, i.e., a photon being scattered one or more times before being absorbed. Both the red and the blue histograms sum to 100%. As shown on the previous page, the mean distance traveled between interactions is , or in meters. For this particular simulation the actual average was (or 0.9974 m for these IOPs). The small diﬀerence is statistical noise resulting from the ﬁnite sample size of the numerical simulation. Likewise, the average distance traveled until the photons are absorbed is . For the present case of this gives 5 m. The average for this simulation was 4.9949 m. Since for this simulation, another way to view this is that the photons were on average scattered four times before being absorbed on the ﬁfth interaction.
We next compare results for the three diﬀerent ways of tracking photons. Because oceanic phase functions scatter much more light at small scattering angles than at large angles, most photons that are scattered just a few times will hit the detector near its center. To even out the numbers of photons (or power) detected by each ring, an annular target was deﬁned with a logarithmic spacing for the radii of the detector rings. A logarithmic spacing is often used in instruments so that each detector ring receives roughly equal amounts of power, which reduces the dynamic range needed for the instrument design. The detector simulated here had rings with the smallest ring radius being and the largest being . This detector is placed in a target plane some distance from the source and centered on the optical axis of the source, as shown in Figs. 1 and 2. The photons crossing the detector plane at some distance from the detector center are tallied in bins as follows:
- Bin :
- Unscattered photons that hit the detector at .
- Bin :
- Scattered photons that hit the target plane inside the ﬁrst detector ring, i.e. at .
- Bins :
- Photons that hit detector rings 1 to .
- Bin :
- Photons that hit the detector plane outside the outer ring, i.e. at .
Simulations were made with photons (for type 1) or photon packets (types 2 and 3) emitted from the collimated source. Figure 4 shows the distribution of photons or packets received anywhere in the detector plane as a function of the number of scatterings, for the three types of tracing and for target plane distances of and 15. The left panel shows that for one or two percent of photons (depending on the tracing type) reach the target plane without being scattered. Most photons are scattered 3 or 4 times, and very few photons are scattered more than 10 times. The right panel shows that for , almost no photons reach the target plane unscattered, and most undergo 5 to 25 scatterings, with a peak around 10 or 15, depending on the way the photons are traced. For type 1 photon tracing, almost no photons are scattered more than 30 times. Note that the probablity of surviving 30 scatterings is . For tracing types 2 and 3, which never have photon packets absorbed, there are broad tails in the number of scatterings, although for type 2 a packet being scattered 40 times has an almost negligible weight of .
The left panel of Fig. 5 shows the intersection points of the photons reaching a area of the target plane for the ﬁrst emitted photons, tracing type 1, and the target plane at . Note that of the emitted photons, only 2973 reached the target plane (of which 2966 are in the area plotted). This is less than 30% of the photons making a contribution to the answer of how much power is detected; i.e., 70% of the calculations were wasted. The photon-intersection dots are color-coded to show the number of times each plotted point was scattered. The two black circles show oﬀ-axis angles of and . Very few photons were scattered more than about oﬀ of the optical axis. The right-hand panel shows the distributions of the photons reaching the target plane by number of scatterings and total distance traveled. Every photon must of course travel a distance of at least before reaching the target plane, and most photons were scattered several times.
Figure 6 shows the same distributions for type 2 scattering. Note than now 93% of the emitted photon packets eventually intersect the target plane. Only 7% of the emitted photons were wasted. Those photons ended up being scattered into directions away from the target plane (either by backscattering or by multiple large-angle forward scatterings). Figure 7 shows the results for type 3 tracing. The distributions are similar to those for type 2, but with over 94% of the photons reaching the target plane.
Figures 8-10 show the corresponding results when the target plane is at . Now, for type 1 tracing, only 179 of emitted photons ever reached the target plane. Over 98% of the photon-tracing calculations were wasted! Note also that there is obvious statistical noise in the distributions of the right panel, due to the small number of photons used to computed the statistics. For types 2 and 3 about 71% and 82%, respectively, of the emitted photons eventually reach the target plane. The statistical noise is now much smaller (but still noticeable) because of the larger number of photon packets.
The total optical distance distributions for types 2 and 3 show broad tails. In both cases, fewer than one fourth of the photons made it to the target plane after traveling a total distance of . About 30% of the photons underwent 30 or more scatterings and traveled a distance of . These broad tails illustrate the phenomenon of pulse stretching in time-dependent problems. If we think of all photons being emitted simultaneously, then the longer distances traveled correspond directly to later arrival times at the target. Pulse stretching is an important limiting factor in time-dependent applications such as lidar bathymetry or communications with high-frequency light pulses.
Figure 11 shows the distributions of numbers of photons at the target plane and the corresponding power (or energy) for the three tracing types and the detector at . The left panel shows the distributions as a function of the detector ring radii, and the right panel is the same information as a function of the bin number deﬁned above. Recall that the ﬁrst abscissa point is unscattered photons, the second is scattered photons inside the ﬁrst detector ring, and the last plotted point is for photons outside the last detector ring. The solid-line histogram represents the 10 detector rings for .
There are several things to notice in this ﬁgure. First, the distributions of the numbers of photons (open circles) are diﬀerent for the three tracing types. For type 1, the distribution of the number of photons is the same as the power distribution because the photons all retain their initial weight of . Thus power detected is simply the number of photons detected. For types 2 and 3, more photons are detected, but each is weighted less to account for absorption along the way. Finally—and most importantly—the power distributions for these three tracing types are identical to within a small amount of statistical noise, which is not visible in these plots.
Figure 12 shows the power distributions for the detector at but only photons emitted. There are obvious diﬀerences in the distributions for the three tracing types. However, if photons are traced, as in Fig. 13 these diﬀerences almost disappear. This indicates that the three ways to trace photons all give identical predictions of the detected power, to within some amount of statistical noise, which can be reduced by tracing more photons.
However, the computation time required by the three tracing types can vary greatly. Recall from Fig. 5 that about 30% of the photons reached the target plane at for type 1 tracing, but that over 90% reached the target plane for types 2 and 3 (Figs. 6 and 7). Thus, if we require a certain number of detected photons to achieve some desired level of statistical noise, we would have to emit and trace over three times as many photons for type 1 tracing (hence three times the computer time) as for type 2 or 3. For the detector at and type 1 tracing, fewer than 2% of the emitted photons reached the target plane (Fig. 8), whereas about 80% of the photons reached the target plane for types 2 and 3 (Figs. 9 and 10). Getting the same number of photons on target would thus require emitting over 40 times as many photons (hence 40 times the computation time) for type 1 as for types 2 or 3.
We have now shown that several ways of tracing photons can be devised and that each gives the same distribution of power or energy at a detector some distance away from the source. However, an intelligent choice of the photon tracing algorithm can greatly reduce the needed computations. Moreover, the computation diﬀerences depend on the particular problem, e.g., on detector distance from the source as shown here (or on the IOPs, not shown here). However, we can do still more to reduce computation times, which leads us to the next topic of variance reduction techniques.