## Abstract

This study is devoted to the development and application of a Monte Carlo ray-tracing model to simulate light scattering when a colloid suspension droplet passes through a highly focused Gaussian laser sheet. Within this study, a colloidal suspension droplet refers to a spherical droplet containing multiple spherical inclusions. Such scattering scenarios arise when using the time-shift measurement technique for particle sizing. The incident laser sheet is treated as a large number of polarized light rays: the Stokes vector of each light ray is tracked, achieved by multiplication of the rotation matrix and the Mueller matrix after each scattering event. For the Monte Carlo simulation of light scattering, a very important issue is to generate the deflection angle and azimuthal angle after each scattering event. The scattering from embedded inclusions is computed using the Lorenz-Mie theory and by employing the rejection sampling technique to update the new propagation direction. Multi-reflection and refraction within the droplet is accounted for, as is total reflection at the drop interface. For this, the Mueller matrix formulation is invoked at the drop surface to update the Stokes vector. To validate this simulation code, the scattering diagram from a nanoparticle is computed with this Monte Carlo method and compared with the scattering diagram computed with the Lorenz-Mie theory, the agreement is excellent. This Monte Carlo code is then applied to simulate signals arising from a time-shift device, when a colloid suspension droplet passes through a focused polarized laser sheet, with the objective of measuring the concentration of colloidal particles in the droplet. Measurements verify the ability of the code to properly simulate this light scattering scenario.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Colloid suspension droplets are widely encountered in industrial processes, for instance fuels with dispersed, metallic nanoparticles, metallic paints with embedded aluminum flakes, colloid suspension droplets (e.g. ink-jet printing) or encapsulation of particles, in which a liquid layer covers the particle, typical in pharmaceutical applications. It is of interest to characterize such droplets using not only the drop size, but also the colloidal concentration. Several existing optical measurement techniques can measure the size of a colloid droplet [1,2]; however, few techniques have been introduced to measure the colloidal concentration [1,3,4,5].

The first study for the estimation of the colloidal concentration using the time-shift technique was made by Rosenkranz et al. [2], in which a two-dimensional Monte Carlo model for the light ray propagation in turbid media was described. The current study is an extension of this work into a fully three-dimensional ray tracing model for arbitrary phase matrices, similar in capability to the model published by Stegmann et al. [6]. In the work of Rosenkranz et al., the simulation was limited to a two-dimensional plane and the polarization of the laser beam was not considered. Therefore, for the scattering between the light ray and the colloidal particles, the Henyey-Greenstein (HG) phase function could be used to generate the new deflection angle for updating the propagation direction of the light ray. However, in practice the time-shift instrument uses a polarized laser beam.

In the present study, the model is therefore extended to three dimensions and the polarization effect of the light ray is considered in order to simulate the time-shift signal more precisely. In this case, instead of the HG phase function for updating the new propagation direction, rejection sampling is used to generate the new deflection angle and azimuthal angle after the scattering event with internal colloidal particles [7–13]. For the scattering event between the light ray and the drop surface, the Mueller matrix is used to calculate the Stokes vector for the reflected light ray and transmitted light ray: a detailed description of the Mueller matrix formalism can be found in [8,14]. However, when total internal reflection occurs, because of the phase shift of the electric field, the Mueller matrix elements *S _{34}* and

*S*can take values other than 0; therefore, these formulations are not suitable to calculate the Stokes vector. In section 2.2 therefore, the Mueller matrix is derived with consideration of total internal reflection.

_{43}This polarized Monte Carlo method is then used to examine signal generation in a time-shift device as a function of colloidal concentration. The goal is to identify signal characteristics which may be used to estimate the colloidal concentration from individual droplets, in addition to their size and velocity.

Before embarking on the description of the simulation, some remarks are appropriate about assumptions and limitations of the approach taken. One main assumption is that the droplets are exactly spherical. Although in practice droplets are typically deformed and become spheroidal, this assumption is not too restrictive, since in previous work the time-shift technique has been shown to be rather accurate for droplets with oblate aspect ratios in the range 0.9 -1.1 [15]. Direct imaging of the droplets studied in the present study indicate that the oblateness is much lower, typically 0.99-1.01. A second assumption is that the colloid inclusions are uniform in size and spherical. While this also is not likely representative of typical industrial situations, it is the logical first step to studying the problem and future work will investigate the influence of non-uniform size distributions. Non-spherical inclusion particles would be much more difficult to treat; however, some first work about the response of the time-shift technique to plate-like inclusions, typical of metallic paints, has been presented by the authors [16]. Finally, we note that only dilute colloidal suspensions are studied, up to a volume concentration of 3%. This limitation allows the use of the Beer-Lambert law for extinction of individual light rays, as will be addressed in section 4.

