## Abstract

In this paper, a general method to calibrate the absorption coefficient of an absorber and the reduced scattering coefficient of a liquid diffusive medium, based on time-resolved measurements, is reported. An exhaustive analysis of the error sources affecting the estimation is also performed. The method has been applied with a state-of-the-art time-resolved instrumentation to determine the intrinsic absorption coefficient of Indian ink and the reduced scattering coefficient of Intralipid-20%, with a standard error smaller than 1% and 2%, respectively. Finally, the results have been compared to those retrieved for the same compounds by applying a continuous wave method recently published, obtaining an agreement within the error bars. This fact represents a cross validation of the two independent calibration methods.

©2007 Optical Society of America

## 1 Introduction

In the last years, the Optical Community has done a lot of work to the aim of validating the use of near infrared (NIR) light for biological tissue probing and disease diagnosis. This approach has been stimulated by its complete non-invasiveness and relative cheapness with respect to other medical investigation techniques.

Applications involving NIR light are based on the property of wavelengths in this range to deeply propagate inside biological tissues (therapeutic windows), allowing to investigate large tissue volumes and deep biological structures (breast, brain, bone, skin, etc.). To this aim, different instrumental set-ups have been developed, namely Continuous Wave (CW), time and frequency domain systems. As for data interpretation, the main information recovered by these techniques is summarized by the optical properties of the explored tissues. Due to the heterogeneity of various biological tissues, also a great effort has been devoted to the development of suitable theoretical models for data analysis.

In this framework, we are observing the moving of this technique to its extreme limits of applicability, trying to extract reliable information from any minimal recovered difference in the optical properties of the investigated tissues, and to assign them to some tissue alterations, both of physiological or pathological origin. To get this task, however, we can no more avoid to evaluate the instrument capability of accurately estimating the optical properties of explored tissues. Even if some interesting results can be assessed by performing differential measurements and analysis, it is out of doubts that, for keeping and improving the skill reached by this field of research, many instruments in a near future more and more will base their performances on the evaluation of the absolute values of the optical properties of probed samples.

To accomplish this need, many methodologies based on measurements both in the time, frequency, and CW domain, have been reported that, in principle, can give accurate measurements of *μ*′_{s} and *μ _{a}* coefficients, and many tissue-like phantoms have been proposed to assess the performance of NIR instrumentation [1]. However, by looking at recently published papers, it seems difficult to find measurements with an error less than 10%, even when they are carried out on homogeneous media in the most favorable experimental conditions. On this line, there are the results of recent papers [2, 3] showing that the calibration of the optical properties of Intralipid (one of the most widespread scatterer compound used to produce calibrated diffusive phantoms) is subjected to an error never less than 8%.

We stress that the availability of calibrated phantoms is of great help for measuring precise absolute values of the optical properties. In order to prepare a calibrated liquid phantom it is important to have both a method to calibrate the reduced scattering coefficient of a liquid suspension of a diffusive medium and a method to calibrate the absorption coefficient of an absorber. For this purpose, we present here a method - for time-resolved instruments - for calibrating the reduced scattering coefficient of a liquid suspension of diffusive medium, namely Intralipid-20%, and the absorption coefficient of an absorber (Indian ink) with a standard error smaller than 2%. In the companion paper recently published [4], a calibration method for CW instrumentations has been proposed and similar performances have been demonstrated. In the procedures proposed here, the absorption coefficient is retrieved by the logarithmic ratio of two measurements at different concentrations of the absorber [5, 6], while the reduced scattering coefficient is retrieved by the logarithmic ratio of two measurements at different source-detector distances [7, 8]. Both inversion procedures used are based on linear regressions. The simplicity of the analysis resulting in such a way has its drawback in the sophistication and complexity of the instrumental set-up. In particular, to get the performances stated above on the error, the time resolution has to be pushed to the state-of-the-art limit. On the other hand, by relaxing the instrumental performances, we get less accurate results.

The proposed method to calibrate absorber and scatterer compounds by means of time-resolved diffusive measurements is presented in Sec. 2. In Sec. 3 a complete analysis of error sources is reported and their entity evaluated. In Sec. 4 we have calibrated with the proposed method Indian ink, as an absorber, and Intralipid-20%, as a scatterer, and compared the results obtained to those retrieved in [4], where the same diffuser and absorber have been calibrated. Finally, some conclusions are reported in Sec. 5.

## 2. The method

In this section, we present the method proposed in order to calibrate absorbers and scatterers by means of time-resolved diffusive measurements.

#### 2.1 The absorber

The light propagation in a scattering and absorbing medium can be described by considering the Radiative Transfer Equation (RTE), that expresses the energy balance inside a volume element of the medium. The time-dependent form of RTE can be cast as follows [9]:

where *I*(*x*⃗,*t*,*s*̑) is the radiance in the direction *s*̑ ; *v* is the speed of light in the medium; *μ _{s}* and

*μ*are the scattering and absorption coefficient of the medium, respectively;

_{a}*p*(

*s*̑,

*s*̑′) is the scattering function representing the probability for a photon traveling in the direction

