## 1.

## Introduction

Recently, short-pulse lasers have been widely used in noninvasive optical tomography.^{1, 2, 3, 4, 5, 6} As compared with continuous wave imaging, more information can be achieved in short-pulse imaging by studying the temporal distribution of reflectance/transmittance.^{7, 8, 9, 10, 11} The temporal distribution of a reflected/transmitted signal is broadened due to multiple scattering of photons into tissues.^{2, 12} The level of broadening depends on optical properties of the tissue, i.e., absorption and scattering; consequently, variations of those parameters can be investigated by studying the temporal distribution of the scattered light. Thus, the temporal distribution of the reflectance/transmittance can be exploited to detect normal or malignant biological tissue, because the optical properties of normal tissues vary during malignant progress.^{13}

Photon transport through biological tissues can be described by a diffuse equation. Sometimes, this equation can be solved analytically.^{14} But in most cases, numerical methods such as the Monte Carlo (MC) method, the finite-element method (FEM), and the finite-difference time-domain (FDTD) method are most often used to solve that equation.^{15, 16, 17, 18, 19, 20, 21, 22} In those methods, the whole sample should be discritized and the calculation is time consuming.^{15}

The boundary integral method (BIM) can also be used to solve a diffuse equation and to simulate a reflected pulse on the surface of biological tissues. Since this method requires surface tessellation, the computation time is reduced and the accuracy of the results increases as compared to other numerical methods.^{23, 24, 25, 26} The effects of optical properties of biological tissues on the intensity of a reflected pulse were investigated and reported in Ref. 26, but the behavior of a pulse penetrating through biological tissues has not yet been studied by this method.

In this paper, pulse penetration into biological tissue is studied using the BIM. First, the appropriate Green's function is obtained to convert the diffuse equation to integral form using Green's second theorem. Then, the surface integral is discretized using the boundary element method (BEM). In this method, the boundary of the sample is first discretized to elements. Then, the observation point is located on the surface of tissue and an equation containing fluence at that point is achieved. Locating the observation point on different nodes, a system of equations is obtained that gives the fluence at those points.^{23, 24, 25} By using this technique, the photon intensity inside the sample and the intensity of the diffusely penetrated pulse inside tissue are calculated. To investigate the accuracy and precision of results, they are compared with those obtained analytically and by other numerical methods. Furthermore, the effects of the absorption and scattering coefficients and the anisotropic factor on the diffusely penetrating pulse are also studied. Finally, the penetration of the short pulse with a arbitrary time shape is studied.

## 2.

## Review of Theory

Short-pulse laser propagation into biological tissue Ω with boundary
[TeX:]
$\Gamma _s$
${\Gamma}_{s}$
can be studied by a diffuse equation and Robin's boundary condition, which are, respectively, given by^{16}

## Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{\partial }{{c\partial t}}\varphi \left({{\bm r}\!,t} \right) - D\nabla ^2 \varphi \left({{\bm r}\!,t} \right) + a\varphi \left({{\bm r}\!,t} \right) = S\left({{\bm r}\!,t} \right)\quad{\bm r} \in \Omega, \end{equation}\end{document} $$\frac{\partial}{c\partial t}\varphi \left(r,t\right)-D{\nabla}^{2}\varphi \left(r,t\right)+a\varphi \left(r,t\right)=S\left(r,t\right)\phantom{\rule{1em}{0ex}}r\in \Omega ,$$## Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \hspace*{25pt}\varphi \left({{\bm r}\!,t} \right) - 2C_R D\frac{{\partial \varphi \left({{\bm r}\!,t} \right)}}{{\partial z}} = 0\quad {\bm r} \in \Gamma _s, \end{equation}\end{document} $$\phantom{\rule{25.0pt}{0ex}}\varphi \left(r,t\right)-2{C}_{R}D\frac{\partial \varphi \left(r,t\right)}{\partial z}=0\phantom{\rule{1em}{0ex}}r\in {\Gamma}_{s},$$where
[TeX:]
$\varphi \left({{\bm r}\!,t} \right)$
$\varphi \left(r,t\right)$
and
[TeX:]
$S\left({{\bm r}\!,t} \right)$
$S\left(r,t\right)$
are, respectively, the fluence and the isotropic source term at position
[TeX:]
${\bm r}$
$r$
and at moment *t*. The velocity of light is shown by *c*. The parameter
[TeX:]
$D = {1 / {3({a + \sigma^ \prime })}}$
$D=1/3\left(a+{\sigma}^{\prime}\right)$
is the diffusion coefficient, where *a* and
[TeX:]
$\sigma^ \prime = \sigma ({1 - g})$
${\sigma}^{\prime}=\sigma \left(1-g\right)$
are the absorption and reduced scattering coefficients, respectively; and σ and *g* are also the scattering coefficient and the anisotropic factor, respectively. In Eq. 2,
[TeX:]
$C_R = {{({1 + R})} / {({1 - R})}}$
${C}_{R}=\left(1+R\right)/\left(1-R\right)$
, where *R* is the Fresnel reflection coefficient.