## 2. Code description

#### 2.1 Geometric model of the colloidal droplet

The spherical droplet considered has a diameter of *D* and the colloidal particles have the uniform diameter of *d*, as pictured in Fig. 1. The volume of the droplet and the particles are denoted *V _{D}* and

*V*respectively. The extinction cross section

_{P}*σ*of the particle is

_{e}*Q*

_{ext}_{·}

*πd*, in which

^{2}/4*Q*is the extinction efficiency factor. The refractive indices of the surrounding medium, droplet and the colloid particles are

_{ext}*n*,

_{1}*n*and

_{2}*n*respectively. The colloidal particles are randomly distributed in space in the droplet.

_{3}The volume concentration *C _{v}* of the particles is specified using an optical mean path length

*E*(

*L*). The relationship between the volume concentration and the optical mean free path length is given by [2]

#### 2.2 Light ray intersects with drop interface

The laser sheet is treated as a very large number of light rays, as in [16,17]. The polarization state of individual light rays is specified with the Stokes vector *S*, given in Eq. (2) [18]. The propagation direction of the light rays is also specified, which can account for Gaussian beam properties, e.g. divergence. Nevertheless, since the measurement volume of the time-shift instrument is in the Rayleigh range of the Gaussian beam, the propagation direction of all light rays has been kept constant in the present study.

The Stokes vector is tracked during its propagation by each scattering event. After a scattering event, the light ray is scattered from the propagation direction *k* to a new direction *k’*, as depicted in Fig. 2. The scattered Stokes vector is updated by multiplication of the rotation matrix *R*(*φ*) and the Mueller matrix *M*(*θ*), according to Eq. (3).

The rotation matrix *R*(*φ*) transforms the polarization state to a new incident reference frame and is written as

When the light ray intersects with the drop surface, the amplitude of electric field of the reflected ray and transmitted ray is calculated using the Fresnel equation and Snell’s law. The computational procedure can follow either a vector decomposition, as in [16,23], or by multiplying the Stokes vector with the Fresnel equations, which should be in the form of Mueller-Stokes matrix. The Mueller-Stokes matrix at the dielectric surface is well described in the literature; however, it does not consider the case of total internal reflection. Here the Mueller matrix at the dielectric surface is derived with consideration of the total internal reflection, to update the Stokes vector.

By treating the incident electric field and the scattered electric field separately, the scattering is described in matrix form [18]

*R*is the radial distance from the scattering particle;

*S*is called the complex scattering matrix. It contains the scattering properties of the particle and can be written as

*A*and

_{11}*A*are zero [24]. For reflection at the drop surface, the reflected electric field for the parallel and perpendicular polarization is written as

_{12}*r*and

_{p}*r*are the Fresnel coefficients for reflection. In a similar manner, for the transmission, the transmitted electric field is given by

_{s}*t*and

_{p}*t*are the Fresnel coefficients for transmission. When total internal reflection occurs, because the phase shift is no longer 0 or

_{s}*π*, the elements of the

*S*can have imaginary parts [25]. Instead of using the complex scattering matrix

*S*, the scattering process is described with the real Mueller matrix. In the case of spherical droplets, only 6 elements in the Mueller matrix are independent and can be written as

*R*and

*T*are the Mueller matrixes for the reflection and transmission at the droplet surface respectively. For the case of light ray transport from optically thin medium to an optically denser medium, the Fresnel coefficients are real and given by [27]:

*S*and

_{34}*S*are zero. For the case of light ray transport from optical denser medium to optical rarer medium, the Fresnel coefficients are given by Eqs. (22)–(23). When the incident angle

_{43}*θ*is larger than the critical angle

_{i}*θ*, total internal reflection occurs. Because of the phase shift of the reflected and transmitted electric field [28,29], which can take values other than 0 or π, the Fresnel coefficients will be complex [30]:

_{c}*i*is the imaginary number,

*m = n*, and

_{t}/n_{i}*n*and

