## Abstract

The effect of dependent scattering on the bulk scattering properties of intralipid phantoms in the 600-1850 nm wavelength range has been investigated. A set of 57 liquid optical phantoms, covering a wide range of intralipid concentrations (1-100% *v/v*), was prepared and the bulk optical properties were accurately determined. The bulk scattering coefficient as a function of the particle density could be well described with Twersky’s packing factor (R^{2} > 0.990). A general model was elaborated taking into account the wavelength dependency and the effect of the concentration of scattering particles (R^{2} = 0.999). Additionally, an empirical approach was followed to characterize the effect of dense packing of scattering particles on the anisotropy factor (R^{2} = 0.992) and the reduced scattering coefficient (R^{2} = 0.999) of the phantoms. The derived equations can be consulted in future research for the calculation of the bulk scattering properties of intralipid dilutions in the 600-1850 nm range, or for the validation of theories that describe the effects of dependent scattering on the scattering properties of intralipid-like systems.

© 2014 Optical Society of America

## 1. Introduction

Propagation of Visible (Vis) and Near Infrared (NIR) light through turbid media, like biological tissues, can be described by the radiative transfer theory (RTT), taking into account the medium’s bulk optical properties (BOP): bulk absorption coefficient *µ _{a}*, bulk scattering coefficient

*µ*, reduced scattering coefficient

_{s}*µ*, and scattering phase function

_{s}’*p(θ)*or derived anisotropy factor

*g*. In biomedical optics and spectroscopic applications, knowledge about the product’s BOP is, therefore, essential for a correct understanding of measured optical signals [1,2]. Accordingly, accurate calibration, validation and comparison of measurement setups with optical phantoms, for which the BOP’s are known accurately, is crucial. Optical phantoms with variable BOP’s can be constructed by mixing an absorbing dye, a scattering agent and a (neutral) matrix in different ratios.

*Intralipid*(IL) is a fat-in-water emulsion which is frequently used in biomedical optics as scattering component for the preparation of liquid optical phantoms [3]. Since the majority of the fat globules in IL have a diameter smaller than 500 nm, scattering is high in the Vis and decreases with increasing wavelength [4–12]. Recent studies demonstrated that the batch-to-batch variability in scattering properties of IL is very small and that IL is highly stable over a long time and at different temperatures [4,5,13]. This confirms that IL is a very suitable scattering component for liquid optical phantoms.

^{®}Variability in scattering levels (*µ _{s}* and

*µ*) of the optical phantoms can be created by diluting IL with water (scattering neutral matrix) [3]. If the scattering events are independent and dilution with the matrix (and absorbing dye) does not influence the optical properties of the single scattering particles,

_{s}’*µ*and

_{s}*µ*follow a linear decrease (through the origin) with reducing IL concentration [6]. Moreover, under these conditions, the normalized scattering phase function and the anisotropy factor are independent of the dilution factor. If these conditions are met, the production process for optical phantoms is straightforward as the mixing ratios required to obtain the desired BOP’s can be easily calculated from the known optical properties of the individual components [8,9,12]. Di Ninni

_{s}’*et al.*[6] already reported that the microphysical properties of the single scattering particles in IL are not affected by dilution with water. Furthermore, the results of their experimental investigation showed that dilution only affects absorption according to the linear Beer-Lambert relation. However, if the concentration of scattering particles is high, individual scattering events can influence each other, especially at longer wavelengths [6,11,12,14]. This phenomenon is often referred to as dependent or correlated scattering, which can be explained by coherent addition of the radiation amplitudes of single scattering events (far-field interaction). For small volume fractions, the arrangement of the scattering particles can be assumed random and at sufficient distance, but this is no longer true at higher volume fractions [15].

In the analytical or multiple scattering theory, which starts from Maxwell equations or wave equations, various multiple scattering processes in a medium can be described with inclusion of dependent scattering, diffraction and interference effects. Although these theories are considered the most rigorous, they are too complex for practical use. Therefore, theories like Twersky’s theory, Foldy-Lax formulation and *T*-matrix computations have been derived [16–21]. However, the computations for these theories are still very intensive and knowledge of the dielectric properties at the microscale is required [2,22]. This information is often not available and difficult to determine.

RTT, on the other hand, deals directly with the transport of energy [23]. It is a conceptually simpler approach, approximating the tissue as a homogeneous layer defined by its BOP’s. Nevertheless, it has been found to accurately describe light propagation in turbid media [1]. Since no correlation between fields is assumed in RTT, dependent scattering is not directly accounted for [24]. Therefore, the effect of dependent scattering is included in the bulk scattering properties (*µ _{s}*,

*µ*,

_{s}’*g*and

*p(θ)*). For this reason,

*µ*and

_{s}*µ*deviate from the linear relation with the volume concentration of scattering particles

_{s}’*Φ*for a specific scattering agent. This deviation typically increases with increasing

_{p}*Φ*. Additionally, the normalized

_{p}*p(θ)*and

*g*become dependent on