The boundary integral method is based on using of Green's function.^{23} The Green's function of Eq. 1 in domain Ω is the solution of

## Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} &&\frac{\partial }{{c\partial t}}g({{\bm r}\!,{\bm r'};t,t'}) - D\nabla ^2 g({{\bm r}\!,{\bm r'};t,t'}) + ag({{\bm r}\!,{\bm r'};t,t'})\nonumber \\ && \quad= - \delta ({{\bm r} - {\bm r}'})\delta ({t - t'}). \end{eqnarray}\end{document} $$\begin{array}{ccc}& & \frac{\partial}{c\partial t}g\left(r,{r}^{\prime};t,{t}^{\prime}\right)-D{\nabla}^{2}g\left(r,{r}^{\prime};t,{t}^{\prime}\right)+ag\left(r,{r}^{\prime};t,{t}^{\prime}\right)\hfill \\ & & \phantom{\rule{1em}{0ex}}=-\delta \left(r-{r}^{\prime}\right)\delta \left(t-{t}^{\prime}\right).\hfill \end{array}$$*t*in Eq. 3 and rearranging the resulting equation, gives

## Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} &&\hspace*{-6pt} \nabla ^2 G({{\bm r}\!,{\bm r'};s,t'})-\kappa ^2 G({{\bm r}\!,{\bm r'};s,t'}) = \frac{{\delta ({{\bm r} - {\bm r'}})}}{D}\,\exp\, ({ - st'}),\phantom{00} \end{eqnarray}\end{document} $$\begin{array}{ccc}& & {\nabla}^{2}G\left(r,{r}^{\prime};s,{t}^{\prime}\right)-{\kappa}^{2}G\left(r,{r}^{\prime};s,{t}^{\prime}\right)=\frac{\delta \left(r-{r}^{\prime}\right)}{D}\phantom{\rule{0.16em}{0ex}}\mathrm{exp}\phantom{\rule{0.16em}{0ex}}\left(-s{t}^{\prime}\right),\phantom{00}\hfill \end{array}$$## Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \tilde G({k,{\bm r}';s,t'}) = \frac{{\exp\, ({ - ik.{\bm r}'})}}{{D({k^2 + \kappa ^2 })}}\exp\, ({ - st'}), \end{equation}\end{document} $$\stackrel{\u0303}{G}\left(k,{r}^{\prime};s,{t}^{\prime}\right)=\frac{\mathrm{exp}\phantom{\rule{0.16em}{0ex}}\left(-ik.{r}^{\prime}\right)}{D\left({k}^{2}+{\kappa}^{2}\right)}\mathrm{exp}\phantom{\rule{0.16em}{0ex}}\left(-s{t}^{\prime}\right),$$*G*.

The inverse Fourier and Laplace transforms of Green's function stated in Eq. 5 gives the Green's function as