_{t}*n*are the refractive index of the optical denser medium and the optical rarer medium respectively. When the incident angle

_{i}*θ*is smaller than the critical angle

_{i}*θ*,

_{c}*r*and

_{p}*r*are still real numbers. Hence, the element

_{s}*S*and

_{34}*S*are zero, however, when the incident angle is larger than the critical angle, total internal reflection occurs and

_{43}*r*and

_{p}*r*become complex numbers and the element

_{s}*S*and

_{34}*S*are no longer zero.

_{43}#### 2.3 Light ray intersects with suspended colloidal particles

For the scattering within the droplet, a random walk model is utilized to predict the light ray path [2]. The procedure used follows that described by the flow chart given in [7]. The essential task is to generate the polar deflection angle *θ* and azimuthal angle *φ*, to yield the new propagation direction of the light ray. The scattered intensity by the isotropic sphere is a function of *θ*, *φ* and the polarization state. Using Eq. (3), the scattered intensity can be expressed by Eq. (24)

*θ*and

*φ*; for the non-polarized beam, whose Stokes vector is [1 0 0 0]

^{T}, the scattered intensity is only a function of

*θ.*For this case, the Henyey-Greenstein phase function can be used to predict the angle

*θ*, and the sampling of angle

*φ*can be treated as isotropic.

Figure 3 shows the three-dimensional scattered intensity from a colloidal particle, which is illuminated by plane waves exhibiting vertical linear polarization (s) and right-hand circular polarization respectively. The intensity is plotted with color in logarithmic scale. The incident wave impinges on the particle in the positive z direction. As the figure shows, the polarization state has a detectable effect on the scattered phase function of the particle, primarily in the backward direction. This influence of polarization on the time-shift signals will be addressed again in section 3.

The left picture in Fig. 3 shows the scattering diagram over the solid angle 4*π* when the particle is illuminated by a plane wave with the polarization state [1 -1 0 0]^{T}. This is used to generate the probability density function *ρ* (*θ*, *φ*) for sampling *θ* and *φ*, whose integration over the solid angle 4π should be unity. Therefore, when a large number of the light rays with the polarization state [1 -1 0 0]^{T} are incident on the particle, the scattered light ray density in each direction specified by angle *θ* and *φ* should be proportional to the scattered intensity *I _{sca}* (

*θ, φ*). Since

*θ*and

*φ*are not separable, and the accumulative method cannot be used for sampling, the rejection sampling method is used [7–9]. With Eq. (24), the scattered intensity over the solid angle 4π could be expressed with Eq. (25) [9]

*θ*,

*φ*) could be expressed with Eq. (26) For using the rejection sampling method, the target scattering direction (

*θ*,

*φ*) has a probability density function

*f*(

*θ*,

*φ*). A probability distribution G (

*θ*,

*φ*) needs to be found, corresponding to the probability density function

*g*(

*θ*,

*φ*), for which we already have an efficient generation algorithm. The ratio

*f*(

*θ*,

*φ*)/

*g*(

*θ*,

*φ*) should be bounded by a constant C, and C should be as close to 1 as possible to save computational time. For the embedded nanoparticle with the diameter 320nm and the refractive index 1.6268, the maximum scattered intensity occurs in the forward direction (

*θ*=0) [32]. Therefore, the probability distribution G can be constructed by utilizing the maximum scattered intensity at

*θ=0*, and its

*g*(

*θ*,

*φ*) could be expressed with Eq. (27): For the probability density function

*g*(

*θ*,

*φ*), both of the scattering angles

*θ*and

*φ*are uniformly distributed within [0, π] and [0, 2π]. The angle sampling can be achieved by using Eqs. (28)–(32):

*θ*, a lookup table having a resolution of 0.001° was created for the Mueller matrix elements

*S*, in order to speed up the computation.

_{11}and S_{12}*ξ*,

_{1}*ξ*and

_{2}*ξ*are random numbers. If

_{3}*I*is smaller than

_{com}*I*(

_{sca}*θ*), then the (

_{0}, ϕ_{0}*θ*) values are accepted; otherwise, these steps need to be repeated until the condition is met that

_{0}, ϕ_{0}*I*<

_{com}*I*(

_{sca}*θ*). After the generation of

_{0}, ϕ_{0}*θ*and

*φ*, the Stokes vector needs to be updated with Eq. (33), to modify the probability of the light ray scattering to the next direction.