*Φ*[11].

_{p}Several approximations have been derived to demonstrate the effect of dependent scattering on the bulk scattering coefficient [14,15,20,25–27]. However, a closed form solution was only derived by Twersky for small particles compared to the considered wavelength [28]. A correlation-packing factor *W _{p}* was determined by reducing the integral of the statistical pair correlation function for impenetrable particles. The pair correlation function characterizes the radial distribution of particles in a three-dimensional space. Hence, the derived packing factor is a simple rational function of the fraction of space occupied by the correlated statistical mechanics particles (volume concentration

*Φ*):

_{p}This packing factor *W _{p}* decreases monotonously from unity to zero as

*Φ*increases from zero to unity, representing the deviation from the linear relation between

_{p}*µ*and

_{s}*Φ*[29]. Apart from

_{p}*Φ*, Twersky’s packing factor only depends on the packing dimension

_{p}*p*[29–31]. The packing dimension describes the rate at which the empty space between the scatterers reduces as the total number density increases [31].

The packing dimension for monodisperse scattering slabs, cylinders and spheres, much smaller than the wavelength, was found to be 1, 2 and 3 respectively [29,31]. Additionally, Twersky derived a correction factor for the case of polydisperse scatterers, resulting in a packing factor closer to one and a reduction of the falloff with increasing *Φ _{p}* [32,33]. However, only size distributions governed by the Poisson distribution were considered. Bascom and Cobbold proposed a new theoretical model with the same formulation as Eq. (1), but with a packing dimension representing both the particle shape and polydispersity [30]. An increased degree of polydispersity results in a reduction in packing symmetry and, therefore, a reduction of the packing dimension. The key feature of this approach is its simplicity as it describes dependent scattering with a single-parameter equation. However, the packing factor does not uniquely define the pair correlation function as it is related to the integral of this function and not to the function itself. This means that two different pair correlation functions could result in the same packing factor [30]. Bascom and Cobbold successfully fitted the packing factor with variable packing dimension to the acoustic backscattering coefficient of red blood cells as a function of the cell count. Moreover, a decreased particle symmetry was related to a reduction in the packing dimension [30]. In several other studies, the Twersky equation with a packing dimension of 3 was successfully consulted to account for dependent scattering in Mie-calculations of the Vis and NIR bulk scattering coefficient

*µ*for IL dilutions and biological tissues [31,34].

_{s}The effect of dependent scattering on the scattering phase function and the anisotropy factor has, however, not been studied intensively. Mishchenko used the Percus-Yevick approximation [35], while Mackowski and Mishchenko consulted their *T*-matrix computations to investigate the effect of the compaction state of monodisperse scattering spheres on the anisotropy factor *g* [18]. Both studies showed a reduction in *g* with an increase of *Φ _{p}*, especially if the particle radii were in the range of the wavelength.

Dependent scattering in IL dilutions has been studied at a few discrete wavelengths in the Visible and short-wave NIR [6,11,14]. Di Ninni *et al.* measured the reduced scattering coefficient (*µ _{s}’*) at 751 and 833 nm for IL dilutions with a

*Φ*up to 0.032 ml/ml. The relation between

_{p}*Φ*and

_{p}*µ*could be well explained by second-order polynomials (through the origin) and showed a deviation within 2% from the linear relation (through the origin) for

_{s}’*Φ*< 0.023 ml/ml [6]. Zaccanti

_{p}*et al.*and Giusto

*et al.*studied the bulk scattering properties (

*µ*,

_{s}*µ*and

_{s}’*g*) at 633 nm for IL dilutions with a

*Φ*up to 0.227 ml/ml (pure IL) [11,14]. Both

_{p}*µ*and

_{s}*µ*followed a second-order polynomial (through the origin) as a function of

_{s}’*Φ*, while

_{p}*g*followed a nearly linear decrease. Unfortunately, no information is available about the effect of dependent scattering on the scattering properties of IL dilutions for wavelengths above 850 nm. However, for spectroscopic applications, the longer wavelength range is more valuable because of the presence of many important molecular overtones and combination vibrations in the NIR [36,37]. Since scattering by IL in the NIR is relatively low, high concentrations are needed in the optical phantoms to mimic the relatively high

*µ*and/or

_{s}*µ*of biological tissues or other turbid samples. Additionally, the effect of dependent scattering is even more dominant for longer wavelengths [6,11,12,14,27,31]. For these reasons, dependent scattering is expected to be more significant in IL phantoms in the NIR range and has to be taken into account in the calculation of the BOP for a specific mixing ratio of the different phantom components [38].

_{s}’The main objective of this research was to study and describe the effect of dependent scattering on the bulk scattering properties of IL phantoms in the 600-1850 nm wavelength range. Semi-empirical equations were derived, describing the bulk scattering properties (*µ _{s}*,

*µ*and

_{s}’*g*) as a function of the wavelength

*λ*and the volume concentration of scattering IL particles

*Φ*.

_{p}## 2. Materials and methods