## Eq. 6

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \hspace*{-16pt}g({{\bm r},{\bm r}';t,t'}) &=& H({t - t'})\frac{c}{{{\{ {[ {4\pi {\rm }Dc{\rm }({t - t'})} ]} ^3\}}^{1/2}}} \nonumber\\ && \times\, \exp[ { - ca({t - t'})} ]\exp \left[ { - \frac{{\left| {{\bm r} - {\bm r}'} \right|^2 }}{{4Dc{\rm }\left({t - t'} \right)}}} \right], \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill g\left(r,{r}^{\prime};t,{t}^{\prime}\right)& =& H\left(t-{t}^{\prime}\right)\frac{c}{{\left\{{\left[4\pi Dc\left(t-{t}^{\prime}\right)\right]}^{3}\right\}}^{1/2}}\hfill \\ & & \times \phantom{\rule{0.16em}{0ex}}\mathrm{exp}\left[-ca\left(t-{t}^{\prime}\right)\right]\mathrm{exp}\left[-\frac{{\left|r-{r}^{\prime}\right|}^{2}}{4Dc\left(t-{t}^{\prime}\right)}\right],\hfill \end{array}$$For a short-point pulse where
[TeX:]
$S\left({{\bm r}\!,t} \right) = \delta \left({\bm r} \right)\delta \left(t \right)$
$S\left(r,t\right)=\delta \left(r\right)\delta \left(t\right)$
, the fluence at any arbitrary point inside an infinite medium for *t* > 0 is equal to^{26}

## Eq. 7

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\rm }\varphi ({{\bm r}\!,t}) = \frac{c}{{ {[ {({4\pi {\rm }Dct})} ^3]^{1/2} } }}\exp ({ - cat})\exp \left({ - \frac{{| {\bm r} |^2 }}{{4Dct}}} \right). \end{equation}\end{document} $$\varphi \left(r,t\right)=\frac{c}{{\left[{\left(4\pi Dct\right)}^{3}\right]}^{1/2}}\mathrm{exp}\left(-cat\right)\mathrm{exp}\left(-\frac{{\left|r\right|}^{2}}{4Dct}\right).$$Assuming the sample as a slab and by using Eq. 7, the transmited intensity from it is given by^{27}

## Eq. 8

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} T\left({d,t} \right) &=& \left({4\pi Dct} \right)^{{{ - 1}/2}} t^{{{ - 3}/2}} \,{\rm exp}\,\left({ - act} \right) \nonumber \\ && \times\, \left\{ {\left({d - z_0 } \right)\,\exp \,\left[ { - \displaystyle\frac{{\left({d - z_0 } \right)^2 }}{{4Dct}}} \right]} \right. \nonumber \\ && -\, \left({d + z_0 } \right)\,\exp \,\left[ { - \displaystyle\frac{{\left({d + z_0 } \right)^2 }}{{4Dct}}} \right] \nonumber\\ && +\, \left({3d - z_0 } \right)\,\exp \,\left[ { - \displaystyle\frac{{\left({3d - z_0 } \right)^2 }}{{4Dct}}} \right] \nonumber \\ && -\, \left.\left({3d + z_0 } \right)\,\exp \,\left[ { - \displaystyle\frac{{\left({3d + z_0 } \right)^2 }}{{4Dct}}} \right] \right\}, \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill T\left(d,t\right)& =& {\left(4\pi Dct\right)}^{-1/2}{t}^{-3/2}\phantom{\rule{0.16em}{0ex}}\mathrm{exp}\phantom{\rule{0.16em}{0ex}}\left(-act\right)\hfill \\ & & \times \phantom{\rule{0.16em}{0ex}}\left\{\left(d-{z}_{0}\right)\phantom{\rule{0.16em}{0ex}}\mathrm{exp}\phantom{\rule{0.16em}{0ex}}\left[{\displaystyle -\frac{{\left(d-{z}_{0}\right)}^{2}}{4Dct}}\right]\right.\hfill \\ & & -\phantom{\rule{0.16em}{0ex}}\left(d+{z}_{0}\right)\phantom{\rule{0.16em}{0ex}}\mathrm{exp}\phantom{\rule{0.16em}{0ex}}\left[{\displaystyle -\frac{{\left(d+{z}_{0}\right)}^{2}}{4Dct}}\right]\hfill \\ & & +\phantom{\rule{0.16em}{0ex}}\left(3d-{z}_{0}\right)\phantom{\rule{0.16em}{0ex}}\mathrm{exp}\phantom{\rule{0.16em}{0ex}}\left[{\displaystyle -\frac{{\left(3d-{z}_{0}\right)}^{2}}{4Dct}}\right]\hfill \\ & & -\phantom{\rule{0.16em}{0ex}}\left.\left(3d+{z}_{0}\right)\phantom{\rule{0.16em}{0ex}}\mathrm{exp}\phantom{\rule{0.16em}{0ex}}\left[{\displaystyle -\frac{{\left(3d+{z}_{0}\right)}^{2}}{4Dct}}\right]\right\},\hfill \end{array}$$*d*is the thickness of the slab, and [TeX:] $z_0 = 1/\sigma ^\prime $ ${z}_{0}=1/{\sigma}^{\prime}$ . For an extended source, the fluence is obtained by using Green's second theorem:

## Eq. 9

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \varphi ({{\bm r}\!,t}) &=& \displaystyle\int\nolimits^t_\Omega\!\! {\int\nolimits_0 {S({{\bm r}',t'})} } g ({{\bm r}\!,t;{\bm r}',t'})\,{\rm d}{\bm r}'dt' \nonumber\\ && - \displaystyle\int\nolimits^t_{\Gamma _s }\! \int\nolimits_0 \Bigg\{\left[ \frac{c}{C_R D}\varphi ({\bm r}_s, t')g({\bm r}\!,t;{\bm r}_s, t') \right]\nonumber \\ && - \left[ \varphi ({\bm r}_s, t')\frac{\partial }{\partial n}g ({\bm r}\!,t;{\bm r}_s, t') \right]\Bigg\} {\rm d}{\bm r}_s\, dt', \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill \varphi \left(r,t\right)& =& {\displaystyle {\int}_{\Omega}^{t}{\int}_{0}S\left({r}^{\prime},{t}^{\prime}\right)g\left(r,t;{r}^{\prime},{t}^{\prime}\right)\phantom{\rule{0.16em}{0ex}}\mathrm{d}{r}^{\prime}d{t}^{\prime}}\hfill \\ & & {\displaystyle -{\int}_{{\Gamma}_{s}}^{t}{\int}_{0}\{\left[\frac{c}{{C}_{R}D}\varphi ({r}_{s},{t}^{\prime})g(r,t;{r}_{s},{t}^{\prime})\right]}\hfill \\ & & -\left[\varphi ({r}_{s},{t}^{\prime})\frac{\partial}{\partial n}g(r,t;{r}_{s},{t}^{\prime})\right]\}\mathrm{d}{r}_{s}\phantom{\rule{0.16em}{0ex}}d{t}^{\prime},\hfill \end{array}$$The intensity of a diffused short pulse can be calculated from Eq. 9, which is the integral form of diffuse equation. The BEM can be used to solve Eq. 9. In this method, the boundary of the sample is first discretized to elements. Then, observation point [TeX:] ${\bm r}_s$ ${r}_{s}$ is located on the surface of tissue and an equation containing fluence at that point is achieved. Locating observation point on different nodes, a system of equations is obtained which gives the fluence at those points. Thus, the boundary [TeX:] $\Gamma _s^{}$ ${\Gamma}_{s}^{}$ is discretized to square elements and the fluence φ is approximated as

## Eq. 10

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \varphi \left({{\bm r}\!,t} \right) = \sum\limits_{k = 1}^n {\sum\limits_{l = 1}^m {N_k \left({\bm r} \right)\varphi _{kl} \left({{\bm r}\!,t} \right)} }, \end{equation}\end{document} $$\varphi \left(r,t\right)=\sum _{k=1}^{n}\sum _{l=1}^{m}{N}_{k}\left(r\right){\varphi}_{kl}\left(r,t\right),$$*k*and

*l*refer to node

*k*and time step

*l*, and [TeX:] $N_i ({{\bm r}_j }) = \delta _{ij}$ ${N}_{i}\left({r}_{j}\right)={\delta}_{ij}$ is the Kronecker symbol. If [TeX:] ${\bm r}_j$ ${r}_{j}$ spans all the nodes on the surface of the boundary, Eqs. 1, 2 will give, respectively, the following set of algebraic equations:

## Eq. 11

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{l} H{\bm u} + \Gamma {\bm v} = \bar {\bm S}, \\[2pt] {\bm v} = - R{\bm u} + P, \\ \end{array} \end{equation}\end{document} $$\begin{array}{c}Hu+\Gamma v=\overline{S},\hfill \\ v=-Ru+P,\hfill \end{array}$$*q*, the prescribed boundary flux

*p*, and the volume source

*S*, respectively, which are given by

## Eq. 12

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} {\bm u}_{n \times 1} &=& \left[ \begin{array}{l} \phi _1 \\ {\cdots} \\ {\cdots} \\ \phi _n \\ \end{array} \right], \quad {\bm v}_{n \times 1} = \left[ \begin{array}{l} q_1 \\ {\cdots} \\ {\cdots} \\ q_n \\ \end{array} \right] , \nonumber\\ {\bm P}_{n \times 1} & =&\left[ \begin{array}{l} p_1 \\ {\cdots} \\ {\cdots} \\ p_n \\ \end{array} \right] , \quad\bar {\bm S}_{n \times 1} = \left[ \begin{array}{l} s_1 \\ {\cdots} \\ {\cdots} \\ s_n \\ \end{array} \right]{\rm, } \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill {u}_{n\times 1}& =& \left[\begin{array}{c}{\phi}_{1}\hfill \\ \cdots \hfill \\ \cdots \hfill \\ {\phi}_{n}\hfill \end{array}\right],\phantom{\rule{1em}{0ex}}{v}_{n\times 1}=\left[\begin{array}{c}{q}_{1}\hfill \\ \cdots \hfill \\ \cdots \hfill \\ {q}_{n}\hfill \end{array}\right],\hfill \\ \hfill {P}_{n\times 1}& =& \left[\begin{array}{c}{p}_{1}\hfill \\ \cdots \hfill \\ \cdots \hfill \\ {p}_{n}\hfill \end{array}\right],\phantom{\rule{1em}{0ex}}{\overline{S}}_{n\times 1}=\left[\begin{array}{c}{s}_{1}\hfill \\ \cdots \hfill \\ \cdots \hfill \\ {s}_{n}\hfill \end{array}\right],\hfill \end{array}$$## Eq. 13

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} s_j = \int\nolimits_{\Gamma _s }\!\! {\int\limits_{} {g({\rho; t,t'})S({\bm r},{t})\ {\rm d}rdt',} } {\rm } \end{equation}\end{document} $${s}_{j}={\int}_{{\Gamma}_{s}}\underset{}{\int}g\left(\rho ;t,{t}^{\prime}\right)S(r,t)\phantom{\rule{4pt}{0ex}}\mathrm{d}rd{t}^{\prime},$$## Eq. 14

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} h_{i,j,l} &=& \delta _{ij} I + \displaystyle\int\nolimits^t_{\Gamma _s }\! {\int\nolimits_0 {\frac{{\partial g\left({\rho; \tau } \right)}}{{\partial n}}N_{kl} \left({{\bm r}\!,t} \right)\,{\rm d}{\bm r}\,{\rm d}t'} }, \nonumber \\[-9pt] \\[-9pt]\nonumber \xi _{i,j,l} &=& -\displaystyle \int\nolimits^t_{\Gamma _s }\! {\int\nolimits_0 {g\left({\rho; \tau } \right)N_{kl} \left({{\bm r}\!,t} \right){\rm d}{\bm r}\,{\rm d}t'}, } \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill {h}_{i,j,l}& =& {\displaystyle {\delta}_{ij}I+{\int}_{{\Gamma}_{s}}^{t}{\int}_{0}\frac{\partial g\left(\rho ;\tau \right)}{\partial n}{N}_{kl}\left(r,t\right)\phantom{\rule{0.16em}{0ex}}\mathrm{d}r\phantom{\rule{0.16em}{0ex}}\mathrm{d}{t}^{\prime},}\hfill \\ \\ \hfill {\xi}_{i,j,l}& =& {\displaystyle -{\int}_{{\Gamma}_{s}}^{t}{\int}_{0}g\left(\rho ;\tau \right){N}_{kl}\left(r,t\right)\mathrm{d}r\phantom{\rule{0.16em}{0ex}}\mathrm{d}{t}^{\prime},}\hfill \end{array}$$## Eq. 15

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} %\begin{array}{l} \varphi ^* ({{\bm r}\!,t}) &=&\displaystyle \int\nolimits^t_\Omega \!\!{\int\nolimits_0 {S ({{\bm r}',t'})} } g ({{\bm r}\!,t;{\bm r}',t'})\,{\rm d}{\bm r}'\,{\rm d}t' \nonumber \\ && - \displaystyle\int\nolimits^t_{\Gamma _s } \!\int\nolimits_0 \Bigg\{\left[ \frac{c}{{C_R D}}\bar \varphi ({{\bm r}_s, t'})g ({{\bm r}\!,t;{\bm r}_s, t'}) \right] \nonumber \\ &&- \left[ \bar \varphi ({{\bm r}_s, t'})\frac{\partial }{{\partial n}}g ({{\bm r}\!,t;{\bm r}_s, t'}) \right]\Bigg\} \,{\rm d}{\bm r}_s \,{\rm d}t', % \end{array} \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill {\varphi}^{*}\left(r,t\right)& =& {\displaystyle {\int}_{\Omega}^{t}{\int}_{0}S\left({r}^{\prime},{t}^{\prime}\right)g\left(r,t;{r}^{\prime},{t}^{\prime}\right)\phantom{\rule{0.16em}{0ex}}\mathrm{d}{r}^{\prime}\phantom{\rule{0.16em}{0ex}}\mathrm{d}{t}^{\prime}}\hfill \\ & & {\displaystyle -{\int}_{{\Gamma}_{s}}^{t}{\int}_{0}\{\left[\frac{c}{{C}_{R}D}\overline{\varphi}\left({r}_{s},{t}^{\prime}\right)g\left(r,t;{r}_{s},{t}^{\prime}\right)\right]}\hfill \\ & & -\left[\overline{\varphi}\left({r}_{s},{t}^{\prime}\right)\frac{\partial}{\partial n}g\left(r,t;{r}_{s},{t}^{\prime}\right)\right]\}\phantom{\rule{0.16em}{0ex}}\mathrm{d}{r}_{s}\phantom{\rule{0.16em}{0ex}}\mathrm{d}{t}^{\prime},\hfill \end{array}$$*t*. Note that to calculate of the fluence at each internal points, we only need to use the solution of the system of equations that were obtained on the boundary, the values of [TeX:] $\bar \varphi _{}$ ${\overline{\varphi}}_{}$ , which we now know, so Eq. 15 can be used as many times as desired, with only the surface integration required to be calculated for every new internal point.

## 3.

## Numerical Results

The diffuse equation is used to study propagation of a short-pulse laser in biological tissues. We use the BIM to solve diffuse equation and to study behavior of light inside the biological tissues.

First, to confirm precision of this method, the results obtained by the BIM are compared with those obtained analytically using Eq. 8 and the MC and FDTD methods. An interesting parameter *in vivo* time domain tomography is 〈*ct*〉, i.e., the mean distance traversed by photons before exiting the tissues. Using Eq. 8, 〈*ct*〉 can be analytically found to be equal to 77.6 mm for a 10-mm slab tissue with
[TeX:]
$a = 0.00434\ {\rm mm}^{{\rm - 1}}$
$a=0.00434\phantom{\rule{4pt}{0ex}}{\mathrm{mm}}^{-1}$
,
[TeX:]
$\sigma = 6\ {\rm mm}^{{\rm - 1}}$
$\sigma =6\phantom{\rule{4pt}{0ex}}{\mathrm{mm}}^{-1}$
, and *g* = 0.72. Patterson calculated this quantity as 80.6 mm, an error of 3.8% by the MC method.^{27} The same quantity is calculated as method 81.6 mm, an error of only 5.1%, and 89.3 mm, an error of 15%, using the BIM and FDTD method, respectively. However, the computational time of the FDTD method is more 4 time longer than the BIM, whereas the mesh resolution in the FDTD method and the BIM are 2.00 pS × 0.45 mm and 12.19 pS × 1.36 mm, respectively.

Note that Eqs. 7, 8 can only be used for point sources. However, it is not appropriate to study propagation of a pulse originating from sources that are temporally extended. We assume a semi-infinite sample with
[TeX:]
$a = 0.02\ {\rm mm}^{{\rm - 1}}$
$a=0.02\phantom{\rule{4pt}{0ex}}{\mathrm{mm}}^{-1}$
,
[TeX:]
$\sigma = 10\ {\rm mm}^{{\rm - 1}}$
$\sigma =10\phantom{\rule{4pt}{0ex}}{\mathrm{mm}}^{-1}$
, and *g* = 0.9, which is illuminated by a Gaussian pulse that is presented by

## Eq. 16

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} I\left({{\bm r}\!,t} \right) = I_0 \,{\rm exp}\,\left[ { - \left({\frac{{t - t_0 }}{\tau }} \right)^2 } \right]{\rm }\delta \left({\bm r} \right), \end{equation}\end{document} $$I\left(r,t\right)={I}_{0}\phantom{\rule{0.16em}{0ex}}\mathrm{exp}\phantom{\rule{0.16em}{0ex}}\left[-{\left(\frac{t-{t}_{0}}{\tau}\right)}^{2}\right]\delta \left(r\right),$$Next, we study the effect of different parameters on the temporal evolution of a diffusely penetrating pulse into a phantom like breast tissue with optical properties similar to that represented in Ref. 28. Temporal broadening of the penetrating pulse is due to scattering and absorption of diffused photons in tissue. Therefore, it is important to study such an effect in more detail. To study the effect of multiple scattering on the diffused pulse, we assume a semi-infinite sample with [TeX:] $a = 0.02\ {\rm mm}^{{\rm - 1}}$ $a=0.02\phantom{\rule{4pt}{0ex}}{\mathrm{mm}}^{-1}$ to be illuminated by a Gaussian pulse with duration of τ = 10.0 ps. Figure 2 shows that the peak value of the pulse decreases for larger scattering coefficients. This is because by increasing the scattering coefficient, more photons are scattered from the illumination direction, and consequently, the intensity of penetrating pulse decreases. The results are in agreement with those reported in Ref. 29.

