To make a tradeoff between computational efficiency and accuracy, a hybrid model of delicate transport physics and its diffusion approximation has attracted great attention in many applications that require more accurate characterization of near-field photon migration, such as mesoscopic and microscopic imaging, spectroscopy, and so on.¹ According to the coverage of the respective computational domains of the transport and diffusion theories, the established hybrid transport-diffusion models can be categorized into three implementations, as briefly discussed below.

The predominant method is to partition the whole computational domain ($\mathrm{\Omega }$) into two neighboring subdomains: the near-field region that necessitates the transport theory (${\mathrm{\Omega }}_{\mathrm{T}}$) and the far-field region that suits the diffusion approximation (${\mathrm{\Omega }}_{\mathrm{D}}$). In these hybrid models, the transport theory is in general described by the radiative transport equation (RTE), with its solution found either by using the discrete-ordinates method ($S_{\mathrm{N}}$)^2,3 or solving its high-order spherical harmonics (${\mathrm{P}}_{\mathrm{N}}/{\mathrm{SP}}_{\mathrm{N}}$, $\mathrm{N}\ge 3$) approximation.^4,5 The physical interdependence between the two subdomains causes neither of the two equations to be solved independently. By formulating the interface conditions, e.g., the continuity of the photon-density and its high-order derivatives (flux and so on) on the crossover interface, one straightforward scheme is to deal with a coupled set of the RTE and the diffusion equation (DE), abbreviated as the coupled RTE-DE model.^2–5 However, solving large-scale simultaneous equations might be mathematically intractable and numerically unstable, in particular for three-dimensional scenarios. In contrast, another roundabout but low-efficiency way takes an iteratively approximating strategy to solve the RTE and DE alternately until convergence of the successive approximations is reached.⁶

Some other divisions, referred to as the overlapped RTE-DE model, were first proposed by Degond and Jin.⁷ The scheme makes the subdomains overlap through an artificially defined buffer zone, over which a transition function is introduced that smoothly passes from 1 in ${\mathrm{\Omega }}_{\mathrm{T}}$ to 0 in ${\mathrm{\Omega }}_{\mathrm{D}}$ or vice versa to degenerate the RTE and DE at their respective ends of the buffer zone. Instead of using computationally intensive or incomplete interface conditions, the coupling inside the buffer-zone is naturally achieved through solving the two degenerated simultaneous equations. As an noniterative alternative to the simultaneous RTE and DE solutions, a buffered RTE-DE model was introduced, where the RTE solution is independently obtained in an empirically extended ${\mathrm{\Omega }}_{\mathrm{T}}$ to take into account the backscattered photons from ${\mathrm{\Omega }}_{\mathrm{D}}$ to ${\mathrm{\Omega }}_{\mathrm{T}}$. The DE in ${\mathrm{\Omega }}_{\mathrm{D}}$ is successively solved according to the Dirichlet interface condition defined by the previous RTE calculation.^8,9

Very recently, a whole-domain–based hybridization, known as the RTE-corrected DE model, has been proposed.^10,11 To suppress the propagation error of the ray effect in the RTE calculation that arises from discretizing the diffuse radiation in limited directions,¹² a more efficient and accuracy-retaining method is adopted in this scheme that uses the coarsely direction-discretized RTE for the non–fully diffusive part and the DE for the fully diffusive part. However, the modeling fails near the boundary¹¹ and might substantially lose the intrinsic advantage of hybrid strategies in computational efficiency due to the whole-domain RTE solution.

It is clear that the previously discussed hybrid models are limited to physical inconsistency and/or mathematical inefficiency due to the difficulty of explicitly decoupling the interdependence between ${\mathrm{\Omega }}_{\mathrm{T}}$ and ${\mathrm{\Omega }}_{\mathrm{D}}$, or the requirement of necessarily solving the whole-domain RTE. To combine the computational simplicity of the subdomain RTE solution and the physical rigorousness of the whole-domain photon-field decomposition, a novel hybrid scheme, abbreviated as the TPK model, is proposed here by analyzing the transient photon kinetics (TPK) in a turbid medium. In such a TPK scheme, the fully diffusive photon field at a transition time $t_{\mathrm{tr}}$ is determined by a computation-saving auxiliary time-domain (TD) DE solution and is utilized in a physically consistent way to decompose the whole-domain RTE solution into a temporally integrated subdomain RTE solution over the early-stage (before $t_{\mathrm{tr}}$) and a temporally integrated whole-domain DE solution over the late-stage (after $t_{\mathrm{tr}}$).