#### 2.1 Liquid optical phantoms

The set of 57 liquid optical phantoms consulted in this study is the same set as used in [12], covering a wide range of IL concentrations (1-100% *v/v*). Methylene Blue (M9140, Sigma-Aldrich, Missouri, USA), *Intralipid ^{®} 20%* (batch 10FH1726, expiring date 07/2014, Fresenius Kabi, Germany) and deionized water were respectively used as absorbing, scattering and dilution agent and mixed in different ratios. Methylene Blue (MB) was chosen as absorber since it has a sharp absorption peak in the 550-750 nm range where absorption by water and IL is minimal (both

*µ*< 0.03 cm

_{a}^{−1}) [8]. Taking into account the density of the different IL components, pure IL contains 22.7% (

*v/v*) scattering particles (soybean oil + egg lecithin) [8–11]. 56 phantoms were created by combining 8 absorption levels, corresponding to MB concentrations of 0, 4, 8, 16, 32, 64, 128 and 136 µM, and 7 scattering levels, corresponding to IL concentrations of 1, 2, 4, 8, 16, 32 and 64% (

*v*/

*v*) or a

*Φ*of respectively 0.00227, 0.00454, 0.00908, 0.01816, 0.03632, 0.07264 and 0.14528 ml/ml. Additionally, pure IL (

_{p}*Φ*= 0.227 ml/ml) was also measured.

_{p}#### 2.2 Optical characterization of phantoms

The diffuse reflectance and total transmittance for all phantoms was measured in a DIS setup, while the unscattered transmittance was measured in a separate measurement path. The measurement setup was especially designed to obtain high signal-to-noise spectra in the 500 – 2250 nm wavelength range for very turbid media. An extensive description about the measurement setup and calibration and measurement procedure can be found in [12].

All liquid phantoms were freshly prepared and measured at room temperature (22 ± 0.5°C) for two sample thicknesses (*d*): 0.55 and 1.1 mm. Measurement of the same sample at two thicknesses should finally result in the same estimated BOP values. Therefore, it provides an interesting means to evaluate the measurement setup and the BOP estimation routine [39]. First, the sample was measured from 500 to 2250 nm with an interval step of 5 nm in the DIS measurement path, followed by measurement in the unscattered transmittance path. All different phantoms and both thicknesses were measured randomly in time. The cuvette was cleaned and dried thoroughly before loading the next sample.

The inverse adding doubling (IAD) program developed and optimized by Prahl [39] was consulted for the estimation of the BOP values from the diffuse reflectance and total and unscattered transmittance measurements. Apart from the three measurements, also the real refractive index (*n*) of the sample has to be provided to the algorithm and was therefore calculated based on the volume concentration of scattering particles (*Φ _{p}*) in the sample:

*n*=

_{sample}*n*+ 0.14

_{water}*Φ*[8,9]. The wavelength-dependent real refractive index of water at room temperature was adopted from Hale and Querry [40]. The effect of MB on the refractive index of the sample was neglected because of the very low concentrations.

_{p}The bulk absorption (*µ _{a}*) and reduced scattering coefficient (

*µ*) could be estimated for the entire range. However, separation between

_{s}’*µ*and

_{s}*g*could not be established when the scattering depth (

*µ**

_{s}*d)*was very high. The reason for this was the presence of a significant amount of scattered photons in the measurement of the unscattered transmittance [12]. In general, good estimations for

*µ*and

_{s}*g*were obtained for

*λ*≥ 560 + 2200

*Φ*in the case of 0.55 mm sample thickness and

_{p}*λ*≥ 600 + 3000

*Φ*for 1.1 mm sample thickness (

_{p}*λ*in nm and

*Φ*in ml/ml). Moreover, because of low signal levels in the diffuse reflectance measurements for wavelengths above 1200 nm (high absorption by water) for phantoms with very low scattering levels (

_{p}*Φ*≤ 0.02 ml/ml), no accurate estimates for the anisotropy factor could be established for that specific range [12]. As for the remaining range, the estimated BOP’s for the phantoms measured at two different sample thicknesses were nearly the same, it was concluded that the measurement and BOP estimation procedure were accurate [12]. Hence, the BOP’s of the same phantom for both thicknesses were treated here as replicates. Moreover, as the separation between scattering and absorption was very successful, with nearly no cross-talk, the phantoms with different absorption levels for a fixed scattering level were also treated as replicates for the bulk scattering properties (

_{p}*µ*,

_{s}*g*and

*µ*).

_{s}’#### 2.3 Equation fitting to the bulk scattering properties

For fitting the equations, only the wavelength range 600-1850 nm was considered. Below 600 nm, no estimates for *µ _{s}* and

*g*were available for phantoms with

*Φ*≥ 0.018 ml/ml since the unscattered transmittance measurements were incorrect. Above 1850 nm, the estimated BOP’s for low scattering levels (

_{p}*Φ*≤ 0.03 ml/ml) showed large variability. This can be explained by the low signal levels in the diffuse reflectance measurements, due to combination of the high absorption by water and the low scattering for low IL concentrations [12].