Figure 3 illustrates the temporal evolution of a penetrating pulse for two different values of the anisotropic factors *g*. One can see that the peak value of the pulse increases for larger anisotropic factor values. This is because for a larger anisotropic factor, the majority of launched photons are scattered along the illumination direction, and therefore, its intensity increases.

Furthermore, to study effect of the absorption coefficient on the penetrating pulse, the temporal distribution of penetrating pulse is calculated and the results for three different absorption coefficients values are calculated and depicted in Fig. 4. The results show that by increasing the absorption coefficient, more diffused photons are absorbed in tissue and the density of the remaining photons decreases.

As mentioned earlier, using this technique, the temporal evolution of the source can take arbitrary shape, and we also consider a source emitting pulses with the following temporal evolution:

## Eq. 17

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} I ({{\bm r}\!,t}) = I_0 \Bigg\{ {\rm exp}\!\left[ {-} \!\left(\frac{t - t_0 }{\tau } \right)^2 \right] + {\rm exp}\!\left[ - \left (\frac{t - 2t_0 }{2\tau } \right) \!\right]\!\Bigg\} \delta ({\bm r}),\nonumber\\ \end{eqnarray}\end{document} $$\begin{array}{c}\hfill I\left(r,t\right)={I}_{0}\left\{\mathrm{exp}\left[-{\left(\frac{t-{t}_{0}}{\tau}\right)}^{2}\right]+\mathrm{exp}\left[-\left(\frac{t-2{t}_{0}}{2\tau}\right)\right]\right\}\delta \left(r\right),\end{array}$$## 4.