The proposed TPK strategy is shown schematically in Fig. 1, and will be theoretically interpreted below based on the RTE-DE hybridization. For the proof-of-concept survey, we considered a cubic and optically homogeneous geometry of $30l_t^{\prime }\times 30l_t^{\prime }\times 15l_t^{\prime }$, impinged by an infinitely thin collimated beam perpendicular to the boundary, i.e., along the $z$-direction, where $l_t^{\prime }=1/[{\mu }_{\mathrm{a}}+(1-g){\mu }_{\mathrm{s}}]$ is the transport mean free pathlength (mfp) with $g$, ${\mu }_{\mathrm{a}}$, and ${\mu }_{\mathrm{s}}$ being the anisotropy factor, absorption, and scattering coefficients, respectively, and no photons travel in an inward direction at the domain boundary $\partial \mathrm{\Omega }$. For the TPK analysis, the temporal dirac-delta source is considered, i.e., $q(\mathbf{r},\widehat{\mathbf{s}},t)=\delta (\mathbf{r}-{\mathbf{r}}_0,\widehat{\mathbf{s}}-\widehat{\mathbf{z}},t-t_0)$, where $\mathbf{r}$ is the positional vector in the medium with ${\mathbf{r}}_0=\{x_0,y_0,z_0\}$ denoting the incident position on the surface, $\widehat{\mathbf{s}}$ is the unit directional vector, $t_0(=0)$ is the time origin, and $\hat{\mathbf{z}}$ is the inward normal vector on the surface. According to the mean anisotropy degree of the incident photons, the transition time point can be defined as $t_{\mathrm{tr}}=k{l}_t^{\prime }/c$, with $k$ being the optical pathlength and $c$ the light speed in the turbid medium, after which the photon migration is considered to be fully diffusive.

Figure 1 calculates the transient photon-density distributions using the time-domain Monte Carlo (MC) simulation to demonstrate the four critical photon migration behaviors at the transition time point $t_{\mathrm{tr}}$ and during the early ($t_0\le t\le t_{\mathrm{tr}}$), late ($t_{\mathrm{tr}}\le t\le +\infty$) and whole ($t_0\le t\le +\infty$) stages, respectively. The whole-stage (i.e., the total steady state) photon-density $\mathrm{\Phi }(\mathbf{r})$ is in principle expressed, as shown in Figs. 1(b)–1(d), as a sum of the early- and late-stage portions, which are the temporal integrations of the TD photon-density $\mathrm{\Phi }(\mathbf{r},t)$ in the early and late stages, respectively, i.e.,

The first term on the right-hand side of the above equation can be expressed as ${\int }_{4\pi }{\varphi }^{\text{early}}(\mathbf{r},\widehat{\mathbf{s}})\mathrm{d}\widehat{\mathbf{s}}$, where the early-stage radiance ${\varphi }^{\text{early}}(\mathbf{r},\widehat{\mathbf{s}})$ is theoretically equal to the temporal integration of the time-domain RTE solution $\varphi (\mathbf{r},\widehat{\mathbf{s}},t)$ over time from $t_0$ to $t_{\mathrm{tr}}$, i.e., ${\varphi }^{\text{early}}(\mathbf{r},\widehat{\mathbf{s}})={\int }_{t_0}^{t_{\mathrm{tr}}}\varphi (\mathbf{r},\widehat{\mathbf{s}},t)\mathrm{d}t$, and can be calculated within the hemispheric near-field domain ${\mathrm{\Omega }}_{\mathrm{T}}$ with a radius of ${\mathrm{R}}_{\mathrm{T}}=c{t}_{\mathrm{tr}}$, surrounded by an interface Г inside the whole-domain $\mathrm{\Omega }$, outside which the photon concentration is essentially zero before $t_{\mathrm{tr}}$, as shown in Figs. 1(a) or 1(b). Consequently, the following steady-state RTE equation is derived:

with the physically rigorous boundary condition of for the modeling, where $p({\widehat{\mathbf{s}}}^{\prime },\widehat{\mathbf{s}})$ is the Henyey–Greenstein phase function specifying the scattering probability of a photon from ${\widehat{\mathbf{s}}}^{\prime }$ to $\widehat{\mathbf{s}}$ and $Q(\mathbf{r},\widehat{\mathbf{s}})$ is the distributed source term: Since the photons are fully diffused after $t_{\mathrm{tr}}$, $\varphi (\mathbf{r},\widehat{\mathbf{s}},t_{\mathrm{tr}})$ can be expressed as¹³ where $D=1/\{3[{\mu }_a+(1-g){\mu }_s]\}$ is the diffusion coefficient and $\mathrm{\Phi }(\mathbf{r},t_{\mathrm{tr}})$ can be obtained by solving an auxiliary TD-DE [${\mathrm{\Phi }}_{\mathrm{TD}-\mathrm{DE}}(\mathbf{r},t_{\mathrm{tr}})$] merely in ${\mathrm{\Omega }}_{\mathrm{T}}$, with the Dirichlet and Robin boundary conditions applied on the inner boundary Г and the outer boundary $\partial \mathrm{\Omega }$, respectively.¹³ Without loss of generality, $q(\mathbf{r},\widehat{\mathbf{s}},t)$ is converted into an isotropic point source at a typical depth of $l_t^{\prime }$ in the auxiliary TD-DE, with its intensity attenuated according to the Beer’s law:¹⁴

For the late-stage of $t_{\mathrm{tr}}\le t\le +\propto$, the integration term in Eq. (1), ${\mathrm{\Phi }}^{\text{late}}(\mathbf{r})={\int }_{t_{\mathrm{tr}}}^{+\infty }\mathrm{\Phi }(\mathbf{r},t)\mathrm{d}t$, can be similarly obtained by solving the temporally integrated DE over the late stage, expressed in the following steady-state form:

Finally, the total steady-state photon-density solved by the proposed hybrid model can be expressed by summing the two components:

The predefined optical pathlength $k$ is first evaluated for the typical optical properties of tissues in the visible or NIR wavelength range.¹⁵ Since the MC simulations have justified that, after the transition time point $t_{\mathrm{tr}}$, the light is considered fully diffused with its transient photon density nearly approaching the Gaussian typed distribution along the incident direction and centered at a fixed depth,¹⁴ it is reasonable in practice to calculate $t_{\mathrm{tr}}$ (or $k$) according to a critical time after which the photon density distribution centroid is unchanged. For the effective calculation of $t_{\mathrm{tr}}$, we introduce the temporal change rate of the transient depth $z_c(t)$ of the photon-density distribution centroid, i.e., ${\mathbf{\gamma }}_z(t)=\mathrm{d}{z}_c(t)/\mathrm{d}t$, as a quantitative criterion. To be concise, Fig. 2 illustrates that the MC calculated normalized ${\mathbf{\gamma }}_z(t)$ for only three different optical property sets, representing the normal (${\mu }_{\mathrm{a}}=0.01{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}=10{\mathrm{mm}}^{-1}$), low-scattering (${\mu }_{\mathrm{a}}=0.01{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}=1{\mathrm{mm}}^{-1}$), and high-absorption (${\mu }_{\mathrm{a}}=0.1{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}=10{\mathrm{mm}}^{-1}$) situations, respectively, all with $g=0.8$. It can be seen from Fig. 2 that all three curves tend to be stabilized as $k\ge 5$, where the change rate sharply drops from 0.7 at $k=1$ to $<0.25$ after $k\ge 5$. Consequently, the transition time point $t_{\mathrm{tr}}$ should be set to at least $5l_t^{\prime }/c$ in principle.