_{p}The Twersky equation [Eq. (1)] describes the deviation from the linear ‘independent scattering’ relation between *µ _{s}* and

*Φ*. In order to describe the direct relation between

_{p}*µ*and

_{s}*Φ*, including dependent scattering, the Twersky equation needs to be multiplied with this linear relation. An accurate equation for the

_{p}*µ*spectrum of pure IL in the 500-2250 nm range, in absence of dependent scattering, was derived in [12], with

_{s}*µ*in cm

_{s}^{−1}and

*λ*in nm:

The linear relation between *µ _{s}* and

*Φ*for independent scattering can be obtained by correcting Eq. (2) with 0.227 ml/ml (

_{p}*Φ*for pure IL) and multiplying with the actual

_{p}*Φ*. The final equation, describing

_{p}*µ*as a function of

_{s}*λ*and

*Φ*, taking into account dependent scattering by including the packing factor, can be written as follows:

_{p}As *p* is unknown for IL in the Vis-NIR range, Eq. (3) was evaluated by fitting a *p*-parameter at each wavelength. The fitted *p*-parameters were then used in the further analysis.

Since no closed form solution describing the effect of dependent scattering on the anisotropy factor for densely packed scattering media was found in literature, an empirical approach was followed. For very low IL concentrations (*Φ _{p}* → 0 ml/ml), independent scattering is valid and the

*g*-value only depends on the wavelength (

*λ*). In [12], the following equation was derived for this relation in the 500-2250 nm range, with

*λ*in nm:

A more general equation for the anisotropy factor *g*(*λ*, *Φ _{p}*) in the 600-1850 nm range, also valid for the case of dependent scattering (

*Φ*> 0.02 ml/ml), will be determined empirically after visual inspection of the data, starting from the formula for

_{p}*g*(

^{indep}*λ*) [Eq. (4)].

The derived equations were fitted to the measured BOP spectra in Matlab (version 7.10, The Mathworks Inc., Massachusetts, USA) with the curve fitting tool (cftool) for the case of 1 independent variable or the surface fitting tool (sftool) for 2 independent variables.

To describe the effect of dependent scattering on *µ _{s}’*, the derived general equations for

*µ*[Eq. (3)] and

_{s}*g*were substituted in the similarity relation [41]:

## 3. Results and discussion

This manuscript presents a more thorough analysis of the effect of dependent scattering on the estimated BOPs of IL phantoms presented in [12]. For the measured diffuse reflectance and total and unscattered transmittance spectra, as well as an overview and detailed description on the estimated BOP’s, the reader is referred to that publication. Moreover, in the results presented in [12] no effect of dependent scattering on the bulk absorption coefficient (*µ _{a}*) can be observed, which confirms the conclusion presented in [6].

#### 3.1 Bulk scattering coefficient

In Figs. 1(a) and 1(b), the estimated bulk scattering coefficient (*µ _{s}*) values are presented as a function of the volume concentration of scattering particles

*Φ*for respectively 3 wavelengths below and 3 wavelengths above 1000 nm.

_{p}The average of 16 replicates (2 sample thicknesses x 8 absorption levels) for each scattering level is plotted as a cyan dot surrounded by error bars, indicating the standard deviation. The error bars are, however, rather small, indicating the high repeatability of the measurements and BOP estimation procedure. For phantoms with high IL concentrations, the scattering (depth) was very high, especially for short wavelengths (Vis). As a result, scattered photons were captured during the unscattered transmittance measurements and no correct separation between *µ _{s}* and

*g*could be established. Therefore,

*µ*data for short wavelengths were only available for low IL concentrations [Fig. 1(a)]. Equation (3), with a variable packing dimension

_{s}*p*, was fitted successfully (all R

^{2}> 0.996) to the

*µ*data for all 6 wavelengths [Figs. 1(a) and 1(b), solid lines]. In Fig. 1(a), the second-order polynomial derived for IL at 633 nm by Zaccanti

_{s}*et al.*[11] has also been plotted. A very good agreement between this model, the data and the fitted Twersky model can be noticed.

The fitting procedure was repeated for all wavelengths individually, resulting in the *p*- and R^{2}-values demonstrated in respectively Figs. 2(a) and 2(b). Overall, Eq. (3) fitted very well at every wavelength (all R^{2} > 0.990). The dip in R^{2} around 1450 nm is probably caused by a minor water absorption bump in the bulk scattering spectra for almost all phantoms [12], indicating some minor cross-talk between *µ _{a}* and

*µ*(or

_{s}*µ*). The baseline in the

_{s}’*µ*values for this wavelength range introduces an underestimation of the actual

_{s}*p*-values as the Twersky equation is forced to go through the origin, resulting in lower R

^{2}values. Taking this into consideration, the packing dimension follows a nearly linear increase with increasing wavelength above 850 nm [Fig. 2(a)]. However, for the wavelengths below 850 nm, the fitted