## Conclusion

Recently, the diffuse equation was solved by the BIM (Refs. 24 to 26), and note that another paper^{30} solved the transient diffuse equation for turbine component by means of the BEM.^{30} In Ref. 26, the BIM solution of transient diffuse equation for biological tissues was presented and the accuracy of this method was inspected and the showed resulting emphasizing that the BIM method is faster than the MC and FDTD methods and more precise than the FDTD method. In Ref. 26, only the reflected pulse on the surface of the tissue was calculated.

In this paper, the propagation of short pulses inside biological tissue such as breast was studied. The BIM was used to study the diffusion of short pulses penetrating inside biological samples and to calculate the photon density at any arbitrary point. The effects of the optical properties of tissues on diffusely penetrating pulses inside it were studied and it was observed those parameters can alter the intensity and temporal distribution of the penetrating pulse. Since the optical properties of malignant tissues gradually changes, this offers an effective technique to trace the progress of cancer.

Using this method, the propagation of pulse emitted from extended sources with an arbitrary temporal evolution can be studied, where in this paper, sources with Gaussian shape were studied.

This method is suitable to study the interaction of a short-pulse laser with biological tissue, especially for regions where shock waves are produced.

## References

**,” Appl. Opt., 32 554 –558 (1993). https://doi.org/10.1364/AO.32.000554 Google Scholar**