After the intersection with the colloidal particle, the local reference frame needs to be updated from ${\left[ {\begin{array}{ccc} {e_x^0}&{e_y^0}&{e_z^0} \end{array}} \right]^{T}}$ to ${\left[ {\begin{array}{ccc} {e_x^1}&{e_y^1}&{e_z^1} \end{array}} \right]^T}$, i.e.

## 3. Validation and comparison between the simulation and experiment results

First, the generation of the new scattering direction of the light ray was validated, since this is one of the most essential parts of this Monte Carlo code. To validate this step, a scattering diagram for the particle is computed with the rejection sampling method and compared to results given by the Lorenz-Mie theory. Thirty billion light rays were launched for the computation of the scattering diagram. With the Mueller matrix computed from the Software MiePlot [32] and Eq. (24), the bivariate scattering diagram is plotted in Fig. 4(a), when the incident plane wave is linearly polarized with the Stokes vector [1 -1 0 0]^{T}. Then the scattering diagram is generated with the rejection sampling method for comparison, as shown in Fig. 4(b). The agreement can be seen to be very good.

To underline this good agreement the difference between the two scattering diagrams from the Lorenz-Mie theory and the rejection sampling method is shown in the Fig. 5. This difference is normalized with the local value; hence, given as percent deviation. The maximum normalized deviation occurs where the absolute intensity becomes very low (due to the normalization), but for the most part the deviation remains extremely low. However the deviation can be reduced to any arbitrary level by increasing the number of the incident light rays. For a detector located in the direction of lower scattered intensity, more incident light rays will needed to obtain a converged signal.

With the methodology described in section 2, time-shift signals are simulated for the situation when a colloid suspension droplet passes through a focused laser light sheet [33]. The laser beam from the time-shift instrument is elliptical Gaussian beam; the beam waist in the X-axis is 1mm and the beam waist in the Y-axis is 14µm. The measurement volume of the time-shift instrument is in the Rayleigh range; therefore, the divergence of the Gaussian beam is neglected in the simulation. The detector is located in the backward direction at *θ*=165° and covers a solid angle of 0.0324sr with an opening aperture of 11.6°. The light rays were uniformly detected over the acceptance angle; the scattered intensity is the intensity sum of all rays impinging onto the detector surface. All rays are treated completely incoherent, and the intensity of individual rays is first computed [16].

The measurement principle of the time-shift instrument is illustrated in Fig. 6 [16]; when a pure drop moves along the direction of the black arrow through the highly focused Gaussian beam, each detector located above and below the laser beam detects three signal peaks. By using the time difference between the signal peaks, such as the time difference Δt220 between the signal peak from p = 0 and the signal peak p = 2.2, and knowing the velocity of the droplet, the size of the droplet can be computed. More details about the measurement principle can be found in [34,35].

The simulated time-shift signals have been compared with experimental time-shift signals, in which a colloidal drop with certain size and colloidal concentration passes through the perpendicular polarized laser beam of a time-shift measurement instrument [33]. The colloidal drops were generated using an FMP droplet generator [36], with a size of 100 µm; the diameter of the colloidal particles was 320 ± 30 nm (polystyrene). When the wavelength of the incident light is 405 nm, the extinction efficiency factor *Q _{ext}* is 0.9145. The comparison is made for two volume concentrations: 0.05% and 0.14%, which corresponds to the optical mean free path lengths

*E(L)*466.76 µm and 166.7 µm respectively. The volume concentrations were kept intentionally very low for comparison, because for higher volume concentrations, only the peak of reflection scattering can be clearly recognized; the second-order refractive peaks are no longer identifiable [2], as will also be addressed below.

Figure 7(a) shows the comparison between the simulated and the measured time-shift signal; the amplitude of the signals was normalized with the peak value. The overall agreement between the simulated and measured signal is satisfactory. At this concentration, the three peaks of the time-shift signal could still be recognized; the signal peaks from left to right are p = 0, p = 2.1 and p = 2.2 corresponding to reflection, and the two existing modes of second-order refraction at this refractive index [35].

Figure 7(b) shows the comparison, when the volume concentration of the droplet is 0.14%. Under this concentration, only the signal peak p = 2.1 could be clearly recognized, verifying the stronger internal scattering from the colloidal particles when the volume concentration is higher.