*p*-values are highly unstable and the relation with

*λ*is far from linear. This is because no

*µ*data is available in this wavelength range for the phantoms with high IL concentrations [Fig. 1(a)]. Hence, the relation between

_{s}*µ*and

_{s}*Φ*is fairly linear for the remaining points (

_{p}*Φ*< 0.08 ml/ml) and a neighboring

_{p}*p*-value will fit almost equally good. Accordingly, it was assumed that the linear relation between

*p*and

*λ*could be extended for wavelengths below 850 nm. The fitted linear function (R

^{2}= 0.827) is given by the equation and solid line in Fig. 2(a). Substitution of this linear function and Eq. (2) in Eq. (3) provides a general equation for

*µ*, taking into account the effect of both

_{s}*λ*and

*Φ*. This general equation can be represented by a surface in a 3D space, as shown in Fig. 3. Next to the surface, also the measured

_{p}*µ*values are presented (grey dots) and a very good agreement was found between the data and the model (R

_{s}^{2}= 0.999).

The reason why the packing dimension was not constant for all wavelengths is most likely because the scattering particles are only slightly smaller (< 500 nm) compared to the considered wavelengths (600-1850 nm). The main assumption in the derivations by Twersky is, however, that they should be much smaller [29]. Nevertheless, very good fits could be established at all wavelengths [Fig. 2(b)]. Moreover, it seems that the ratio of particle size to wavelength can be taken into account in the packing dimension, next to the effect of particle shape and polydispersity [30]. For a decreasing wavelength, becoming closer to the size of the scattering particles, the packing dimension decreases [Fig. 2(a)]. This indicates a reduced effect of dependent scattering, characterized by a packing factor closer to one and a reduced falloff of the linear relation between *µ _{s}* and

*Φ*. This decreasing effect with decreasing wavelength was also noticed in other studies [6,11,12,14,27,31]. Therefore, it is expected that for decreasing wavelengths below 600 nm, the packing dimension will converge to one. For monodisperse spheres much smaller than the wavelength, the packing dimension should theoretically be equal to 3 [31]. However, as the fat globules in IL are polydisperse, a smaller packing factor can be expected [30,32,33]. As the wavelength increases towards 1850 nm, it becomes much larger than the scattering particles in IL such that the main assumption for Twersky’s equation is close to being realized. This can be noticed from the increase of the packing dimension (towards 3) for increasing wavelength [Fig. 2(a)]. However, since Twersky’s theory proposes a packing dimension smaller than 3 for polydisperse spherical particles much smaller than the wavelength, it is expected that the fitted linear relation between the packing dimension and wavelength will not continue for wavelengths above 1850 nm and that the packing dimension will converge to a stable value between 2.3 (

_{p}*p*at 1800-1850 nm) and 3.

#### 3.2 Anisotropy factor

In order to describe the effect of dependent scattering on the anisotropy factor *g* of IL, a conceptually similar approach was followed. Inspection of the measured values of *g* for dilutions of IL in [11] and [12] shows a nearly linear decrease of *g* with increasing *Φ _{p}* in the 600-1850 nm range. This linear relation can be clearly seen in Figs. 4(a) and 4(b) for 6 wavelengths. The intercept of this linear relation (

*Φ*→ 0 ml/ml) is related to the

_{p}*g*-value for the case of independent scattering [Eq. (4)]. Given a wavelength dependent slope

*k*, the empirical general equation for

*g*as a function of

*Φ*and

_{p}*λ*, taking into account dependent scattering can be written as:

For each individual wavelength, the linear function [Eq. (6)] with predetermined intercept and variable slope *k* was fitted to *g* as a function of *Φ _{p}*. The results for 3 wavelengths below and 3 wavelengths above 1000 nm are presented in respectively Figs. 4(a) and 4(b). It can be concluded that the fitted linear equations describe the effect of dependent scattering on

*g*well. Comparison of our results at 635 nm with the results reported by Zaccanti

*et al.*[11] for

*g*as a function of

*Φ*at 633 nm shows a good agreement. Nevertheless, one might argue that for

_{p}*Φ*going to zero, a convergence to a stable value would be more appropriate. Moreover, following these linear equations for

_{p}*Φ*above 0.227 ml/ml (pure IL),

_{p}*g*would become negative which is very rare for biological tissues or fluids [35]. So, the underlying relation is probably more complex, but not enough data were available to fit a more complex function.

In Fig. 5(a) the fitted slope parameters at each wavelength are presented, while the R^{2} for each fit is illustrated in Fig. 5(b). The R^{2} values corresponding to the linear equations for wavelengths outside the 800-1400 nm range indicate a rather poor fit. This can be explained by the relatively large variation in *g* for a specific *Φ _{p}* compared to the variability of the average

*g*between different levels of

*Φ*, characterized by a slope close to 0. However, the linear function tends to describe the slightly decreasing trend fairly good. The slope, and with that the effect of dependent scattering on