*Ultrafast laser-pulse transmission and imaging through biological tissues***,” Opt. Lett., 27 336 –338 (2002). https://doi.org/10.1364/OL.27.000336 Google Scholar**

*Analysis of short-pulse laser photon transport through tissues for optical tomography***,” Opt. Express., 16 10434 –10448 (2008). https://doi.org/10.1364/OE.16.010440 Google Scholar**

*``Improved accucuracy in timedomain diffuse reflectance spectroscopy***,” Phys. Med. Biol., 52 2827 –2843 (2007). https://doi.org/10.1088/0031-9155/52/10/013 Google Scholar**

*Solution of time-dependent diffuse equation for a three-layer medium: application to study photon migration through a simplified adult head model***,” Appl. Opt., 39 4411 –4417 (2000). https://doi.org/10.1364/AO.39.004411 Google Scholar**

*Equitant isotropic scattering formulation for transient short-pulse radiative transfer in anisotropic scattering planar media***,” J. Quantum Spectrosc. Radiat. Transf., 93 337 –348 (2005). https://doi.org/10.1016/j.jqsrt.2004.08.033 Google Scholar**

*Time-resolved optical tomography using short pulse laser for tumor detection***,” Opt. Rev., 10 609 –610 (2003). https://doi.org/10.1007/s10043-003-0609-3 Google Scholar**