Returning to the remark concerning polarization in section 2.3, Fig. 8 shows the comparison between time-shift signals computed using parallel and perpendicular polarizations. For these simulation conditions the polarization has a minor effect on the simulated time-shift signal.

The results of this section and especially the good agreement achieved between simulated and measured time-shift signals, indicates that the polarized Monte Carlo code reliably captures the light scattering when a colloid suspension droplet passes through a highly focused Gaussian beam. The question now arises, whether this understanding of the light scattering can be exploited in the time-shift signal processing, to estimate not only the size of the droplet, but also the colloid concentration in the droplet. To examine this question further, simulated time-shift signals are examined in more detail and over a larger range of parameter variations, in particular in dependence of drop size, beam waist and colloid concentrations.

## 4. Analysis the time-shift signal from a colloidal suspension droplet

First, the dependence of the time-shift signal on the incident beam waist size is shown in Fig. 9. As expected, the peaks of the individual scattering orders begin to overlap more as the beam waist increases. This confirms the general design rule for the time-shift technique, that the incident beam should be as focused as highly as possible to maintain good peak separation; hence, accurate measurement of the time shifts. The signal level outside the peak regions is rather insensitive to the beam waist, suggesting that the secondary scattering from internal colloidal particles remains about constant: while the illuminated volume in the drop becomes larger with increasing beam waist, the intensity decreases for a given overall light power (intensity = power/area). This lower intensity leads also to a decrease of peak amplitude of the reflection peak. The second-order refractive peak (p = 2.2) decreases to a much larger extent with increasing beam waist. This is understandable, since the second-order refraction scattering undergoes extinction when passing through the droplet, so that both the incident intensity decreases and the extinction must be considered. Note that for the conditions simulated in Fig. 9, the peak for second-order refraction second mode (p = 2.2) is completely masked by internal, secondary scattering, as discussed next.

To separate the effect of internal, secondary scattering from colloid particles from the light scattering from a pure liquid droplet passing through the beam of a time-shift device, a direct comparison is drawn between two such simulated signals, in which the size and refractive index of the host droplet remained constant. The result of these simulations is illustrated in Fig. 10. For this diagram, the code was modified to separate the contribution to the signal coming only from the internal scattering, which is plotted separately with the blue line.

For the colloidal concentration shown in Fig. 10 (0.14% vol., *E*(*L) * =166.7 um) the signal peak of the mode p = 2.1 can still be recognized. This would allow estimation of the droplet size, since the time shift between the p = 2.1 peaks on the two detectors of the time-shift device could be used for this purpose [34]. However, a comparison of the time-shift signal from a pure liquid droplet with the scattered signal from the second-order refraction scattering indicates that the signal is significantly attenuated at this droplet position in the laser beam, again, arising from extinction through internal, secondary scattering. This attenuation will now be quantitatively estimated.

The intensity of the p = 2.1 peak above the internal scattering background is designated as *I* in Fig. 10. The length of the optical mean path (*OL*) within the droplet for second-order refraction scattering is related to the droplet size, the relative refractive index between the medium and the droplet, and the incident angle for second-order refraction scattering (Fig. 1). This length *OL* is given by

*m*is the relative refractive index between the droplet and surrounding medium,

*θ*is the incident angle for the mode p = 2.1 and

_{i}*R*is the radius of the droplet. The attenuation of the signal is expected to follow the Beer-Lambert law [7,37,38], which reads as in which

*I*is the intensity of the signal peak from second-order refraction scattering from a pure droplet with the same size and the refractive index. Under the condition for the simulation in the Fig. 10, the ratio

_{0}*I/I*should be 0.3088 according to the Beer-Lambert law. This ratio from the simulation results, 0.3149, is consistent with the theoretical value; the deviation is less than 2%. This agreement holds true over a range of simulated colloid concentrations 0.05%<vol. conc. < 0.25% (80 µm <

_{0}*E*(

*L*)

*<*466.7 μ

*m*). For higher volume concentrations, the signal peak from second-order refraction scattering attenuates strongly, and could no longer be separated from the simulated time-shift signal, as will be shown in Fig. 12. For lower colloidal concentrations, by knowing the attenuation of the second-order refraction scattering, the optical mean free path length as well as the colloidal concentration could be estimated. The attenuation of the signal peak is used in the work of Onofri et al. [3] to estimate the colloidal concentration of the colloidal particle as well. Note that in the simulations shown in both Fig. 9 and Fig. 10 the p = 2.2 peak is masked by internal, secondary scattering, whereas the p = 2.1 remains identifiable. This effect can be understood by considering the incident areas on the drop, over which these scattering order peaks are generated. This is illustrated in Fig. 10, showing regions of the particle, in which the light rays of the laser beam propagate in the Z direction. The detector is located in the lower backward direction at