_{p}*g*, is more negative for wavelengths in the 800-1400 nm window with a minimum around 1000 nm. It was noticed by Mishchenko [35] that the effect of dependent scattering on

*g*is maximum if the wavelength is between 0.5 and 20 times larger than the radii of the scattering particles. As most of the fat globules in IL have a diameter smaller than 500 nm, with a mean particle radius around 50 nm [9], our results confirm these observations.

In order to define a general model for *g* as a function of both *λ* and *Φ _{p}*, an equation for the slope parameter

*k*had to be derived, while the wavelength dependency of

*k*cannot be described by a simple function. From Fig. 5(a), it can be observed that the

*k*-parameter decreases with increasing wavelength until it reaches a minimum around 1000 nm. For longer wavelengths it increases again until it tends to converge to a constant value. A fifth order polynomial was fitted (R

^{2}= 0.934) to the

*k*-data and plotted as a solid line in Fig. 5(a). The combination of the equation for

*k*and Eqs. (4) and (6) provides a general equation for

*g*as a function of both

*λ*and

*Φ*. This general model is presented as a surface in Fig. 6 and compared (R

_{p}^{2}= 0.992) with the measured

*g*-data (grey dots). It should, however, be emphasized that the general equations for

*µ*and

_{s}*g*were only validated against data for which

*λ*≥ 560 + 2200

*Φ*,

_{p}*λ*≤ 1850 nm and

*Φ*≤ 0.227 ml/ml, and that extrapolation of the models is not advised.

_{p}#### 3.3 Reduced scattering coefficient

If the general equations for *µ _{s}* and

*g*are combined with Eq. (5), a general model is obtained describing the reduced scattering coefficient (

*µ*) as a function of

_{s}’*λ*and

*Φ*. In Figs. 7(a) and 7(b), the effect of dependent scattering on the measured

_{p}*µ*values (cyan markers) for respectively 3 wavelengths below and 3 wavelengths above 1000 nm is presented. Moreover, the results of the general equation for

_{s}’*µ*at the 6 discrete wavelengths are plotted as solid lines. A very good agreement was found between the measured values and the model (all R

_{s}’^{2}> 0.993). At 635 nm, the model was in better agreement with the measured data than the results at 633 nm presented in [11] [Fig. 7(a)]. Moreover, a very good agreement was found with the empirical equations derived by Di Ninni

*et al.*[6] at 751 and 833 nm [Fig. 7(a) detail].

In Fig. 8, the general model for *µ _{s}’* (surface) is plotted together with the measured

*µ*data (grey dots). A very good correspondence (R

_{s}’^{2}= 0.999) can be observed. The model even performs well in the region for which no

*µ*and

_{s}*g*data could be extracted from

*µ*because of erroneous unscattered transmittance measurements (

_{s}’*λ*< 560 + 2200

*Φ*,

_{p}*λ*≥ 600 nm and

*Φ*≤ 0.227 ml/ml). This result strengthens the general models for

_{p}*µ*and

_{s}*g*, even in the regions where only limited data was available.

## 4. Conclusion

A set of 57 liquid optical phantoms was designed, prepared and optically characterized in order to study the effect of a dense distribution of scatterers on the scattering properties. Under the assumption of independent scattering, the bulk and reduced scattering coefficients follow a linear relation with the volume concentration for a specific type of scattering particles, while the anisotropy factor remains constant. However, for densely packed particles, the scattering events are correlated and the scattering properties deviate from those simple relations. From our study, it can be concluded that the bulk scattering coefficient as a function of the particle concentration for IL dilutions corresponds well with the packing factor proposed by Twersky in combination with a variable packing dimension. It was observed that in the case of constant particle sizes, the effect of dependent scattering is more significant for longer wavelengths. The developed general model for the bulk scattering coefficient describes the relation with wavelength and particle concentration very well (R^{2} = 0.999). From observations of anisotropy factors of IL dilutions, both in literature and measured, it was noticed that they follow a nearly linear decrease with increasing particle density. Moreover, the slope of the linear decrease was more steep for wavelengths around 1000 nm. This information was combined in a general model which describes the anisotropy factor as a function of the wavelength and particle concentration very well (R^{2} = 0.992). Additionally, by merging the models for the bulk scattering coefficient and the anisotropy factor, very accurate predictions of the reduced scattering coefficient values were obtained (R^{2} = 0.999). The equations or models given in this manuscript were derived following a semi-empirical approach and do not rely on a rigorous theory. Nevertheless, very good agreement was found with measured data and models derived by other researchers. Therefore, they provide a valuable tool for the calculation of the bulk scattering properties of IL dilutions, and the validation of theories that describe dependent scattering for IL-like systems.

## Acknowledgments

Ben Aernouts is funded as aspirant of the Research Foundation-Flanders (FWO-Flanders, grant 11A4813N). Rodrigo Watté and Robbe Van Beers are funded by the Institute for the Promotion of Innovation through Science and Technology in Flanders (IWT-Flanders, respectively grants 101552 and 131777). The authors gratefully acknowledge IWT-Flanders for the financial support through the GlucoSens project (SB-090053).