*Time-resolved optical tomography using short pulse laser for tumor detection***,” J. Phys. D Appl. Phys., 38 1714 –1721 (2003). https://doi.org/10.1088/0022-3727/36/14/310 Google Scholar**

*Short pulse laser propagation through tissues for biomedical imagin***,” Appl. Opt., 38 188 –196 (1999). https://doi.org/10.1364/AO.38.000188 Google Scholar**

*Development and comparison of models for light pulse transport through scattering absorbing media***,” Appl. Opt., 36 7270 –7276 (1997). https://doi.org/10.1364/AO.36.007270 Google Scholar**

*Imaging very-low-contrast objects in breastlike scattering media with a time resolved method***,” Opt. Lett., 15 320 –322 (1990). https://doi.org/10.1364/OL.15.000320 Google Scholar**

*Time resolved coherent and incoherent component of forward light scattering in random media***,” Appl. Opt., 45 6270 –6282 (2006). Google Scholar**

*Time-resolved optical tomography using short-pulse laser for tumor detection***,” Waves Rand. Med., 9 551 –560 (1999). https://doi.org/10.1088/0959-7174/9/4/307 Google Scholar**

*Optical pulse propagation through a slab of random medium***,” Opt. Express., 10 159 –170 (2002). Google Scholar**

*Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head***,” Appl. Opt., 30 4515 –4520 (1991). https://doi.org/10.1364/AO.30.004515 Google Scholar**