*θ*=165° and covers a solid angle of 0.0324 sr with an opening aperture of 11.6°. In the left picture of the Fig. 11, the angular incident regions and exit regions for the modes p = 0, p = 2.1, p = 2.2 are illustrated with thicken lines on the circumference of the droplet, with the colors red, blue and black respectively. When the light rays in the YZ plane are incident on these angular glare regions, the exit light rays will fall on the detector.

As the left picture shows, the angular incident region for the mode p = 2.1 is much larger than for p = 2.2. The right picture in the Fig. 11 shows the incident areas for the laser sheet in the XY plane for the three signal peaks. Also in this plane, the incident area for the mode p = 2.1 is the largest. When the light rays impinge on a colloidal particle, as Fig. 4 depicts, the forward (*θ*=0) scattering exhibits the highest amplitude. Therefore, in the light ray path for second-order refraction scattering, when a light ray interacts with a colloidal particle, the scattered ray is most likely to follow the original light ray path and be registered by the detector. This explains why the highest amplitude of internal scattering occurs at the same position as the signal peak p = 2.1.

Finally, the effect of increasing the colloidal concentration will be considered. Figure 12 shows simulated time-shift signals for a higher volume concentration, 2% (*E*(*L*) = 10 µm). In this signal, only the signal peak from p = 0 can be observed. This signal resembles closely the signal structure observed by Rosenkranz et al. [2] using a two-dimensional simulation, in which the *E*(*L*) is also set to 10µm. The signal structure is significantly different than the signal observed for lower volume concentrations (Fig. 7), since neither of the second-order refractive peaks are identifiable. This is a result of two counteracting effects: the peak amplitudes are decreased due to extinction and simultaneously the internal, secondary scattering at these particular drop positions in the beam have increased. However, the reflective peak is not susceptible to attenuation through extinction, so that the peak amplitude consists of the undisturbed reflection plus the background level from internal, secondary scattering, as marked in Fig. 12.

A ratio *η* is defined to characterize the relative intensity of the internal scattering

*η*is therefore dependent on the volume concentration of the colloidal particle; for a droplet with same size, as the volume concentration of the colloidal particle increases, the internal scattering from the colloidal particle becomes stronger and

*η*will increase. This then represents a viable method of determining colloid suspension concentrations at higher loading values. The main remaining difficulty is that the ratio will also depend on the drop size, since the internal secondary scattering intensity will change according to the path length through the drop. However, with the simulation tool presented in this study, this dependency can be computed beforehand. Indeed, it is not even necessary to simulate the entire time-shift signal, but only the signal amplitude just before and during the drop position for reflective scattering. This remains a topic of future investigation.

## 5. Summary and conclusions

In this study, a Monte Carlo approach is used to simulate the polarized laser light scattering from a colloidal suspension droplet. The Stokes vector is tracked throughout the entire scattering process. For internal scattering at the drop surface, total internal reflection can occur and for this situation, the Mueller matrix formulation is derived and implemented. After the scattering event with each colloidal particle, a new propagation direction is generated using the rejection sampling method. The code was applied to simulate time-shift signals from colloidal droplets, resulting in very good agreement with experimentally obtained results.

Further efforts will be made to estimate the volume concentration of the colloidal particle. For the lower colloidal concentration, it is possible to use the attenuation ratio of second-order refraction scattering to estimate the colloidal concentration. For higher colloidal concentration, the ratio *η* may be used to characterize the colloidal concentration; however, *η* is dependent also on the droplet size. Therefore, the relationship between *η* and the droplet size still needs to be investigated.

Several extensions of the present code are planned in order to address practical measurement situations. For one, time-shift signals obtained for polydispersed colloid particles will be examined, starting with bimodal size distributions. The aim would be to identify signal features which may be useful in estimating also the size of the colloidal particles. A second extension foresees a non-uniform spatial distribution of the colloidal particles in the droplet. This situation occurs in droplets with metallic nano-particles, in which the settling velocity of the particles is significant. In this case the aim would be to identify signal features alluding to a non-homogeneous spatial distribution of the colloidal particles in the droplet.