*s*̑′ to be scattered in the direction

*s*̑ ;

*q*(

*x*⃗,

*t*,

*s*̑) is the source term, representing the light power per unit of volume and solid angle.

If the time dependence of the source term is a Dirac delta function, Eq. (1) has the important property that, if *I*
_{0}(*x*⃗,*t*,*s*̑) is the solution corresponding to a non-absorbing medium (*μ _{a}*=0),

*I*(

*x*⃗,

*t*,

*s*̑) =

*I*

_{0}(

*x*⃗,

*t*,

*s*̑)exp(-

*μ*) is the solution when a uniform absorption coefficient

_{a}vt*μ*is present, that is the absorption coefficient appears only in the exponential factor. This fact can be exploited to calibrate an absorber.

_{a}Let’s suppose to prepare a liquid diffusive medium with defined but not known values of *μ _{s}* and

*μ*

_{a0}; the detected power in a defined point

*x*⃗

_{d}of the medium can be written as:

where *A*
_{0} accounts for the source power and the detection efficiency, while in the term *S* are gathered all the geometry factors, the boundary conditions and the scattering properties of the medium. Now, if we add to the medium a defined quantity of the absorber to be calibrated, the new power detected at the same position is:

where we have assumed that the added absorber does not change the scattering properties of the medium and that the new absorption coefficient is *μ _{a}*. We have also allowed that the source power and/or the detection efficiency can change from the previous measure.

Then, if we consider the ratio of Eq. (2) and (3), the term *S* simplifies and, by taking the natural logarithm of both sides, we recover a linear relation between the measured power and time:

where Δ*μ _{a}* =

*μ*-

_{a}*μ*

_{a0}is the increase of the absorption coefficient due to the added absorber. Now, we can extract the value of Δ

*μ*from the slope of the linear regression between the independent variables

_{a}*X*=

*t*and

*Y*≡ ln(

*M*

_{0}/

*M*).

Then, if we assume to be in a linear regime, by repeating the measure above with different quantities of the added absorber, we can calculate the intrinsic absorption of the absorber *ε _{a}*, that is the absorption coefficient

*per*unit of absorber concentration

*ρ*, expressed as the volume of the absorber over the volume of the diffusing mudium:

_{a}#### 2.2 The scattering compound

In order to calibrate a scattering compound, we can not exploit a simple dependence from *μ _{s}* in the RTE, as we did for the absorption coefficient. However, it is well known that, if some simplifying assumptions are considered in Eq. (1) [9], the light propagation in the diffusive medium can be described by the Diffusion Equation (DE) for the average diffuse intensity

*U*(

*x*⃗,

*t*):

where *D* = 1/(3*μ*′_{s}) is the diffusion coefficient, *μ*′_{s} = *μ*
_{s}(1-*g*) being the reduced scattering coefficient and *g* the anisotropy factor, and *Q*(*x*⃗,*t*) is the isotropic source term.

If we assume *Q*(*x*⃗,*t*) = 1/4*π*
*δ*(*x*⃗-*x*⃗_{s})*δ*(*t*), representing a point-like isotropic source located at *x*⃗_{s}, and consider an infinite slab of diffusing medium with optical properties *μ _{a}* and

*μ*′

_{s}, according to Eq. (6), the measured power in both the reflectance and the transmittance geometry is given by [7, 8]:

where *A* accomplishes for the source power and the detection efficiency, while the factor *S* accounts for the boundary conditions and depends on only the slab thickness *d* and the scattering properties of the medium. Finally, *r* is the distance in the slab plane between the injection and detection points.

Since *r* appears only in the exponential factor, if we perform two measurements at two different source-detector lateral distances (*r*
_{0} and *r*) and take the natural logarithm of the ratio of the measured power in the two conditions *M*
_{0} and *M*, respectively, we get a linear relation between such logarithm and the inverse of time:

where we have assumed that the source power and/or the detection efficiency can change in the two measurements. From the slope of the linear regression between the independent variable *X* = 1/*t* and *Y* ≡ ln(*M*
_{0}/*M*) one can determine the reduced scattering coefficient *μ*′_{s}. We note that expression (7) and, then, (8) are strictly valid only for point-like source and detector, located at *x*⃗_{s} and *x*⃗_{d} , respectively. For finite dimensions of source and/or detector some distortions of the linear dependence foreseen in (8) can be expected.

Next, one can characterize the scatterer by giving the intrinsic reduced scattering coefficient *ε*′_{s}, that is the reduced scattering coefficient *per* unit of the scatterer concentration *ρ _{s}*. In a linear regime, one has simply:

where *ρ _{s}* is expressed as the volume of the scatterer over the volume of the diffusing suspension.

As a concluding remark on the proposed method, we point out that both the absorber and the scatterer calibrations finally refer to the intrinsic values of the optical properties. As a matter of fact, the intrinsic absorption coefficient and the intrinsic reduced scattering coefficient are the values that are actually needed to prepare calibrated phantoms.

## 3 Error analysis

In this section we describe the possible error sources that can affect the previous estimation of the absorption and scattering coefficients. Some of the error sources are the same for the two coefficients, while others are peculiar to only one of them. In order to give a clear description, we have preferred to itemize the error sources and for each of them to describe the effect they have on the two coefficients.

#### 3.1 Time scale

Because the method relies on time-resolved measurements, it is important to consider the time characteristics of the system one uses. The errors that can affect the time measure are of two main kind: i) the origin of time *t*
_{0} can fluctuate from one measure to another; ii) the time scale can be expanded or shrunk by a factor *α*.

As for point i), such kind of error does not affect at all the estimation of the absorption coefficient. From Eqs. (2) and (3), if the origin of time is *t*
_{0} and *t*′_{0} for the measurements *M*
_{0} and *M*, respectively, the expression of these two quantities becomes:

$$M({\overrightarrow{x}}_{d},t)=AS({\overrightarrow{x}}_{d},t,{\mu}_{s})\mathrm{exp}\left[-{\mu}_{a}v\left(t-{\mathrm{t\prime}}_{0}\right)\right].$$

and the natural logarithm of their ratio is:

Then, the slope of the linear regression between *X* ≡ *t* and *Y* ≡ ln(*M*
_{0}/*M*) is not affected by the two time shifts *t*
_{0} and *t*′_{0}.

Now we consider the reduced scattering coefficient. If the origin of time is *t*
_{0} and *t*′_{0} for the measurements *M*
_{0} and *M*, respectively, the Eq. (8) becomes:

where we have retained only the first order of the quantities *t*
_{0}/*t* and *t*′_{0}/*t*. This assumption is reasonable because one expects only small fluctuations in time. From Eq. (12), it is possible to foresee a quadratic deviation from the linear dependence of *Y* ≡ ln(*M*
_{0}/*M*) on *X* ≡ 1/*t*. This deviation is of the order of |*t*
_{0}|/*t* (or |*t*′_{0}|/*t*) and obviously depends on the time range where the linear regression is performed. Then, one simple way for reducing the uncertainty due to a systematic time shift *t*
_{0} is to avoid in the retrieval procedure the use of the first part of the temporal range. As an example, we consider a slab 20 mm thick, *μ*′_{s} =1.5 mm^{-1} and source-detector lateral distances *r*
_{0} = 0 and *r*
_{1} = 25 mm: for a *t*
_{0} = ±20 ps, we obtain a systematic error on the retrieved *μ*′_{s} that, assuming a time range for the linear regression with a lowest time of 1000 ps and a maximum time always greater than 3000 ps, is less than 3%, while for a *t*
_{0} = ±10 ps the error is less than 1.5%.

As for point ii), the time scale of the measurements determined by the instrumentation could suffer from some uncertainty (*e*.*g*. with 0.1% error, a time interval of 2 ns could results in a measure of 2.002 ns). We took into account this fact by putting such a kind of uncertainty into a scale factor *α* so that the time *t* of the measurements can be expressed as *t* = *α*
*t*̃, where *t*̃ is the real time. Then, because in the proposed method we are interested to the slope of a linear regression where the independent variable is proportional to time (*X* ≡ *t* and *X* ≡ 1/*t* for absorption and scattering, respectively), *α* enters this slope as a factor and then its eventual relative error must be added to the relative error of the slope.

#### 3.2 Instrument Response Function

In real measurements we have to face the problem of the Instrument Response Function (IRF). The linear expression obtained in Eqs. (4) and (8) are relative to the case of an infinitely short injected light pulse (Dirac delta function in time). Actually, the IRF of any instrument has a finite width and specific shape that has the effect to deform the time-resolved curves. Then, we have studied how different IRFs affect the estimation of the absorption and scattering coefficients obtained by applying the method reported above.

To this aim, we simulated the measured light distribution *T*(*r*,*t*;*μ _{a}*,

*μ*′

_{s}) considering the theoretical time-resolved curve, solution of the DE for a slab in transmission geometry with extrapolated boundary conditions, convolved with a given IRF.

In particular, for exploring the effect of the IRF on the absorption coefficient, we have calculated the error, as a function of time *t*, in the estimation of Δ*μ _{a}* by applying the method described in Sect. 1.1. Due to the fact that we have the same amplitude factor

*A*for all curves, the estimated absorption increment (Δ

*μ*)

_{a}^{est}(

*t*) results:

The estimated absorption increment (Δ*μ _{a}*)

^{est}(

*t*) differs from the true Δ

*μ*of a quantity that depends on the IRF, on the geometric characteristic and the optical properties of the medium,

_{a}*μ*

_{a0}and

*μ*′

_{s}, and also on the same Δ

*μ*. In order to view the behavior of the uncertainty on (Δ

_{a}*μ*)

_{a}^{est}(

*t*) due to the IRF, Figs. 1(b) and 1(c) show the ratio (Δ

*μ*)

_{a}^{est}/Δ

*μ*as a function of time. In these figures, a couple of IRFs are presented, for different values of Δ

_{a}*μ*, all referred to the light transmitted through a non-absorbing slab, 40 mm thick, with

_{a}*μ*′

_{s}= 0.5cm

^{-1}. The first IRF considered (see curve IRF_1 in Fig. 1(a)) is related to the instrument we have employed for the experimental measurements (see Sec. 4.1) and has a full-width at half-maximum (FWHM) of about 50 ps. The second IRF considered (see curve IRF_2 in Fig. 1(a)) is related to the multi-channel time-resolved tissue oximeter described in [10] and is characterized by a FWHM of about 500 ps. The investigated values of Δ

*μ*are 0.001, 0.004 and 0.008 mm

_{a}^{-1}. The figure shows that the distortion introduced by the IRF is strongly depending on the FWHM of the IRF and on the Δ

*μ*used for the calibration procedure. Furthermore, excluding the first 1000 ps, the distortion introduced by the convolution mainly results in a systematic underestimation of Δ

_{a}*μ*. As an example, for the IRF_1 and Δ

_{a}*μ*= 0.001mm

_{a}^{-1}we have that (Δ

*μ*)

_{a}^{est}is affected by an error smaller than 1.5%. This error is also decreasing as the

*μ*′

_{s}of the medium is increased.

As for the scattering coefficient, for applying the method described in Sect. 1.2, we considered simulated measurements for two source-detector lateral distances: *r*
_{0} = 0 and *r*
_{1} = 25 mm. Two examples of simulated curves for such lateral distances are reported in Fig. 1(d). From Eq. (8), considering that the amplitude factor is the same, we can calculate the estimated reduced scattering coefficient (*μ*′_{s})^{est} as a function of time as follows:

Also the estimated reduced scattering coefficient (*μ*′_{s})^{est} differs from the true *μ*′_{s} of a quantity that depends on the IRF, on the geometric characteristics and optical properties of the medium, and on the two distances *r*
_{0} and *r*
_{1}. In order to evaluate the behavior of the uncertainty affecting (*μ*′_{s})^{est}(*t*) due to the instrumental response function, in Figs. 1(e) and 1(f) it is shown the ratio (*μ*′_{s})^{est}/*μ*′_{s} as a function of time. Three different values of *μ*′_{s} are considered (*μ*′_{s} = 1, 1.5 and 2 mm^{-1}) all related to the light transmitted through a 20 mm thick slab with *μ*
_{a} = 0. The figures show that the error introduced by the IRF on the assessed (*μ*′_{s})^{est}(*t*) is strongly depending on the FWHM of the IRF and on the *μ*′_{s} used for the calibration procedure. Furthermore, excluding the first 1000 ps, the distortion due to the convolution prevalently results in a systematic overestimation of *μ*′_{s}. As an example, for the IRF_1 and *μ*′_{s} = 2 mm^{-1} we have that the error affecting (*μ*′_{s})^{est} remains less than 1.1%. This error increases as the absorption of the medium is increased.

In general, we can conclude that the key factors for keeping the systematic distortion due to the IRF sufficiently low are the optical parameters and the thickness of the diffusing medium: one must adjust these parameters in such a way to increase as much as possible the temporal width of the diffused curve with respect to the FWHM of the IRF. In the experimental results that will be shown in Sec. 4 this uncertainty will be evaluated and added to all the other uncertainties affecting the measurements.

#### 3.3 Model adopted for light propagation

Generally speaking, also the approximations adopted in the description of the light propagation in diffusive media are sources of errors in the estimation of the optical properties. Because the method described in Sect 2.1 for retrieving the absorption properties of the medium relies on a general property of the RTE, we do not expect any influence on the estimation of Δ*μ*
_{a} as long as the RTE gives an accurate description of the stream of photons flowing inside the medium.

On the contrary, the method described in Sect. 2.2 for estimating the scattering properties of the medium assumes that the RTE can be approximated by the DE. Furthermore, when we use Eq. (8), we also assume that source and detector can be considered point-like. These approximations imply some constraints on the optical properties of the diffusive medium on the geometry of the measurements and on the time range of the TPSF. In order to have a quantitative evaluation of the influence of the diffusion approximation, we have performed some Monte Carlo simulations of the light propagation and compared the resulting TPSFs with those calculated from the DE, for different values of the optical parameters. In particular, a diffusing slab of thickness 20 mm, *μ*
_{a}=0 and *μ*′_{s} = 1.5mm^{-1} has been considered. Moreover, the source and receiver have been assumed with a dimension of 3 mm. The method described in Section 2.2 has been used for retrieving the *μ*′_{s} of the medium, considering the two distances *r*
_{0} = 0 and *r*
_{1} = 25 mm. Provided that all the times with *t* < 2.5*t*
_{B} (*t*
_{B} being the ballistic time between source and receiver) were disregarded from the retrieval procedure, the systematic error due to the intrinsic approximations of the DE shows to be always less than 0.2%. By changing some of the parameters, also the error changes. In particular, if the source-detector lateral distance is decreased, the error increases: for instance, with *r*
_{1} = 10 mm we
have an error always less than 1.2%. On the other hand, increasing *μ*′_{s} and keeping all the others parameters unchanged the error decreases.

#### 3.4 Source-detector lateral distance

The method suggested for the determination of the increase in the absorption coefficient due to an addition of absorber is really independent on geometrical aspects. For this reason, the knowledge of the source-detector lateral distance is not necessary.

On the contrary, the method for the scattering properties is based on a comparison between two measurements at two different values of the source-detector lateral distance *r* (see Eq. (8)). Then, errors in the determination of such distances contribute straightforwardly to the indetermination of the scattering coefficient. More precisely, from Eq. (8) we have:

where *s* is the slope of the linear relation appearing in the equation. Then, if we call σ_{r} the standard deviation for *r*, by using the formulas for the propagation of stochastic errors, the contribution to the relative error of the scattering coefficient *ε*
_{μ′s} due to the uncertainty of *r* is given by:

#### 3.5 Preparation of the liquid diffusive media

In the calibration method proposed above, different diffusive suspensions must be prepared. In particular, in order to estimate the intrinsic absorption and reduced scattering coefficients of the used absorber and scatterer, one need to know the concentrations of such compounds in the prepared suspensions. From Eqs. (5) and (9) we note that an uncertainty affecting the concentrations *ρ _{a}* and

*ρ*contribute accordingly to the relative errors of

_{s}*ε*and

_{a}*ε*′

_{s}, respectively.

#### 3.6 Concluding remarks on error analysis

In order to make clearer the role of possible contributions to the indetermination of the optical properties determined by means of the method proposed in this work, we sum up here the considerations stated above for the absorption and the scattering coefficient separately.

As for the absorption, the relative error *ε*
_{Δμa} affecting the increase in the absorption
coefficient Δ*μ*
_{a} due to all the possible stochastic error sources can be estimated as (see Sec. 3.1):

where *ε _{s}* and

*ε*are the relative fluctuations of the slope

_{α}*s*of the linear relationship present in Eq. (4) and of the time-stretching factor

*α*introduced in Sec. 3.1, respectively. Moreover, to this value we have to add the systematic error due to the IRF of the measurement system (see Sec. 3.2).

As for the scattering property, the relative error of the scattering coefficient *ε*
_{μ′s} due to all the possible stochastic error sources can be estimated as (see Sec. 3.1, Eqs. (15) and (16)):

where *ε*
_{t0} is the relative error due to the stochastic shift of the time origin *t*
_{0}. As before, to
expression (18) we have to add the systematic errors introduced by the IRF and the diffusion model approximation (see Secs. 3.2 and 3.3).

Finally, we note that the uncertainties affecting the intrinsic absorption coefficient *ε _{a}* and the reduced scattering coefficient

*ε*′

_{s}of the absorber and the scatterer compounds, respectively, can be smaller than those of Δ

*μ*and

_{a}*μ*′

_{s}calculated in Eqs. (17) and (18), if they are retrieved from a series of measurements by varying the compound concentrations, according to Eqs. (5) and (9).

## 4 Experimental results

Following the indication suggested by the error analysis performed above, we set-up a time-resolved system for the calibration as accurate as possible of the absorber and scatterer compounds already calibrated with a CW system at the wavelength of 750 nm, as described in the companion paper [4].

#### 4.1 System set-up

The system we used is an adaptation of an automated system for time-resolved spectroscopy measurements described in previous publications [11]. It allows fully spectroscopic measurements in the range from 700 to 1000 nm, exploiting an actively mode-locked picosecond Ti:sapphire laser, as a source, and a double microchannel plate photomultiplier tube (MCP) for detection. The temporal resolution is obtained by means of the time-correlated single-photon counting technique, implemented by a personal computer board. To our aim, the Ti:sapphire laser emission wavelength was tuned and optimized to 750 nm, in order to compare the calibration measurements with those obtained with the CW system. The output beam is coupled to a 50/125 *μ*m graded-index telecommunication optical fiber, that guarantees minimal temporal dispersion, and delivered to the sample. It consists of a black modular tank with lateral dimension of 300×300 mm, while the thickness can be adjusted in a wide range (from few millimeters to 150 mm). The front wall of the tank has a series of transparent windows of 3 mm in diameter, spaced by steps of 5 mm, that allow off-axis measurements. The back wall has a similar transparent window that has been located just in front of the collecting lens system of the MCP. In this way, we have reduced to the minimal extent the temporal width of the IRF of the system, obtaining a full-width at half-maximum of about 50 ps. Furthermore, we adjusted the thickness of the tank, the optical parameters of the diffusing medium to be measured and the acquisition times, in order to obtain time-resolved curves filling as much as possible the time window, that we fix to about 8 ns, with the dynamical range limited only by the background. In this way we can reduce significantly the errors due to the IRF. Fig. 2 shows, as an example of the quality of the data, the IRF and a time-resolved diffused curve, recorded by our system. We can see that on the tail of the curve the measurement extends over about three decades. From this figure we can also evaluate the entity of the background present in our measurements, that is always subtracted in the analysis procedure.

For what concerns the time stability of the system, during the measurements we often recorded the IRF in order to monitor the time shift of our system. On the other hand, the error on the time-scale factor *α*, considered in Sec. 3.2, for our system consists in the indetermination of the conversion factor from MCA channel to picoseconds. For integrated systems this factor is usually calibrated by the factory: for our system, *ε _{α}* can be assumed to be a fraction of a percent [12].

For the preparation of the diffusive suspensions, we used the same batches of Intralipid-20% and Indian ink calibrated in [4]. We determined the concentration of the Intralipid in the diffusive suspension by means of a precision balance (error of ±1 mg). Due to the quantities of Intralipid considered, the uncertainty affecting its concentration *ρ*
_{sil} is of the order of
2 ∙ 10^{-4}. On the other hand, we pre-diluted in distilled water the Indian ink, in order to reduce the uncertainty affecting its concentration *ρ*
_{aink} in the diffusive suspensions. In this way, the relative error of *ρ*
_{aink} can be estimated as 2 ∙ 10^{-3}.

#### 4.2 Calibration of the absorber

Here we apply the method described in Sec. 2.1 in order to calibrate the absorption properties of the same batch of Indian ink considered in [4].

Firstly, we have prepared a diffusing suspension with a reduced scattering coefficient
*μ*′_{s} = 0.64 mm^{-1}. This will be our starting diffusive medium to which we will add different quantities of the pre-diluted absorber. We put this diffusing suspension in a tank of thickness *d* = 40.8 mm and measure the transmitted on-axis time-resolved curve. In Fig. 2 is reported such a curve together with the IRF.

In order to evaluate the error caused by the IRF for this specific case, we have performed numerical simulations as described in Sec. 3.2, for increasing values of Δ*μ _{a}*. The linear fitting method described in Sec. 2.1 was applied for different time ranges. We decided to fix the time ranges considering, for each curve, temporal points where the counts are a fixed percentage of the counts of the curve peak. In this way the time ranges adapt themselves to different curve shapes. From what we concluded in Sec. 3.2, we have to consider both extremes of the time ranges on the tail of the curves. For the ending time we chose the time corresponding to 0.5% of the curve peak. In this way we are reasonably away from the noise of the final part of the curve (see Fig. 2). Moreover, we varied the initial time taking the time corresponding to different percentages (namely, from 10% to 1%) of the counts of the curve peak. In Fig. 3 we report the relative error of the absorption coefficient, calculated for different values of the absorption increment, as a function of the initial time. It results that a reasonable choice for the percentage corresponding to the initial time is 5%: indeed, for this choice, if we limit to absorption increase of about 0.006 mm

^{-1}, we keep the relative error introduced by the convolution well below 1%.

Then, we performed 10 different addings of the absorbing suspension, measuring each time the transmitted time-resolved diffused curve. Next, we applied the linear method and calculated for each value of the added absorbing suspension the corresponding value of Δ*μ _{a}*. In Fig. 4 we have reported the value of Δ

*μ*as a function of the absorber concentration. Furthermore, we have evaluated the relative error of the slope of the linear regression expressed in Eq. (4), resulting of the order of 1%. For what we stated about the error due to the convolution and the time-scale factor, this is the major source of error for the first five addings, corresponding to a maximum absorption increase Δ

_{a}*μ*= 0.0061

_{a}*mm*

^{-1}. Then, by considering that independent relative errors have to be summed squared (see Eq. (17)), for these first five measurements the relative error affecting the estimated value of Δ

*μ*can be assumed to be 1%. For the higher values of absorber concentrations, corresponding to larger values of absorption, the systematic error introduced by IRF becomes prevalent, resulting in an underestimation of the absorption coefficient (compare Fig. 3).

_{a}Finally, from the previous set of measurements, we can calculate the intrinsic absorption of the absorber. According to our considerations on the errors, we performed the linear regression reported in Eq. (5), considering only the first five measurements. We are led to the following best value for the intrinsic absorption *ε*
_{aink} of Indian ink:

Here, we have estimated a relative error of 0.5% as a result of the linear fitting procedure and concentration uncertainties. It is worth noting how the method proposed allows a very small uncertainty in the intrinsic absorption determination.

#### 4.3 Calibration of the scatterer

In this section we describe the results obtained with our system in the determination of the scattering properties of the batch of Intralipid 20% already calibrated with a CW laser system, as described in the companion paper [4].

We used a black tank of thickness *d* = 20.4 mm and performed transmittance measurements for 6 different off-axis values *r*: 0 (on-axis), 5, 10, 15, 20, 25 mm. Indeed the method for the scatterer calibration relies on measurements at different source-detector lateral distances. Then, we prepared four diffusing suspensions corresponding to the following values for reduced scattering coefficient: 0.97, 1.4, 1.8, 2.2 mm^{-1}, as one can establish adopting the calibration reported in [4]. For each suspension, we applied the linear method described in Sec. 2.2, comparing each curve obtained for *r* = 5, 10, 15, 20, 25 mm with that measured on-axis (*i*.*e*. we assumed *r*
_{0} = 0). We note that we reduced the thickness of the tank (see Sec. 4.2), in order to enhance the differences between the on-axis and off-axis measurements. Then, by applying the linear regression reported in Eq. (8), we get 5 estimations for the reduced scattering coefficient *μ*′_{s}.

As for the estimation of the errors, we note that in our measurements the biggest source of indetermination is given by the uncertainty affecting the lateral distance *r*: we can assume that in centering the 3 mm diameter transparent window with an about 2 mm wide beam, we make a maximum error of σ_{r} =0.5 mm. According to Eq. (16), this means that the contribution to the relative error of *μ*′_{s} ranges from 20%, for *r* = 5 mm, to 4%, for *r* = 25 mm. As for the other sources of errors, the contribution *ε _{s}* due to the indetermination of the slope of the linear regression is always lower than 1%. Furthermore, the error due to the convolution with the IRF can be evaluated by means of simulations. To this aim, we considered different time ranges for the linear regression in Eq. (8), expressed in percentage of the curve peak, considering for the raising edge the curve at

*r*> 0 and for the tail the curve at

*r*= 0, and calculated the relative error of the estimated reduced scattering coefficient. In Fig. 5 are reported the plots resulting from such a calculation as a function of the source-detector lateral distance

*r*, for the values of reduced scattering coefficients used in the measurements. By inspecting this figure, we decided to use as linear fitting range the time interval from 10% in the rising edge of the curve to 1% in the tail: with this choice, we have a relative error less that 2%, at least for the three largest lateral distances (compare Fig. 5(b)). As for the time shift affecting our measurements, we monitored it by recording the IRF after every set of about 5 measures taken over a period of about 10 minutes. The resulting indeterminacy on the time origin for each independent acquisition was about ±2 ps. Since the time range we consider for the linear regression has an initial time of some hundreds of picoseconds, also the

*ε*

_{t0}contribution to the total error can be assumed well lower than 1%. Finally, considering the time ranges adopted, the optical and geometrical properties of the measured suspensions and the simulations reported in Sec. 3.3, the error due to the diffusion approximation is estimable on the order of a fraction of percent.

According to this analysis of the errors affecting the estimation of the reduced scattering coefficient, we decided for each diffusing suspension to retain only the estimations of *μ*′_{s} emerging from the largest lateral distances: *r* = 15, 20, 25 mm. Only for these measurements, in fact, the error can be assumed lower than 10%. In order to calculate the errors affecting these three estimations of *μ*′_{s}, we summed up all the relative errors squared, as reported in Eq. (18), noting that the only relevant contributions come from the lateral distances *r* and the convolution. Finally, we ended up with the best value for *μ*′_{s}, by performing a weighted average of the three estimations of *μ*′_{s} for each source-detector lateral distance [13]. The resulting relative error affecting the best value is of 3%.

In Fig. 6(a) the natural logarithm of the ratios of the measured curves, for *r*
_{0} = 0 and *r* = 25mm, *Y* = ln(*M*
_{0}/*M*), as function of *X* ≡ 1/*t* and resulting linear best fits (line) according to Eq. (8), are reported as examples for four values of Intralipid-20% concentration. Furthermore, in Fig. 6(b) the four best estimations of *μ*′_{s} are plotted as a function of the Intralipid concentration. Also the error bars for each reduced scattering coefficient are reported. Finally, we calculated the intrinsic reduced scattering coefficient *ε*′_{sil} of the Intralipid-20%, by performing the best linear fit of these data, as stated in Eq. (9):

The relative error in *ε*′_{sil} is 1.5%. In comparison, the uncertainty due to preparation of the
suspension (see Sec. 4.1) can be neglected. Then, this is the accuracy that would characterize a diffusive sample prepared with this calibrated batch of Intralipid-20%. On the contrary, with a single multi-distance measurement the reduced scattering coefficient can be measured with a relative error of 3%, as stated above.

#### 4.4 Comparison with CW results

Because for both the measurements in the CW and in the time domain we have used the same Intralipid-20% and Indian ink batches and we have considered the same wavelength (750 nm), we can compare the retrieved intrinsic absorption and reduced scattering coefficients, *ε*
_{aink} and *ε*′_{sil}, respectively, obtained by applying the calibration methods in the two experimental approaches. In table 1 the results obtained with the two independent procedures are summarized. The best estimations of the uncertainty affecting the results are also reported. As one can see from the table, the results obtained with the two calibration methods are in agreement within the error bars. A discussion on how these values compare with similar results reported in literature is reported in [4].

## 5. Conclusions

In this paper, we have presented a method to calibrate accurately the absorption and scattering properties of an absorber and scatterer compound that can be suspended in a liquid solution. The method exploits a time-resolved experimental set-up and simple analysis procedures based on linear regressions. An extended discussion on all the possible sources of uncertainty affecting the calibration procedure has been performed. What is found is that the distortion introduced by the instrumental response function is one of the main sources of uncertainty affecting the estimation of both the absorption and scattering properties.

By applying this method with a time-resolved spectroscopy system characterized by a sharp IRF (FWHM of about 50 ps), we determined the intrinsic absorption coefficient of a sample of Indian ink and the intrinsic reduced scattering coefficient of Intralipid-20% with an error that we evaluated to be less than 0.5% and 1.5%, respectively. On the other hand, in the companion paper recently published [4], the same Indian ink and Intralipid-20% batches have been calibrated at the same wavelength, with a procedure based on CW instrumentation that allows standard errors smaller than 2%. The results obtained by the two calibration methods agree within the error bars and this can be considered as a cross validation of the two proposed methods, due to the complete independence of the two methods.

Intralipid and Indian ink are widely used for the preparation of biological tissue phantoms. The usefulness to have calibrated phantoms to work with is quite clear and deeply discussed in [4]. However, one could argue that both the calibration methods relies on liquid phantoms, that are difficult to use in some cases and not so handy as the solid ones, and, for the time-domain procedure, on an instrumentation that is the state-of-the-art for what concerns the time response. Nevertheless, from one side, the procedure proposed in time domain to determine the reduced scattering coefficient can be applied also on solid phantoms, and, from the other side, with these calibrated liquid phantoms in the laboratory, one can characterize the accuracy performance of his instrument once and for all, in all the situations one meets in applications. Both the data acquisition and the analysis procedure, that will be non linear in the majority of the cases, due to the convolution with the instrument response function, must be taken into account. With a calibrated instrument, for example, one can assign to a series of solid phantoms the correct optical properties with their uncertainties. Following this philosophy, work is in progress in our group with the final goal to accurately calibrate the solid phantoms of the MEDPHOT protocol [14].

## Acknowledgments

The work was partially supported by MIUR under the project PRIN2005 (prot. 2005025333).

## References and links

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

**2. **H. Xu and M. S. Patterson, “Determination of the optical properties of tissue-simulating phantoms from interstitial frequency domain measurements of relative fluence and phase difference,” Opt. Express **14**, 6485 (2006). [CrossRef] [PubMed]

**3. **C. Chen, J. Q. Lu, H. Ding, K. M. Jacobs, Y. Du, and X. -H. Hu, “A primary method for determination of optical parameters of turbid samples and application to intralipid between 550 and 1630 nm,” Opt. Express **14**, 7420 (2006). [CrossRef] [PubMed]

**4. **F. Martelli and G. Zaccanti, “Calibration of scattering and absorption properties of a liquid diffusive medium at NIR wavelengths. CW method,” Opt. Express **15**, 486–500 (2007). [CrossRef] [PubMed]

**5. **Y. Hasegawa*et al*., “Monte Carlo simulation of light transmission through living tissues,” Appl. Opt. **30**, 4515 (1991). [CrossRef] [PubMed]

**6. **Y. Nomura, O. Hazeki, and M. Tamura, “Relationship between time-resolved and non-time-resolved Beer-Lambert law in turbid media,” Phys. Med. Biol. **42**, 1009 (1997). [CrossRef] [PubMed]

**7. **F. Martelli, A. Sassaroli, Y. Yamada, and G. Zaccanti, “Method for measuring the diffusion coefficient of homogeneous and layered media,” Opt. Lett. **25**, 1508 (2000). [CrossRef]

**8. **R. K. Wang and Y. A. Wickramasinghe, “Fast algorithm to determine optical properties of a turbid medium from time-resolved measurements,” Appl. Opt. **37**, 7342 (1998). [CrossRef]

**9. **D. Contini, F. Martelli, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. **36**, 4587 (1997). [CrossRef] [PubMed]

**10. **D. Contini, A. Torricelli, A. Pifferi, L. Spinelli, F. Paglia, and R. Cubeddu, “Multi-channel time-resolved system for functional near infrared spectroscopy,” Opt. Express **14**, 5418 (2006). [CrossRef] [PubMed]

**11. **R. Cubeddu, A. Pifferi, P. Taroni, A. Torricelli, and G. Valentini, “Noninvasive absorption and scattering spectroscopy of bulk diffusive media: An application to the optical characterization of human breast,” Appl. Phys. Lett. **74**, 874 (1999). [CrossRef]

**12. **W. Becker, *Advanced Time-Correlated Single Photon Counting Techniques* (Spinger-Verlag, Berlin Heidelberg, 2005). [CrossRef]

**13. **J. R. Taylor, *An Introduction to Error Analysis, The Study of Uncertainties in Physical Measurements* (University Science Books, 1982).

**14. **A. Pifferi, A. Torricelli, A. Bassi, P. Taroni, R. Cubeddu, H. Wabnitz, D. Grosenick, M. Möller, R. Macdonald, J. Swartling, T. Svensson, S. Andersson-Engels, R. L. P. van Veen, H. J. C. M. Sterenborg, J. M. Tualle, H. L. Nghiem, E. Tinet, S. Avrillier, M. Whelan, and H. Stamm, “Performance assessment of photon migration instruments: the MEDPHOT protocol,” Appl. Opt. **44**, 2104 (2005). [CrossRef] [PubMed]