1. Introduction

When light is incident on a turbid medium, which comprises of suspended particulate matter, it interacts in the form of scattering (and absorption). This light matter interaction is governed by the Maxwell’s equations, and depending on the problem, several theories are used to study a scattering system [1]. Thus, for example, for light scattering by spherical particles, analytical solutions can be obtained by using the Mie theory [2]. The propagation of light through a turbid medium can be described by several scattering parameters, such as the interaction coefficient $\mu _t$, albedo, and the anisotropy parameter $g$ [3,4]. These scattering parameters depend on the properties of the particles, like shape, size, concentration of particles, and also the refractive indices of the particles and the surroundings. The interaction coefficient dictates the probability of light interacting with a particle when propagating through the turbid medium; it can be expressed as the sum of the scattering coefficient $\mu _s$ and the absorption coefficient $\mu _a$, each of which is the respective probability of light being scattered and absorbed by a particle while propagating through a certain pathlength inside the medium. For a collection of identical particles, the interaction coefficient can be expressed as [2]:

The anisotropy parameter $g$ is the average of the cosine of the scattering angles from all the interaction events in the turbid medium, based on which the medium can be identified as primarily forward scattering ($g>0$) or backward scattering ($g<0$); $g=0$ corresponds to the case of isotropic scattering. The actual distribution of the scattering angle is described by the phase function $p(\theta )$, which is a measure of the amount of light scattered at an angle $\theta$ ($g$ being the first moment of the phase function) [5,6]. The phase function depends on the shape and size of the particles; some of the commonly used phase functions are the Henyey-Greenstein phase function [4], the Reynolds-McCormick phase function [7], and the Mie phase function [8]. Therefore, measurement of the scattered light can be used to determine the scattering parameters, and thus the particle properties. Along with the scattered light, measurements on the transmitted– and reflected lights can also provide information about the particles in a turbid medium [9]. Thus, light scattering techniques are ubiquitous when it comes to characterizing a turbid medium, with applications in probing biological tissues [3,10], analyzing the distribution of suspended particles [11–13], and monitoring quality of water in terms of suspended sediment concentration or turbidity [14–16].

The Monte Carlo method is widely used to simulate the propagation of light inside a scattering medium, and various simulation models can be found in the literature [17–20]. Such simulations are useful for verifying experimental results [21], or for optimizing the configuration of an experimental setup without repeatedly performing the experiment. But most importantly, Monte Carlo simulations can be used to estimate the scattering parameters by comparing the results of experiment and simulation [22,23]. In these kind of methodologies, the simulations are performed for different combinations of scattering parameters until the simulated values of scattered power match the experimental measurements, at which point the scattering parameters used for the calculations are concluded to be the scattering parameters of the given sample [4,7,8,24,25]. The simulations can also be used to pre-plot a calibration curve for a particular range of any scattering parameter, and thus results can be arrived at as soon as the experimental measurements are performed [26]. Since Monte Carlo simulations are considered to be "gold standard" [19], they can even be used in lieu of experiments to validate different analytical models [27,28]. Therefore, the Monte Carlo method has a large scope of application in various light scattering studies.

Most of the Monte Carlo models are applicable to turbid media which comprise of just one type of particles, i.e., the medium is described by a single phase function. In practical applications, however, there may be several types of particles present in the turbid medium, which can be considered as a mixture of several turbid solutions, each with a different value of $g$. The interaction coefficient of such a mixture can be expressed in different forms [29–33], all of which are effectively an extension of Eq. (1):

There are a few reports on the characterization of turbid mixtures using light scattering studies: Diehl et al. [35] carried out simultaneous measurements with several detectors to statistically identify and quantify broad categories of particles in mixtures. Green and Boon [36] described a methodology to identify concentrations of constituents in a mixture of sand and silt, using measurements on the backscattered light. Sharma [37] has suggested a method to determine the size distribution of particles in a fractal tissue model by using Monte Carlo simulations, and even mixtures comprising of different liquids have been characterized using light scattering [11,38].