## Funding

Deutsche Forschungsgemeinschaft (SFB-TRR 75); Bundesministerium für Bildung und Forschung (ZIM ZF4415101LT7).

## Acknowledgments

The authors would like to thank the Federal Ministry for Economic Affairs and Energy of Germany for financial support through the Central Innovation Programme for SMEs (ZIM), contract ZF4415101LT7. They are grateful for the very fruitful discussion with the Prof. Feng Xu from the University of Oklahoma and the help for the experiment from Can Li from Zhejiang University. They also appreciate the cooperation of the High-Performance Computing Center of the TU Darmstadt with the project ID 960.

## Disclosures

The authors declare no conflicts of interest

## References

**1. **C. Tropea, “Optical particle characterization in flows,” Annu. Rev. Fluid Mech. **43**(1), 399–426 (2011). [CrossRef]

**2. **S. Rosenkranz, W. Schäfer, C. Tropea, and A. M. Zoubir, “Modeling photon transport in turbid media for measuring colloidal concentration in drops using the time-shift technique,” Appl. Opt. **55**(34), 9703–9711 (2016). [CrossRef]

**3. **F. Onofri, L. Bergougnoux, J. L. Firpo, and J. Misguich-Ripault, “Size, velocity, and concentration in suspension measurements of spherical droplets and cylindrical jets,” Appl. Opt. **38**(21), 4681–4690 (1999). [CrossRef]

**4. **T. Wriedt and R. Schuh, “The inclusion-concentration measurement of suspension droplets based on Monte Carlo ray tracing,” Meas. Sci. Technol. **13**(3), 276–279 (2002). [CrossRef]

**5. **D. Jakubczyk, G. Derkachov, M. Kolwas, and K. Kolwas, “Combining weighting and scatterometry: Application to a levitated droplet of suspension,” J. Quant. Spectrosc. Radiat. Transfer **126**, 99–104 (2013). [CrossRef]

**6. **P. G. Stegmann, B. Sun, J. Ding, P. Yang, and X. Zhang, “Study of the effects of phytoplankton morphology and vertical profile on lidar attenuated backscatter and depolarization ratio,” J. Quant. Spectrosc. Radiat. Transfer **225**, 1–15 (2019). [CrossRef]

**7. **J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Express **13**(12), 4420–4438 (2005). [CrossRef]

**8. **F. Jaillon and H. Saint-Jalmes, “Description and time reduction of a Monte Carlo code to simulate propagation of polarized light through scattering media,” Appl. Opt. **42**(16), 3290–3296 (2003). [CrossRef]

**9. **B. Kaplan, G. Ledanois, and B. Drévillon, “Mueller matrix of dense polystyrene latex sphere suspensions: measurements and Monte Carlo simulation,” Appl. Opt. **40**(16), 2769–2777 (2001). [CrossRef]

**10. **M. Xu, “Electric field Monte Carlo simulation of polarized light propagation in turbid media,” Opt. Express **12**(26), 6530–6539 (2004). [CrossRef]

**11. **M. J. Raković, G. W. Kattawar, M. Mehrűbeoğlu, B. D. Cameron, L. V. Wang, S. Rastegar, and G. L. Coté, “Light backscattering polarization patterns from turbid media: theory and experiment,” Appl. Opt. **38**(15), 3399–3408 (1999). [CrossRef]

**12. **Y. X. Hu, D. Winker, P. Yang, B. Baum, L. Poole, and L. Vann, “Identification of cloud phase from PICASSO-CENA lidar depolarization: A multiple scattering sensitivity study,” J. Quant. Spectrosc. Radiat. Transfer **70**(4-6), 569–579 (2001). [CrossRef]

**13. **I. Lux, * Monte Carlo Particle Transport Methods*. (CRC, 1991).

**14. **E. Collett, “Mueller-Stokes matrix formulation of Fresnel's equations,” Am. J. Phys. **39**(5), 517–528 (1971). [CrossRef]

**15. **L. Li, S. Rosenkranz, W. Schäfer, and C. Tropea, “Sensitivity of the time-shift technique in characterizing non-spherical drops,” In: Proceedings of the nineteenth International Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics, Lisbon; 16-19 July. 2018.