*Monte Carlo simulation of light transmission through living tissues***,” J. Thermophys. Heat Transf., 14 504 –511 (2000). https://doi.org/10.2514/2.6573 Google Scholar**

*``Multidimensional Monte Carlo simulation of short-pulse laser transport in scattering media***,” Phys. Med. Biol., 43 1285 –1302 (1998). Google Scholar**

*Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues***,” J. Quant. Spectrosc. RA, 73 159 –168 (2002). https://doi.org/10.1016/S0022-4073(01)00203-5 Google Scholar**

*Monte Carlo simulation and experiments of pulsed radiative transfer***,” IEEE. Trans. Med. Imaging, 21 181 –184 (2002). https://doi.org/10.1109/42.993136 Google Scholar**

*Finite difference time domain (FDTD) analysis of optical pulse responses in biological tissues for spectroscopic diffused optical tomography***,” Phys. Med. Biol, 44 2747 –2763 (1999). https://doi.org/10.1088/0031-9155/44/11/305 Google Scholar**

*Study on the propagation of ultra-short pulse light in cylindrical optical phantoms***,” Med. Phys., 34 4545 –4557 (2007). https://doi.org/10.1118/1.2795832 Google Scholar**

*Estimation***,” Opt. Laser. Eng., 47 965 –970 (2009). https://doi.org/10.1016/j.optlaseng.2009.04.006 Google Scholar**

*Study of light propagation in Asian and Caucasian skins by means of the boundary element method***,” Google Scholar**

*Study of short pulse laser in biological tissue by means of boundary element method***,” Appl. Opt., 28 2331 –2336 (1989). https://doi.org/10.1364/AO.28.002331 Google Scholar**

*Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties***,” Appl. Opt., 40 176 –184 (2001). https://doi.org/10.1364/AO.40.000176 Google Scholar**

*Measurement of optical transport properties of normal and malignant human breast tissue***,” Numer. Heat Transfer**

*Discrete transfer method applied to transient radiative transfer problems in participating medium**A*, 44 183 –197 (2003). Google Scholar

**,” Appl. Math. Model., 10 107 –113 (1986). https://doi.org/10.1016/0307-904X(86)90080-6 Google Scholar**

*A boundary element method for the solution of the transient diffusion equation in two dimensions*