The simplest case of a mixture is a binary turbid mixture, i.e., a turbid medium consisting of two types of particles, and characterization of such a medium is required in many industrial applications, such as determination of particle concentration [39] and study of flow behaviour [40,41] in a fluidized bed, or determining the effect of particle shape on segregation tendency [42]. In this paper, we present a method to determine the concentrations of constituents in such a binary turbid medium, by using Monte Carlo simulations for light scattering in mixtures. The simulation model is discussed in detail in Section 2, and it is shown that Eqs. (2) and (3) can be obtained by using a model wherein homogeneously distributed scattering particles are randomly moving inside a turbid medium. The methodology to determine the constituent concentrations from light scattering measurements, which also includes a recipe to prepare turbid mixtures and the implementation of the Monte Carlo model to simulate the experiment, is outlined in Section 3. The results are presented and discussed in Section 4, and Section 5 brings out the conclusions. An important aspect of our simulation model is that it accurately simulates the actual experimental setup, and therefore the scattered powers obtained from the simulations and the experiments can be directly compared without any additional calibration or modifications (see Section 3.3). For a binary medium, the methodology presented here achieves the same objective as that of Green and Boon [36]; however, their method requires experimental calibration with various samples, and also the relative concentration of the constituents need to be known beforehand. On the other hand, our method requires a single measurement of the scattered light and the (undeviated) transmitted light, to determine the respective concentrations of the two types of particles present in the turbid mixture.

2. Monte Carlo model for light scattering in turbid mixtures

In a typical Monte Carlo model [18,26], the propagation of light is represented in terms of the propagation of photon packets, where each packet comprises of a large number of photons. Every photon packet incident on the turbid medium undergoes scattering (and absorption, if applicable) in accordance with the scattering parameters of the medium, and a random variable $\xi$, using the inverse distribution method [3]. $\xi$ introduces the stochasticity in the model, which is required since the particles of the turbid medium move in a random manner. $\xi$ also allows the conversion of the angular distribution of the scattered light into the equivalent picture of a photon packet being scattered at an angle $\theta$, the probability of which is determined from a suitably chosen phase function $p(\theta )$. The azimuthal angle $\phi$ is considered to be equally distributed between $0$ to $2\pi$ for symmetric particles. After determining the value of $\theta$ and $\phi$ for a scattering event, the direction cosines $\mu _x$, $\mu _y$, and $\mu _z$ are determined, and the packet is propagated by a distance $s$ until the next point of scattering which is called the step-size [18]:

In this section, we present a model describing light propagation in a turbid mixture which may contain several types of particles (Fig. 1), and verify that the interaction coefficient and the phase function are indeed given by Eqs. (2) and (3) in the case of a mixture. Consider a turbid mixture in which there are $N_i$ number of particles of type $i$, each of which has a total interaction cross-section $C_{ti}$; $i=1, 2, 3$ represent the different types of particles. If the length of the medium is $L$, and the transverse area is $A$ (Fig. 1), then the concentration of each type of particles in the medium, $n_i$, can be written as

 

Fig. 1. a) A photon packet travelling a distance s₁ in a turbid mixture of volume V = LA, which comprises of multiple types of suspended particles. b) The distance $s_1=z_1/\mu _z$, divided into m equal segments of length Δz and volume ΔV . c) Shaded cylindrical volume in which a particle must be present to scatter the photon packet in that particular segment.

Download Full Size | PDF

Now, a photon packet having direction cosines $\mu _x$, $\mu _y$, and $\mu _z$ starts to propagate inside the medium, and say it travels a distance $s=s_1$ before getting scattered. The probability of the photon packet travelling a distance $s_1=z_1/\mu _z$, without being scattered by a particle, is designated as $P(s\geq s_1)$ [18]. To determine this probability, we divide the turbid medium into a large number of equal segments of volume $\Delta V$ and length $\Delta z$, so that $z_1=m\Delta z$, where $m$ is the number of volume segments, and $\Delta V=A\Delta z$ (Fig. 1(b)). The particles are assumed to be stationary inside each segment (in comparison with the velocity of the photon packet), and placed at random coordinates. If only one particle of the $i^{th}$ type having cross-section $C_{ti}$ is present in the first segment (Fig. 1(c)), then the photon packet will be scattered whenever the particle lies within the shaded tilted cylindrical region with cross-sectional area ($C_{ti}/\mu _z$) and length $\Delta s=\Delta z/\mu _z$. This can be geometrically shown by treating the particle as a sphere of radius $\sqrt {(C_{ti}/\pi )}$. Since the particle is positioned randomly, the probability that the particle lies in the shaded cylindrical volume, and hence scatters the photon packet, is equal to the volume of the shaded region divided by the total volume $\Delta V$ of the segment:

From Eqs. (2) and (3) (or equivalently, Eq. (13)), it can be inferred that the step-size of the photon packet, and the effective phase function of the turbid mixture, both depend only on the products of the concentration and the cross-section of the different types of particles. In our Monte Carlo model, instead of using the effective phase function (Eq. (3)), the scattering angle after each interaction is chosen randomly from any one of the different $p_i(\theta )$, depending on the concentrations $n_i$ of the particles (Eq. (13)) [44] (also, see Section 3.3). Although this method requires the choice of an additional random number after each scattering event, it offers more versatility in the following manner: For example, if some types of particles are described by the Mie phase function, which requires lookup table-based sampling [43,45], and some other types of particles are described by analytically invertible phase functions (such as the Rayleigh phase function or the Henyey-Greenstein phase function [3]), then the analytical functions do not have to be converted into a discretized lookup table (which would be necessary in case an equivalent phase function is used). Further, each lookup table (corresponding to the different phase functions) can have a different number of elements in our case, whereas in the other method, the equivalent phase function has to be constructed with the minimum number of elements among all the lookup tables.

In the case of particles with finite absorption, the weight of the photon packet is reduced after each scattering event, which can be calculated from the cross-sections as [2]

3. Methodology

3.1 Experimental setup

The experimental setup is shown in Fig. 2. Light from a He-Ne laser (REO make, model number 30991, wavelength 632.8 nm, power 5 mW) is incident on the turbid sample contained in a quartz cuvette. The scattered light is measured at an off-axis detector (Orion-7Z01803 by Ophir), $D_1$, which has a circular aperture of diameter 2.2 mm, and the transmitted light is measured at a detector, $D_2$, placed far ($\sim$1 m) from the sample so that (almost) no scattered light reaches the detector. In our experiment, the distance between the cuvette and the detector $D_1$ is $l_x=9.2$ cm, while the lateral shift from the direction of the incident beam is $l_y=1.4$ cm. The position of the detector $D_1$ was chosen such that the detector is close to the axis of the incident beam (without obstructing the transmitted light), which ensured that the collected scattered power would be high enough, yet it is sufficiently far from the cuvette that light scattered from the cuvette walls do not reach the detector. In principle, however, the detector can be placed at any arbitrary position, depending on the sophistication of the setup. Measurement of the transmitted light allows us to determine the total interaction coefficient of the sample by using Beer-Lambert’s law [4,11]:

The measurements were performed in a dark room. To remove the effects of scattering due to imperfections in the cuvette walls and loss due to Fresnel reflections the measurements on the transmitted light were normalized with respect to distilled water. The measurements on the scattered power were offset to remove the background noise, which is the power collected at the detector when the cuvette contains only distilled water. The measurements were carried out several times for every sample to check repeatability of the scattered and transmitted powers at detectors $D_1$ and $D_2$, respectively.

 

Fig. 2. Experimental setup to measure the scattered power at an off-axis detector, $D_1$, and to estimate $\mu _t$ by measuring the undeviated transmitted light at detector $D_2$.

Download Full Size | PDF

As discussed in Section 1, the scattered power depends on the scattering parameters of the turbid medium, like $\mu _t$, $\mu _a$, and $g$, which depend on the physical properties of the scattering particles. Since $\mu _t$ depends on the cross-section of the particles, it can be said that the (undeviated) transmitted power also depends on the size, shape and the concentration of the particles. Therefore, the scattered– and transmitted powers contain information about the turbid medium, which can be extracted; in our case we determine the $\mu _t$ and the constituent concentrations, as explained in the following sections.

3.2 Preparation of samples

To experimentally validate our model, it was necessary to prepare turbid mixtures in which the concentrations of the constituent particles were known. Consider two monodisperse solutions which are mixed to form a bidisperse mixture; let the respective volumes of the first and second solutions be $V_1$ and $V_2$. The concentration of particles in the first solution is $n_{m1}$, and the cross-section of each particle is $C_{t1}$; similarly, the concentration of particles in the second solution is $n_{m2}$, and the cross-section of each particle is $C_{t2}$. Therefore, the interaction coefficient of the first solution is $\mu _{t1}=n_{m1} C_{t1}$ and that of the second solution is $\mu _{t2}=n_{m2} C_{t2}$. If $n_1$ and $n_2$ are the respective concentrations of the two types of particles in the mixture, then we get the probability of scattering by the first particle, $p_{1s}$, as (see Eq. (13))

For our experiment, monodisperse solutions of different concentrations were prepared by diluting stock solutions of polystyrene microspheres of sizes 1 µm, 2 µm, and 3 µm. The stock solutions for the 1 µm and 2 µm sizes were procured from Alfa Aesar, while the 3 µm sized stock solution was procured from Sigma-Aldrich. By mixing two solutions of different particle sizes, bidisperse mixtures with different particle concentrations could be obtained.

3.3 Simulation of the experimental system