**16. **L. Li, S. Rosenkranz, W. Schäfer, and C. Tropea, “Light scattering from a drop with an embedded particle and its exploitation in the time-shift technique,” J. Quant. Spectrosc. Radiat. Transfer **227**, 20–31 (2019). [CrossRef]

**17. **P. G. Stegmann, C. Tropea, E. Järvinen, and M. Schnaiter, “Comparison of measured and computed phase functions of individual tropospheric ice crystals,” J. Quant. Spectrosc. Radiat. Transfer **178**, 379–389 (2016). [CrossRef]

**18. **M. Wendisch and P. Yang, * Theory of Atmospheric Radiative Transfer: A Comprehensive Introduction*. (John Wiley & Sons, 2012).

**19. **G. Mie, “Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen,” Ann. Phys. **330**(3), 377–445 (1908). [CrossRef]

**20. **J. P. Briton, B. Maheu, G. Gréhan, and G. Gouesbet, “Monte Carlo simulation of multiple scattering in arbitrary 3-D geometry,” Part. Part. Syst. Charact. **9**(1-4), 52–58 (1992). [CrossRef]

**21. **M. P. Sentis, F. Onofri, O. Dhez, J. Y. Laurent, and F. Chauchard, “Organic photo sensors for multi-angle light scattering characterization of particle systems,” Opt. Express **23**(21), 27536–27541 (2015). [CrossRef]

**22. **K. F. Ren, F. Onofri, C. Rozé, and T. Girasole, “Vectorial complex ray model and application to two-dimensional scattering of plane wave by a spheroidal particle,” Opt. Lett. **36**(3), 370–372 (2011). [CrossRef]

**23. **H. Yu, F. Xu, and C. Tropea, “Simulation of optical caustics associated with the secondary rainbow of oblate droplets,” Opt. Lett. **38**(21), 4469–4472 (2013). [CrossRef]

**24. **H. C. Hulst, * Light Scattering by Small Particles*. (Wiley, 1981).

**25. **R. M. Azzam, “Phase shifts that accompany total internal reflection at a dielectric–dielectric interface,” J. Opt. Soc. Am. A **21**(8), 1559–1563 (2004). [CrossRef]

**26. **C. F. Bohren and D. R. Huffman, * Absorption and Scattering of Light by Small Particles*. (John Wiley & Sons, 2008).

**27. **E. Hecht, * Optics*. (Addison-Wesley, 2002).

**28. **H. Yu, J. Shen, and Y. Wei, “Geometrical optics approximation of light scattering by large air bubbles,” Particuology **6**(5), 340–346 (2008). [CrossRef]

**29. **M. P. Sentis, F. Onofri, L. Méès, and S. Radev, “Scattering of light by large bubbles: Coupling of geometrical and physical optics approximations,” J. Quant. Spectrosc. Radiat. Transfer **170**, 8–18 (2016). [CrossRef]

**30. **F. A. Jenkins and H. E. White, * Fundamentals of Optics*. (Tata McGraw-Hill Education, 1976).

**31. **M. I. Mishchenko and P. Yang, “Far-field Lorenz–Mie scattering in an absorbing host medium: theoretical formalism and FORTRAN program,” J. Quant. Spectrosc. Radiat. Transfer **205**, 241–252 (2018). [CrossRef]

**32. **P. Laven, “Simulation of rainbows, coronas, and glories by use of Mie theory,” Appl. Opt. **42**(3), 436–444 (2003). [CrossRef]

**33. **AOM-Systems, https://www.aom-systems.com/de/

**34. **W. Schäfer and C. Tropea, “Time-shift technique for simultaneous measurement of size, velocity and relative refractive index of transparent droplets or particles in a flow,” Appl. Opt. **53**(4), 588–596 (2014). [CrossRef]

**35. **H. E. Albrecht, N. Damaschke, M. Borys, and C. Tropea, * Laser Doppler and Phase Doppler Measurement Techniques*. (Springer Science & Business Media, 2013).

**36. **FMP-Technology GmbH, https://fmp-technology.com/?lang=en

**37. **C. Li, X. C. Wu, J. Z. Cao, L. H. Chen, G. Gréhan, and K. F. Cen, “Application of rainbow refractometry for measurement of droplets with solid inclusions,” Opt. Laser Technol. **98**, 354–362 (2018). [CrossRef]

**38. **M. F. Modest, * Radiative Heat Transfer*. (Academic, 2013).