Also available in 日本語
This article describes how to simulate polarization-sensitive bulk scattering and fluorescence using a custom DLL in OpticStudio. The bulk scattering model defined in MSP.dll (available for download with this article) considers the polarization of incident non-sequential rays and simulates how the polarization and direction of propagation change with each scattering event. This DLL can also be used to simulate fluorescence in combination with Mie scattering. Both fluorescence and polarization-sensitive scattering are important for modelling biological imaging. This article summarizes seven experiments that use the MSP DLL bulk scatter model.
Authored By Guillem Carles, Research Associate at University of Glasgow
Introduction
The modeling of light diffusion through turbid media is important for a range of fields, including remote sensing, underwater imaging, and imaging through atmospheric turbulence. It is particularly important for biological imaging, because the system design requires the modeling of light propagation in turbid media and scattering samples. Using Non-Sequential Mode in OpticStudio is very convenient for such modeling, provided accurate models of the imaging system are known, and valid and rigorous models for volumetric scattering are available [1]. This article reviews the use of the MSP DLL, a custom DLL that can accurately model in OpticStudio the polarization of rays that undergo bulk scattering with the Mie phase model. For more details on its implementation, see [1].
Modeling scattering
When a beam of light enters a turbid medium, it suffers extinction mainly due to two factors: absorption and scattering [2]. To model this, a Monte Carlo simulation can be used to sample the beam using rays, and then the rays can absorb and scatter as they propagate through the turbid medium. Each ray is absorbed or scattered according to defined probabilities that depend on the optical parameters of the medium.
One option for modelling absorption, which is used in OpticStudio, is to decrease the intensity of the rays according to the Beer-Lambert law as exp(-zµa), where z is the distance travelled inside the medium and µa is the absorption coefficient. The following figure illustrates this. On the left is a macroscopic view of the beam; on the right is a Monte Carlo simulation using traced rays:
The beam is sampled by N rays (in this case, the rays are randomly distributed in the cross-section of the beam), and each ray is assigned an intensity Iray that is equivalent to the intensity of the beam I0 divided by N. Each ray undergoes absorption as it travels through the medium. The extinction rate that is due to absorption is defined by exp(-µaz).
When scattering occurs, the intensity of each ray is preserved, but the direction of propagation is changed during a scattering event. As a result, rays are removed from the beam, also causing extinction. In this case, the scattering coefficient µs (the inverse of the mean-free path, which determines the probability of the occurrence of scattering events) determines the scattering extinction rate at exp(-µsz). This can be illustrated as follows:
The properties of the medium (such as particle size, index of refraction, and particle density) determine the mean-free path, the angular dependence of the scattering, and the polarization. The resulting scatter model changes the direction of light propagation, which means that light is propagated to other directions than the direction of the incident beam. The properties of the scattering medium produce a far-field intensity distribution, but the scattered rays preserve their original intensity if there is no absorption. (Note: For this model to be valid, the scattering events must be mutually independent.) Therefore, beam extinction is the consequence of the reduced number of rays that preserve the incident direction at the output.
On a scattering event, the new direction is determined using a probability distribution function that is set to match the far-field intensity distribution. For historic reasons, it is called the phase function. This means that tracing enough rays reproduces the far-field distribution. Examples of common phase functions are the Henyey-Greenstein phase function, the Mie phase function, and the Rayleigh phase function.
Polarization-sensitive scattering
In general, phase functions are dependent on the polarization of incident rays [2]. The practical and widely used Henyey-Greenstein phase function has very convenient properties for simulating scattering in biological samples, but it does not account for polarization. Some scattering simulations require a more detailed physical model, such as Mie scattering, as well as the consideration of polarization. The MSP DLL can be used to model Mie scattering with the correct tracking of polarization (that is, the correct electric field associated with each ray).
To use the MSP DLL in a non-sequential OpticStudio file, save the DLL locally to the BulkScatter user data folder <…\Documents\Zemax\DLL\BulkScatter>, and then access it in the Object Properties of a non-sequential volume (Object Properties > Volume Physics tab > DLL drop-down menu > MSP.dll). Note that the parameter Index refraction (real) should have a value greater than zero, its default value.
Plots of the phase function are known as scattering diagrams. The following diagrams show examples of Mie scattering for different particle sizes (parameter α) and different polarization states (parameter L; L=0 unpolarized and L=1 linearly polarized). In the following diagrams, a ray is initially propagating along the z axis. Its polarization (whether linear or elliptical) is oriented towards the y axis. The surface denotes the probability of scattering in each direction. Notice how the phase function and therefore the direction of the ray’s propagation after a scattering event are dependent on the polarization of the incident ray, L, just before it scatters.
Scattering events that occur based on the Mie model defined above will also modify the polarization of the ray deterministically, so it is important to update the electric field of the ray as it is being scattered.
Correct treatment and tracking of polarization at a scattering event requires three steps:
- Calculate the degree of linear polarization of the incident ray, L
- Calculate the new direction of propagation using the polarization-sensitive phase function
- Update the electric field (polarization) of the ray after the scattering event
In OpticStudio, the other bulk scattering models that do not consider polarization assume a random polarization for each scattered ray. The electric field is correctly perpendicular to the new direction of propagation that satisfies k*E=0, but the actual polarization is lost.
Experiments
We performed several simulations in OpticStudio using the MSP DLL to explore the scattering dependency on polarization.
Experiment 1
In the first experiment, we configured a setup such that rays from a point source are launched toward a scattering sample, and a detector is placed above the sample, as shown in the following image.
We set the source to allow only one scattering event (in OpticStudio > on the Properties tab of the source), and we set the mean-free path to a very small number (mfp=0.001mm). This configuration ensures that rays scatter only once and do so very closely to the surface of the sample. The sample material has a refractive index of 1, so its shape and size do not affect the experiment. Therefore, the intensity distribution that is recorded at the detector directly depends on the phase function of the scattering properties of the sample material.
The particle size is set to a very small number, such that scattering is in the Rayleigh regime. In this case, a simple closed-form solution is known for Rayleigh scattering [2], and the results can be readily compared with the theory. The following diagram shows the intensity at the detector for incident unpolarized or circularly polarized rays (the incoherent irradiance across the detector is on the left; the corresponding profile from the central column is on the right).
Notice how the vertical direction of the detector, which is labeled Y coordinate value in the graph on the right, corresponds with the propagation direction of the incident rays because the detector is placed parallel to the incoming rays (see the experimental layout above). With this geometry, each pixel of the detector subtends a different solid angle and the signal is reduced toward the edges, as expected.
If we compute the theoretical Rayleigh phase function for unpolarized light, we observe a perfect match:
The left graph above is the theoretical phase function for Rayleigh scattering of unpolarized light. The red-colored portion represents the range of angles that the detector spans in this setup. These angles are therefore captured by the detector (detector length is 1mm and the sample is placed at 0.2mm from the detector). Results projected on the detector are shown in the graph above on the right.
Experiment 2
Setting the polarization of the incident rays to be linear and vertical (in OpticStudio, on the Source tab of the object properties of a source object) produces entirely different results:
As shown in the lower-right graph above, the the MSP DLL correctly simulates Mie scattering and polarization to produce results that match the expected theoretical results.
Experiment 3
The previous experiment highlighted the dependency of the phase function on polarization. We can perform a further experiment to highlight that scattering alters polarization: Consider the same setup as above, with unpolarized incident light, but now we perform polarized detection to only detect light polarized at the detector’s local Y direction. (In OpticStudio > Detector Rectangle, set Polarization parameter to 2.) The following results again show a perfect agreement with the theoretical expectation:
As shown in the lower-right graph above, the the MSP DLL correctly simulates Mie scattering and polarization to produce results that match the expected theoretical results.
Experiment 4
Another way to explore the dependency of the phase function for various states of polarization is to perform direct measurements in OpticStudio using single-pixel detectors placed at a range of scattering angles, and to plot the intensity of each detector:
In the graph on the right above, each detector represents a measurement point. The expected theoretical results for Rayleigh scattering are also plotted for reference. In the graph, blue, red, and yellow correspond to linear polarization in the x direction, linear polarization in the y direction, and unpolarized, respectively.
Experiment 5
So far, we’ve explored the polarization dependency in the Rayleigh regime. The MSP DLL is also valid for general Mie scattering. In this experiment, we use the same setup as in the first three experiments but with larger particles so the scattering is in the Mie regime. We also compare the results against theory (although Mie calculations are not in a simple closed form, they can be numerically computed. For example see [3]). Results are shown in the following diagrams for particles sizes of (left) 0.1λ, (center) 0.2λ, and (right) 0.5λ:
As expected, increasing the size of the scattering particles, forward scattering increasingly dominates, which is reflected in the detected intensity distributions.
Experiment 6
Similarly, we can measure the phase function using a setup with a range of single-pixel detectors that are placed 360° around the sample:
The polar plots on the right were constructed using the intensity values that were measured by each detector. These plots resemble perfectly the known scattering diagrams for general Mie scattering.
Experiment 7
Finally, we perform a more generic polarimetric experiment. We calculate the Mueller matrices (one is associated with each detector pixel) as spatial maps by measuring backscattering coming from a slab of turbid material. A beam of light is launched toward the slab of turbid material. Upon entering the scattering material, light rays scatter several times. Some rays eventually propagate back out of the object and are detected. A detector is placed close to the face of the slab and only detects back-scattered rays:
Using the MSP DLL, you can calculate the polarization-sensitive scattering and relate the polarization of the output rays to the polarization of the input rays. The Stokes vector, a standard way to describe polarization, uses four components to define any general state of polarization. Thus, the effect of any optical component on the polarization of a light beam can be completely defined by the Mueller matrix, which relates the input and output polarization states expressed as a four-component vector. It is therefore a 4x4 matrix:
For this experiment, the slab of turbid material can be considered an optical component that modifies the polarization state of the backscattered rays. The location of the rays on the detector is relative to the point of injection, which is the center. Therefore, it is possible to calculate its associated Mueller matrix. However, this 4x4 matrix depends on the point of exit of the backscattered rays, which is why the Mueller matrix is calculated as a spatial map. Results are plotted as follows:
These results can be considered to be a general case for polarimetric sensing, because they represent a measurement of full polarization of (in this case) the backscattering from a slab of turbid material.
Simulation of fluorescence
In addition to calculating polarization-sensitive Mie scattering, the MSP DLL can also incorporate the simulation of fluorescence.
To simulate fluorescence, the MSP DLL uses the wavelength shift capability that OpticStudio offers for bulk scattering. (In OpticStudio, the Object Properties of the scattering object > Volume Physics tab > Fluorescence box.) By entering text strings “wave1, wave2, prob” it is possible to set the probability such that on a scattering event, the incident wavelength number “wave1” is shifted to wavelength “wave2” with probability “prob.”
This approach makes fluorescence occur at a rate set by the scattering mean-free path (in OpticStudio > Object Properties tab > Volume Physics > DLL Defined Scattering > Mean Path). If a longer mean-free path is required for the fluorescence (the mean path length a ray travels before it fluoresces), it can be done effectively by simply reducing the probability of the wavelength shift, as described above. Using the MSP DLL, it is also possible to use shorter mean-free paths for fluorescence than for scattering. You can do this by entering a fluorescence-associated mean-free path, which determines the mean distance a ray propagates before it fluoresces (with probability as defined above). This distance can be set lower than the scattering mean-free path, but not higher.
If fluorescence occurs, the ray undergoes a wavelength shift, a change in direction of propagation, and polarization. The direction of propagation is selected isotropically, the polarization is set linear at random direction, and the phase is randomized. These properties simulate incoherent fluorescence emission after they are averaged over many rays. This happens only if fluorescence occurs.
If the ray is bulk scattered, the scattering mean-free path applies, the new direction of propagation is calculated based on the material properties and the phase function for Mie scattering, and polarization is updated according to Mie theory. Scattering and fluorescence can therefore be simulated simultaneously. This can be used to simulate both samples and imaging through a scattering material that expresses auto-fluorescence [1].
References
- G. Carles, P. Zammit, and A. R. Harvey, “Holistic modeling of optical systems and photon transport in turbid media using a commercial ray-tracer,” in preparation.
- H. C. van de Hulst, Light scattering by small particles (John Wiley & Sons, 1957).
- W. J. Wiscombe, "Improved Mie scattering algorithms," Appl. Opt. 19, 1505-1509 (1980).