The calculation of ${\mathrm{\Phi }}_{\mathrm{TD}\text{-}\mathrm{DE}}(\mathbf{r},t_{\mathrm{tr}})$ is one crucial step in the TPK to achieve closed forms of Eqs. (2) and (7), respectively. The accuracy of ${\mathrm{\Phi }}_{\mathrm{TD}\text{-}\mathrm{DE}}(\mathbf{r},t_{\mathrm{tr}})$ is evaluated by comparing it with the TD-MC results through their relative deviations, ${ϵ}_{\mathrm{TD}\text{-}\mathrm{DE}}=|{\mathrm{\Phi }}_{\mathrm{MC}}-{\mathrm{\Phi }}_{\mathrm{TD}\text{-}\mathrm{DE}}|/{\mathrm{\Phi }}_{\mathrm{MC}}$, as shown in Fig. 3. As expected, for a small $t_{\mathrm{tr}}$ with $k=3$, the TD-DE calculation with the source in Eq. (6) leads to a considerable error of $\Vert {ϵ}_{\mathrm{TD}\text{-}\mathrm{DE}}\Vert =2.72$, where $||·||$ indicates the two-norm result for the $45\times 45$ matrix covering half of the $y\text{-}z$ cross-section through ${\mathbf{r}}_0$. Once it enters the fully diffusive stage, e.g., for the case of $k=5$, $\Vert {ϵ}_{\mathrm{TD}\text{-}\mathrm{DE}}\Vert$ sharply drops because of the dissolving and diffusive effects of the original impulse source; i.e., the exponential attenuation in the intensity and the rapid loss in the original direction of the photons instantaneously launched at ${\mathbf{r}}_0$, whose extents are essentially decided by the transition time $t_{\mathrm{tr}}$.

Although some analytic solutions to the RTE are available for optically homogeneous backgrounds,¹⁶ the TPK model is, in general, 3-D and can be numerically solved for arbitrary distributions and ranges of the optical properties, as well as for complex geometries. To be concise, the geometry configuration for the TPK model is the same as that for the MC in Fig. 1. For implementation herein, the finite-element method (FEM) and the ${\mathrm{S}}_{\mathrm{N}}\text{-}\mathrm{FEM}$ are employed to solve the DE in $\mathrm{\Omega }$, i.e., Eq. (6), for the late stage and the RTE in ${\mathrm{\Omega }}_{\mathrm{T}}$, i.e., Eq. (2), for the early stage, respectively. Owing to the DE description of the whole-domain photon migration in the late stage, a discrete-ordinates-method with fewer discretized directions, such as the ${\mathrm{S}}_{12}\text{-}\mathrm{FEM}$ scheme, accommodates the highly anisotropic scattering medium for early-stage modeling,^10,11 leading to a total of 168 directions, while the FEM calculation uses tetrahedral elements for the spatial discretization, with their number depending on ${\mathrm{R}}_{\mathrm{T}}$ for a nearly constant element dimension.

Without loss of generality, Fig. 4 compares the TPK model and MC solutions by calculating the relative deviations between the two models, ${ϵ}_{\mathrm{TPK}}=|{\mathrm{\Phi }}_{\mathrm{MC}}-{\mathrm{\Phi }}_{\mathrm{TPK}}|/{\mathrm{\Phi }}_{\mathrm{MC}}$, with different ${\mathrm{R}}_{\mathrm{T}}$, for the optical properties of normal tissue. For an overall assessment, $\Vert {ϵ}_{\mathrm{TPK}}\Vert$ and $\Vert {ϵ}_{\mathrm{DE}}\Vert$ are additionally calculated as a function of ${\mathrm{R}}_{\mathrm{T}}(k{l}_t^{\prime })$ in Fig. 4(c). Basically, in accordance with the evaluation of the optical pathlength $k$ and ${\mathrm{\Phi }}_{\mathrm{TD}\text{-}\mathrm{DE}}(\mathbf{r},t_{\mathrm{tr}})$, $\Vert {ϵ}_{\mathrm{TPK}}\Vert$ drops from 2.13 for $k=3$ to only 1.14 for $k=9$. The evident error for the first scenario can be attributed to the estimation of ${\mathrm{\Phi }}_{\mathrm{TD}\text{-}\mathrm{DE}}(\mathbf{r},t_{\mathrm{tr}})$ by incorrectly introducing the collimated component. Mathematically, the erroneous calculation of $\mathrm{\Phi }(\mathbf{r},t_{\mathrm{tr}})$ exerts adverse influence on the RTE [Eq. (2)] and DE [Eq. (6)] solutions through the source terms $Q(\mathbf{r},\widehat{\mathbf{s}})$ and $\mathrm{\Phi }(\mathbf{r},t_{\mathrm{tr}})$, respectively. Some other deviations within ${\mathrm{\Omega }}_{\mathrm{T}}$ are mainly caused by the error inherent to the RTE calculation, such as the ray effect, numerical smearing, and so on.¹² Moreover, difference in the incident collimation expressions between the MC and RTE calculations also contributes to the near-field deviation.