The scattering system is simulated using the Monte Carlo model described in Section 2, with the same values of $l_x$, $l_y$, and other geometrical parameters as that in the experiment (Section 3.1). In the simulations, 100 million photon packets are launched into the turbid medium, and the step-size between each scattering site is obtained from Eq. (4). Since polystyrene is nonabsorbing at He-Ne laser wavelength [46], the weight of the photon packet is not reduced after each interaction. A basic flowchart showing the various steps involved in the simulation is shown in Fig. 3. When a photon packet is scattered, the phase function corresponding to the first particle, $p_1(\theta )$, is used to calculate $\theta$ if the random number $\xi <p_{1s}$, otherwise the phase function corresponding to the second particle, $p_2(\theta )$, is used. Since the particles in our turbid solutions are spherical, the Mie-phase function is used for the simulations [26,32,43,45]. The phase function and the cross-sections of the polystyrene spheres (refractive index 1.587 at 632.8 nm [47]) of different sizes are calculated from Mie theory, using the MiePlot software by Philip Laven [48]. The azimuthal angle after each interaction is $\phi =2\pi \xi$ since the particles are spherical.

 

Fig. 3. Flowchart for Monte Carlo simulation of the experimental setup.

Download Full Size | PDF

If a photon packet exits the cuvette towards the direction of the detector, it is propagated till the plane of the detector and its weight is collected if it enters the detector area. If the packet does not reach the detector area, or if it exits from any other surface, then its propagation is terminated. Once a photon packet reaches the detector area, or is terminated otherwise, the next packet is launched. After all the packets are launched and propagated, the total weight collected at the detector divided by the total number of photon packets launched corresponds to the experimentally obtained fractional scattered power.

We may mention that when a photon packet exits the cuvette, it is refracted and the weight is reduced by using the Fresnel refraction coefficients [49]; for accurate results, the lateral shift and the possibility of total internal reflections at the cuvette walls (thickness 1 mm) is also taken into account. Thus, the scattered powers obtained experimentally and from the simulation are essentially equivalent to each other. For determining the refraction coefficient, the refractive index of quartz is calculated to be 1.457 from the Sellmeier equation [50], and the refractive index of water is taken to be 1.3317 [51].

3.4 Determination of constituent concentrations

Given a bidisperse mixture with known constituents, i.e., the cross-sections $C_{t1}$ and $C_{t2}$ are known, the experimental setup described in Section 3.1 is used to measure the scattered power, and also to determine the interaction coefficient $\mu _t$ (using Eq. (15)). The experiment is simulated for this value of $\mu _t$, and the scattered power is obtained for different values of $p_{1s}$ ranging from 0 to 1 (Section 3.3). Thus, a reference plot of the fractional scattered power (at the detector) vs. $p_{1s}$ is obtained, which has a monotonic variation, increasing or decreasing depending on the scattering probabilities of the two types of particles. Corresponding to the measured power for the sample mixture, the numerical value of $p_{1s}$ can be determined from the reference plot. Since the values of $\mu _t$ and $p_{1s}$ of the given sample are now known, the respective concentrations of the constituents, $n_1$ and $n_2$, can be determined from Eq. (16):

4. Results and discussion

To verify the applicability of the proposed method, bidisperse turbid mixtures of polystyrene spheres with different values of $p_{1s}$ were prepared (Section 3.2). By measuring the scattered power, the value of $p_{1s}$ for each of the samples was estimated from the corresponding reference plot (Section 3.4), which was obtained by changing $p_{1s}$ in increments of 0.1, and thereafter the constituent concentrations were determined using Eq. (18).

The reference plots corresponding to the simulation results for two different mixtures A and B, one comprising of polystyrene spheres of sizes 2 µm and 3 µm, and the other with spheres of sizes 1 µm and 3 µm, are shown in Fig. 4. For a given pair of particle sizes, the simulated scattered powers and hence the slope of the curve depends on $\mu _t$ and the position of the detector (which is fixed in our setup). It can be seen that for a fixed value of $\mu _t$, the scattered power collected at the detector varies monotonically with $p_{1s}$; the values of $p_{1s}$ of the two samples can be estimated form the reference plot without any ambiguity, as shown in Fig. 4. It should be noted that although the plots appear to be linear, the theoretical analysis does not necessitate the scattered power to vary linearly with $p_{1s}$, and depending on the concentrations and the cross-sections of the particles, the variation can be non-linear as well; in this case, a smaller step size of $p_{1s}$ may be required for improving the accuracy. In the figure, errors bars are not shown because the standard deviation of the measurements were found to be negligible compared to the value of the scattered power. Since $p_{1s}$ of the prepared samples are known, the actual values of $n_1$ and $n_2$ of the samples can also be calculated using Eq. (18). The results for some different samples are shown in Table 1, where the actual values of the constituent concentrations are compared with those obtained by using our method.

 