## References and links

**1. **V. V. Tuchin, *Tissue Optics: Light Scattering Methods and Instruments for Medical Diagnosis*, 2nd ed. (SPIE Press, 2007), p. 840.

**2. **A. Ishimaru, *Wave Propagation and Scattering in Random Media, Volume 1: Single Scattering and Transport Theory*, 1st ed. (Academic Press, 1978), p. 250.

**3. **B. W. Pogue and M. S. Patterson, “Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry,” J. Biomed. Opt. **11**(4), 041102 (2006). [CrossRef] [PubMed]

**4. **B. Cletus, R. Künnemeyer, P. Martinsen, and V. A. McGlone, “Temperature-dependent optical properties of Intralipid measured with frequency-domain photon-migration spectroscopy,” J. Biomed. Opt. **15**(1), 017003 (2010). [CrossRef] [PubMed]

**5. **P. D. Ninni, F. Martelli, and G. Zaccanti, “Intralipid: towards a diffusive reference standard for optical tissue phantoms,” Phys. Med. Biol. **56**(2), N21–N28 (2011). [CrossRef] [PubMed]

**6. **P. Di Ninni, F. Martelli, and G. Zaccanti, “Effect of dependent scattering on the optical properties of Intralipid tissue phantoms,” Biomed. Opt. Express **2**(8), 2265–2278 (2011). [CrossRef] [PubMed]

**7. **S. T. Flock, S. L. Jacques, B. C. Wilson, W. M. Star, and M. J. van Gemert, “Optical properties of Intralipid: a phantom medium for light propagation studies,” Lasers Surg. Med. **12**(5), 510–519 (1992). [CrossRef] [PubMed]

**8. **R. Michels, F. Foschum, and A. Kienle, “Optical properties of fat emulsions,” Opt. Express **16**(8), 5907–5925 (2008). [CrossRef] [PubMed]

**9. **H. J. van Staveren, C. J. Moes, J. van Marie, S. A. Prahl, and M. J. van Gemert, “Light scattering in Intralipid-10% in the wavelength range of 400-1100 nm,” Appl. Opt. **30**(31), 4507–4514 (1991). [CrossRef] [PubMed]

**10. **X. Wen, V. V. Tuchin, Q. Luo, and D. Zhu, “Controling the scattering of intralipid by using optical clearing agents,” Phys. Med. Biol. **54**(22), 6917–6930 (2009). [CrossRef] [PubMed]

**11. **G. Zaccanti, S. Del Bianco, and F. Martelli, “Measurements of optical properties of high-density media,” Appl. Opt. **42**(19), 4023–4030 (2003). [CrossRef] [PubMed]

**12. **B. Aernouts, E. Zamora-Rojas, R. Van Beers, R. Watté, L. Wang, M. Tsuta, J. Lammertyn, and W. Saeys, “Supercontinuum laser based optical characterization of Intralipid® phantoms in the 500-2250 nm range,” Opt. Express **21**(26), 32450–32467 (2013). [CrossRef] [PubMed]

**13. **P. I. Rowe, R. Künnemeyer, A. McGlone, S. Talele, P. Martinsen, and R. Oliver, “Thermal stability of intralipid optical phantoms,” Appl. Spectrosc. **67**(8), 993–996 (2013). [CrossRef] [PubMed]

**14. **A. Giusto, R. Saija, M. A. Iatì, P. Denti, F. Borghese, and O. I. Sindoni, “Optical properties of high-density dispersions of particles: application to intralipid solutions,” Appl. Opt. **42**(21), 4375–4380 (2003). [CrossRef] [PubMed]

**15. **G. Göbel, J. Kuhn, and J. Fricke, “Dependent scattering effects in latex-sphere suspensions and scattering powders,” Waves Random Media **5**(4), 413–426 (1995). [CrossRef]

**16. **M. I. Mishchenko, L. D. Travis, and A. A. Lacis, *Scattering, Absorption, and Emission of Light by Small Particles*, 3rd ed. (Cambridge University, 2006).

**17. **M. I. Mishchenko, L. Liu, D. W. Mackowski, B. Cairns, and G. Videen, “Multiple scattering by random particulate media: exact 3D results,” Opt. Express **15**(6), 2822–2836 (2007). [CrossRef] [PubMed]

**18. **D. W. Mackowski and M. I. Mishchenko, “Direct simulation of extinction in a slab of spherical particles,” J. Quantum Spectrosc. Radiat. Transfer **123**, 103–112 (2013). [CrossRef]

**19. **M. Lax, “Multiple scattering of waves. II. The effective field in dense systems,” Phys. Rev. **85**(4), 621–629 (1952). [CrossRef]

**20. **V. P. Dick, “Applicability limits of Beer’s law for dispersion media with a high concentration of particles,” Appl. Opt. **37**(21), 4998–5004 (1998). [CrossRef] [PubMed]