Performance comparisons among the TPK model and the four others, e.g., the coupled, overlapped, buffered, and corrected models, are listed in Table 1, where $\mathrm{\Delta }{\mathrm{R}}_{\mathrm{T}}$ is the required extension to ${\mathrm{R}}_{\mathrm{T}}$, ${\mathrm{\Omega }}_{\mathrm{T}}\cap {\mathrm{\Omega }}_{\mathrm{D}}$ indicates the overlapping domain of the RTE and DE solutions, and ${\mathrm{\Omega }}_{\mathrm{Err}}$ denotes the domain probably contaminated by the numerical error of the RTE solution. It can be seen that ${\mathrm{R}}_{\mathrm{T}}$ is considerably extended in both the overlapped and buffered models, and therefore greatly degrades their computational efficiencies. Since the ray-effect error in the RTE solution might be suppressed to some extent by leaving the fully diffusive component in ${\mathrm{\Omega }}_{\mathrm{T}}$ to be calculated by the DE,^10,11 the coverage of ${\mathrm{\Omega }}_{\mathrm{T}}\cap {\mathrm{\Omega }}_{\mathrm{D}}$ is compared for the five hybrid models in the table. In the TPK model, due to the late-stage DE calculation for the whole-domain fully diffusive photons, the boundary-value errors are alleviated even with the ${\mathrm{S}}_{12}\text{-}\mathrm{RTE}$, as shown in Fig. 4(b). As indicated in Table 1, ${\mathrm{\Omega }}_{\mathrm{Err}}$ can go beyond ${\mathrm{\Omega }}_{\mathrm{T}}$ through the interdependency between the DE and RTE solutions, resulting in error propagation and amplification in to be subdomains. This can eventually cause the coupled RTE-DE modeling for the near-field even inferior to that of the DE once ${\mathrm{R}}_{\mathrm{T}}$ is less than $\sim 8.5l_t^{\prime }$, as demonstrated in Ref. 2. In contrast, this adversity of error propagation can be avoided in the corrected and TPK models; e.g., the TPK outperforms the DE with the smaller ${\mathrm{R}}_{\mathrm{T}}$ of $5l_t^{\prime }$ or $7l_t^{\prime }$, as shown in Fig. 4.

Table 1

Performance comparisons among the TPK and the other four hybrid RTE-DE models.

In conclusion, we develop a physically consistent novel RTE-DE coupling scheme based on the TPK. This scheme decomposes steady-state photon transport modeling into a near-field localized RTE solution and a full-domain DE solution for the early and late stage, respectively. By properly defining the transition time at which the transient photon-density distribution can be found simply by solving the TD-DE, the early-stage RTE and the late-stage DE can be solved independently for an arbitrary distribution of the optical properties and generalized to complex geometries. Future work will focus on applying the proposed method in experiments and combining the RTE with the other higher-order diffusion approximations (${\mathrm{P}}_{\mathrm{N}}/{\mathrm{SP}}_{\mathrm{N}}$) to further promote computational efficiency with a reduced dimension (e.g., ${\mathrm{R}}_{\mathrm{T}}$) for near-field modeling.