Fig. 4. Reference plots of the fractional scattered power vs $p_{1s}$ for two bidisperse mixtures A and B (comprising of polystyrene spheres), obtained from Monte Carlo simulations. It can be seen that the variation is monotonic, i.e., there is one-to-one correspondence between $p_{1s}$ and the scattered power. The estimated value of $p_{1s}$ corresponding to the measured scattered power is shown by dotted line for both the mixtures.

Download Full Size | PDF

Table 1. Estimated values of $p_{1s}$, along with the calculated constituent concentrations, $n_1$ and $n_2$, for different bidisperse turbid mixtures. The actual values of $p_{1s}$, $n_1$, and $n_2$ of the prepared samples are given in parentheses, with respect to which the relative errors are calculated. The parameters corresponding to the turbid mixtures A (2 µm, 3 µm) and B (1 µm, 3 µm) are shown in bold.

View Table

The simulations were performed on a personal computer (operating system: Microsoft Windows 10 Home, central processing unit: Intel i5-4210U @ 1.70 GHz, memory: 8 GB RAM), and plotting each curve took an average of 40 seconds per million photon packets. The exact time taken depends on the magnitude of the scattering parameters being used in the simulation; for example, increase in the interaction coefficient will result in an increase in the simulation time, because a higher value of $\mu _t$ leads to a larger number of scattering events, each of which involves obtaining a step-size, choosing a phase function, and calculation of the direction cosines. As each iteration of launching a photon packet and propagating it are independent, the calculation speed can be drastically increased by using more powerful computation techniques such as parallel processing, if a faster result is required.

Table 1 also shows the relative error in the estimation of the concentrations of the constituents, which is primarily due to two reasons: Firstly, there may be a small difference in the scattered power obtained from the simulations and the experiments, because of different practical factors. For example, the cuvette may not be exactly perpendicular to the incident light, the measurement of detector position may have slight inaccuracies, the circular aperture of the detector may not be perfect, and so on. The second source of error is not due to any limitation of the method, but due to error in calculation of the actual concentrations of the constituents in a prepared sample. This error is introduced during the estimation of the $\mu _t$ of the monodisperse solutions before mixing, which involves measurement of the transmitted light. A possible source of error may arise due to the fact that the experimental samples did not comprise of particles of only two sizes: there is a very narrow but finite size distribution in the procured polystyrene sphere solutions. However, the simulations were performed for strictly bidisperse samples. If more accuracy is desired, the effect of the size distribution may be included by integrating the phase function over the distribution [31]. This was not done in the present work because the exact size distributions of the polystyrene samples were not known, which is usually the case in practical applications. A small error may also be introduced during mixing of the monodisperse solutions, since measuring exactly equal volumes of the solutions may not be practically possible; in our case, this error was negligible because the volume is measured using micropipettes which have a high accuracy. Some other sources of error, such as fluctuations in the power of the incident light, and local inhomogeneities in the prepared mixture, may also be present.

Although the method is validated for mixtures of particles of the same material (polystyrene), but of different sizes, this is equivalent to a general binary mixture, because the cross-sections and phase functions are different, and hence the particles can be considered to be of different types without any restrictions or approximations. Further, the method can be applied to absorbing particles as well, provided the absorption cross-sections $C_{a1}$ and $C_{a2}$ of the particles are known, in addition to the cross-sections $C_{t1}$ and $C_{t2}$. In such cases, the simulations will include a reduction of the weight after each interaction, which is given by Eq. (14).

The Monte Carlo model described in Section 2 is applicable in the case of independent (incoherent) scattering [2,52], which suffices for our method since it involves static light scattering measurements. However, in the case of coherent light scattering [53], which is important in dynamic measurements [54], the model would require suitable modifications. Another aspect of this study is the possibility of characterizing a turbid mixture having three or more constituents, which may involve measurements at different configurations using multiple detectors. In this case various other factors need to be introduced, such as possible modifications to the Monte Carlo method, and the uniqueness of the solution; efforts in this direction are presently underway.

5. Conclusion

A Monte Carlo model for light propagation through a turbid mixture is presented, and the model is used to formulate a recipe to determine the concentrations of constituents in a turbid mixture with two types of particles. The method is experimentally validated by using in-house prepared bidisperse mixtures of polystyrene spheres. The calculated concentrations of the constituents are found to be in close agreement with the actual constituent concentrations of the samples. The proposed method requires a single measurement of the scattered and the transmitted powers to determine the constituent concentrations, provided the cross-sections of the constituent particles are known.