**21. **D. W. Mackowski and M. I. Mishchenko, “Calculation of the T matrix and the scattering matrix for ensembles of spheres,” J. Opt. Soc. Am. A **13**(11), 2266–2278 (1996). [CrossRef]

**22. **A. Ishimaru, *Wave Propagation and Scattering in Random Media, Volume 2: Multiple Scattering, Turbulence, Rough Surfaces, and Remote Sensing*, 1st ed. (Academic, 1978), p. 319.

**23. **M. I. Mishchenko, “Directional radiometry and radiative transfer: A new paradigm,” J. Quantum Spectrosc. Radiat. Transfer **112**(13), 2079–2094 (2011). [CrossRef]

**24. **M. I. Mishchenko, D. H. Goldstein, J. Chowdhary, and A. Lompado, “Radiative transfer theory verified by controlled laboratory experiments,” Opt. Lett. **38**(18), 3522–3525 (2013). [CrossRef] [PubMed]

**25. **R. West, D. Gibbs, L. Tsang, and K. Fung, “Comparison of optical scattering experiments and the quasi-crystalline approximation for dense media,” J. Opt. Soc. Am. A **11**(6), 1854–1858 (1994). [CrossRef]

**26. **L. Tsang, J. Kong, and T. Habashy, “Multiple scattering of acoustic waves by random distribution of discrete spherical scatterers with the quasicrystalline and Percus–Yevick approximation,” J. Acoust. Soc. Am. **71**(3), 552–558 (1982). [CrossRef]

**27. **A. Ishimaru and Y. Kuga, “Attenuation constant of a coherent field in a dense distribution of particles,” J. Opt. Soc. Am. **72**(10), 1317–1320 (1982). [CrossRef]

**28. **V. Twersky, “Low-frequency scattering by correlated distributions of randomly oriented particles,” J. Acoust. Soc. Am. **81**(5), 1609–1618 (1987). [CrossRef]

**29. **V. Twersky, “Acoustic bulk parameters in distributions of pair-correlated scatterers,” J. Acoust. Soc. Am. **64**(6), 1710–1719 (1978). [CrossRef]

**30. **P. A. Bascom and R. S. Cobbold, “On a fractal packing approach for understanding ultrasonic backscattering from blood,” J. Acoust. Soc. Am. **98**(6), 3040–3049 (1995). [CrossRef] [PubMed]

**31. **J. M. Schmitt and G. Kumar, “Optical scattering properties of soft tissue: a discrete particle model,” Appl. Opt. **37**(13), 2788–2797 (1998). [CrossRef] [PubMed]

**32. **N. E. Berger, R. J. Lucas, and V. Twersky, “Polydisperse scattering theory and comparisons with data for red blood cells,” J. Acoust. Soc. Am. **89**(3), 1394–1401 (1991). [CrossRef] [PubMed]

**33. **V. Twersky, “Low-frequency scattering by mixtures of correlated nonspherical particles,” J. Acoust. Soc. Am. **84**(1), 409–415 (1988). [CrossRef]

**34. **A. Bashkatov, E. Genina, V. I. Kochubey, and V. Tuchin, “Effects of scattering particles concentration on light propagation through turbid media,” Proc. SPIE **3917**, 256–263 (2000). [CrossRef]

**35. **M. Mishchenko, “Asymmetry parameters of the phase function for densely packed scattering grains,” J. Quantum Spectrosc. Radiat. Transfer **52**(1), 95–110 (1994). [CrossRef]

**36. **S. N. Thennadil, H. Martens, and A. Kohler, “Physics-based multiplicative scatter correction approaches for improving the performance of calibration models,” Appl. Spectrosc. **60**(3), 315–321 (2006). [CrossRef] [PubMed]

**37. **E. Zamora-Rojas, B. Aernouts, A. Garrido-Varo, D. Pérez-Marín, J. E. Guerrero-Ginel, and W. Saeys, “Double integrating sphere measurements for estimating optical properties of pig subcutaneous adipose tissue,” Innov. Food Sci. Emerg. Technol. **19**, 218–226 (2013). [CrossRef]

**38. **R. Watté, N. N. Do Trong, B. Aernouts, C. Erkinbaev, J. De Baerdemaeker, B. M. Nicolaï, and W. Saeys, “Metamodeling approach for efficient estimation of optical properties of turbid media from spatially resolved diffuse reflectance measurements,” Opt. Express **21**(26), 32630–32642 (2013). [CrossRef] [PubMed]

**39. **S. A. Prahl, “Everything I think you should know about Inverse Adding-Doubling,” http://omlc.ogi.edu/software/iad/iad-3-9-10.zip.

**40. **G. M. Hale and M. R. Querry, “Optical constants of water in the 200-nm to 200-μm wavelength region,” Appl. Opt. **12**(3), 555–563 (1973). [CrossRef] [PubMed]

**41. **H. C. van de Hulst, *Light Scattering by Small Particles* (John Wiley, 1957), p. 470.