On the geometry dependence of differential pathlength factor for near-infrared spectroscopy. I. Steady-state with homogeneous medium

Abstract. This work analytically examines some dependences of the differential pathlength factor (DPF) for steady-state photon diffusion in a homogeneous medium on the shape, dimension, and absorption and reduced scattering coefficients of the medium. The medium geometries considered include a semi-infinite geometry, an infinite-length cylinder evaluated along the azimuthal direction, and a sphere. Steady-state photon fluence rate in the cylinder and sphere geometries is represented by a form involving the physical source, its image with respect to the associated extrapolated half-plane, and a radius-dependent term, leading to simplified formula for estimating the DPFs. With the source-detector distance and medium optical properties held fixed across all three geometries, and equal radii for the cylinder and sphere, the DPF is the greatest in the semi-infinite and the smallest in the sphere geometry. When compared to the results from finite-element method, the DPFs analytically estimated for 10 to 25 mm source–detector separations on a sphere of 50 mm radius with μa=0.01  mm−1 and μs′=1.0  mm−1 are on average less than 5% different. The approximation for sphere, generally valid for a diameter ≥20 times of the effective attenuation pathlength, may be useful for rapid estimation of DPFs in near-infrared spectroscopy of an infant head and for short source–detector separation.


Introduction
Near-infrared spectroscopy (NIRS) 1 is becoming an increasingly important modality for a number of investigational and clinical needs, including noninvasive functional study of neurophysiological connectivity, 2 bedside monitoring of cerebral hemodynamics, 3 and evaluating responses to cancer treatment. 4 Most NIRS systems for long duration monitoring operate under the continuous-wave (CW) or steady-state condition. 5 Quantitation of the changes of chromophores, most importantly oxyhemoglobin and deoxyhemoglobin, within the highly diffusing tissue typically relies upon a priori estimates of the effective pathlength of the photons that propagate from a source to a detector. Delpy et al. 6 first suggested utilizing the ratio of the effective pathlength of the photon trajectory over fixed separations of source-detector pairs as the differential pathlength factor (DPF) in the routine for quantitative NIRS. A DPF value of 5.3 AE 0.3 across a rat head was experimentally obtained by using time-of-flight techniques. Although the DPF was observed to be approximately constant for source-detector spacings greater than 2.5 cm, frequency-domain NIRS studies by Duncan et al. 7 at source-detector spacings >4 cm yielded DPFs of 6.26 (AE14.1%) for an adult head, 5.51 (AE18%) for an adult leg, 4.99 (AE9%) for the head of a newborn infant, and 4.16 (AE18.8%) for an adult arm. These experimentally obtained values indicate that the DPF may depend on the dimensions and shape of the tissue structure under investigation, providing the same source-detector separations (fixed physical configuration) and tissue optical properties (likely subject condition).
In terms of the analytical form of DPF, Delpy et al. 6 suggested a wavelength-, medium-, and geometry-dependent constant 8 in addition to the term straightforwardly derived from the solution to photon diffusion in an unbounded homogeneous medium. In the case of a homogeneous semi-infinite medium, explicit analytical formation of the DPF was given by reports including Fantini et al. 9 and Boas et al. 10 A semi-infinite geometry is a reasonable representation of the measurement geometry for a set of optodes separated by a few centimeters when externally placed on large tissue areas. However, in cases involving measurement of the brain, different head sizes (infant to adult), different global shapes of the head (Caucasian to Asian), different regions of the head 11 (e.g., prefrontal or parietal), and even different ages 12 may result in different DPFs. Therefore, an analytical form of DPF that indicates the effect of the size and shape of the nonsemi-infinite geometry, even at the simplest case of a homogeneous medium that is the basis for analysis of more realistic geometries, may be informative for the development of better strategies to the quantitation of the chromophore changes when the wavelength dependency of DPF is also taken into account. 13 In addition, the more accurate a priori estimation of the DPF for the subsets of optode configurations with shorter source-detector separations, the more reliable the determination of the shallow-layer contamination to the deeper-layer signal. 14 Whereas estimates of DPFs from complex geometries are readily available with the use of discrete numerical methods [e.g., finite-element method (FEM) [15][16][17] ], access to analytical estimators of DPFs that can be quickly implemented with reasonable tolerance can hold considerable practical and clinical value, especially in situations where near immediate or realtime quantitates of DPFs are required (e.g., acute care environments) to allow timely diagnosis and intervention. Mainly lacking in support of this, an analytically rendered understanding of the DPF in geometries more complex than the semi-infinite case, specifically how the shape and dimension of a spherically curved geometry resembling an infant head affect the DPF, is relevant to NIRS applications for better clinical care.
Mathematical treatments of analytical solutions to the photon diffusion equation (DE) in a medium bounded externally by a circular cylindrical applicator or a spherical applicator have been investigated over a course of nearly two decades in many comprehensive studies. [18][19][20][21][22] A particularly complete set of solutions in time domain and frequency domain in a number of geometries were presented by Arridge et al. 18 in 1992; the DPFs could be obtained from the explicit solutions derived for the mean time of flight in a variety of geometries using the zero boundary condition. In 2001, Sassaroli et al. 22 developed the corresponding solutions of time-resolved DE using the more widely used extrapolated boundary condition (EBC) 21,23,24 for curved geometries including infinite cylinder and sphere; the DPFs for infinite cylinder and sphere can also be derived or computed from the results that were in excellent agreement with Monte Carlo (MC) simulations. While these studies as represented by Refs. 18 and 22 have theorized the computational availability of accurate time-domain DPFs for curved geometries of interest, the computational complexity arising from applying high-precision arithmetic library on infinite-series-summations involving various special functions renders them out of reach for most members of the NIRS community. Solutions for the CW reflectance and subsequently the DPF that preferably involve only the regular functions and do not impose lengthy summations, while perhaps less accurate, may be much more convenient for most clinically oriented CW NIRS applications in rapidly estimating the DPF for timely physiological assessment. This present study aims to develop the first steps toward this unmet need.
In our recent approach to a new routine of solutions of DE for an infinite-length cylinder domain as demonstrated in a series of work applicable to CW, 25 frequency domain 26 and time domain, 27 the solutions for DE were presented using a pair of the modified Bessel functions of the first and the second kinds. The basic form of the DE solution for the unbounded homogeneous medium geometry shown in Refs. 25-27 may be similar to the one used by Ref. 28; however, our method of utilizing the EBC in the case of the medium bounded within an infinitelength cylinder was new in terms of leading to a more intuitive interpretation of the solution. The solution obtained for the infinite-length cylinder domain is presented to consist of two parts: the first part is associated with the "real" isotropic source that is the common approach of treating the light normally incident to the medium from the source fiber or channel, and the second part accounts for the contribution of the "image" source as a result of the medium boundary that is represented by the "real" source term scaled by a factor determined by the radius of the cylinder. This solution is subsequently formatted to a form similar to that for the semi-infinite geometry except a shape-curvature-associated term that makes the solution for the cylinder domain approaching the one for semi-infinite geometry as the radius of the cylinder domain increases without bound. Such treatment allowed obtaining approximating solutions of CW DE in cylinder geometries of relatively large radius using functions more convenient to evaluate than special functions and without involving the infinite series of summations. Additionally, the analytical methodology explored in Refs. 25-27 applies to two opposite medium geometries, one for the medium enclosed by the optode applicator as with the arm or leg, and the other for the medium enclosing the optode applicator as with the prostate via transrectal probing.
In continuation of our demonstrated approaches to the homogeneous infinite cylinder domain, we develop in this present work the solution of steady-state DE in a homogeneous spherical domain (limited to the external applicator case although the approach is extendable to the internal applicator case) to a form that is also composed of two parts: the first part is associated with the "real" isotropic source, and the second part is the contribution of the "image" source term that is represented by the "real" source term scaled by a factor depending upon the radius of the sphere. The format of the solution for the spherical domain, or the subsequent approximation, shall also be similar to that for the semi-infinite geometry except for a shape-curvature-associated term, which nonetheless must be different from its counterpart for the infinite cylinder, to make the solution for the sphere domain approaching the one for semi-infinite geometry, as the radius of the sphere increases without bound.
With the solutions of steady-state DE associated with the geometries of interest becoming available, the DPF can then be evaluated by calculating the partial derivative of the attenuation of the photon fluence with respect to the absorption and dividing that by the source-detector separation. The usefulness of the quantitative patterns of analytically approximated DPF associated with curved geometries can then be examined by using numerical simulation, which in this work is limited to FEM method.

Geometries Considered and Analytic Preparations for Steady-State Photon Diffusion in Infinite and Semi-Infinite Geometries
We consider a scattering domain with an absorption coefficient of μ a and a diffusion coefficient of D ¼ 1∕ð3μ 0 s Þ where μ 0 s is the reduced scattering coefficient. It must be noted that it is yet to reach a consensus regarding whether the diffusion coefficient in steady-state measurement shall include the absorption coefficient. [29][30][31][32][33] We argue that, as the DPF is unbiasedly defined by the use of mean time of flight through time-dependent DE, the definition of DPF in steady-state may agree with the one in the time-resolved condition only if the diffusion coefficient is not dependent upon the absorption coefficient. Therefore, this study uses the diffusion coefficient as D ¼ 1∕ð3μ 0 s Þ. This study intends to analytically compare the DPF in three geometries of the bounded medium that include a semi-infinite geometry, an infinite-length cylinder geometry, and a spherical geometry. Since the solutions for bounded geometries originate from the solution for the infinite medium geometry, it becomes imperative to establish a general analytical routine to the studied geometries as are shown in Fig. 1, and from which to evaluate the DPF of steady-state measurement. For obvious conveniences, the infinite and semi-infinite domains will be represented using the Cartesian coordinates, where we are concerned about photon pathlength from a source atχ 0 ¼ ðx 0 ; y 0 ; z 0 Þ to a detector atχ ¼ ðx; y; zÞ. The infinite-length cylinder domain will be represented with the polar cylindrical coordinates, where we are concerned about photon pathlength from a source at χ 0 ¼ ðρ 0 ; φ 0 ; z 0 Þ to a detector atχ ¼ ðρ; φ; zÞ. The sphere domain will be modeled in the spherical coordinates, where we are concerned about photon pathlength from a source atχ 0 ¼ ðr 0 ; θ 0 ; φ 0 Þ to a detector atχ ¼ ðr; θ; φÞ.

Infinite Geometry and Semi-Infinite Geometry
A homogeneous medium of infinite geometry is illustrated in Fig. 1(a), wherein the source is regarded as being isotropic. In the infinite geometry, it is convenient to set the source at χ 0 ¼ ð0;0; 0Þ, and the detector atχ ¼ ð0; d; 0Þ, where d is the source-detector separation.
A homogeneous medium of semi-infinite geometry with the physical source and the detector located on the medium boundary is depicted in Fig. 1(b). In the semi-infinite geometry, it becomes convenient to set the physical source atχ 0 ¼ ð0;0; 0Þ, and the detector atχ ¼ ð0; d; 0Þ. The physical source at χ 0 ¼ ð0;0; 0Þ launches photon into the medium at an initial direction orthogonal to the medium-applicator interface, and the resulted photon diffusion is treated as being excited by an equivalent "real" isotropic source located at one step of reduced scattering, represented by R a ¼ 1∕μ 0 s , into the medium. The "real" isotropic sourceχ 0 real thus has the coordinates of ð0;0; −R a Þ. The effect of the medium-applicator interface on photon diffusion as modeled by the EBC 21,23,24 sets zero the photon fluence rate on an imaginary boundary located at a distance of R b ¼ 2AD away from the physical boundary, where A is a coefficient 23 related to the refractive index differences across the physical boundary. This EBC is conventionally accommodated by setting a sink or a negative "image" of the "real" isotropic source, with respect to the extrapolated boundary. The "image" source for the semi-infinite geometry has the opposite strength of the "real" isotropic source and locates atχ 0 imag , which has the coordinates of ð0;0; R a þ 2R b Þ. The distances from the detector to the "real" isotropic sourceχ 0 real and the "image" sourceχ 0 imag are denoted by l r and l i , respectively. The notations of l r ¼ jχ −χ 0 real j and l i ¼ jχ −χ 0 imag j also apply to other Fig. 1 Four cases of medium geometry. (a) Infinite geometry, where an isotropic source is at χ 0 ¼ ð0;0; 0Þ and a detector is atχ ¼ ð0; d; 0Þ. (b) Semi-infinite geometry wherein the source and the detector locate atχ 0 ¼ ð0;0; 0Þ andχ ¼ ð0; d; 0Þ, respectively. The equivalent isotropic source locates inwardly with respect to the interface at ð0;0; −R a Þ, and the image of it with respect to the boundary locates at (0; 0; R a þ 2R b ). (c) Infinite-length cylindrical geometry, where the source and the detector locate at ðR 0 ; φ 0 ; z 0 Þ and ðR 0 ; φ; zÞ, respectively. The equivalent isotropic source locates inwardly with respect to the curvature of the interface at ðR 0 − R a ; φ 0 ; z 0 Þ, and the image of it with respect to the boundary locates at a to-be-determined radial position of ðρ ? ; φ 0 ; z 0 Þ, which will reach ðR 0 þ R a þ 2R b ; φ 0 ; z 0 Þ as the radius of the cylinder increases toward infinity. (d) Spherical geometry, where the source and the detector locate at ðR 0 ; θ 0 ; φ 0 Þ and ðR 0 ; θ; φÞ, respectively. The equivalent isotropic source locates inwardly with respect to the curvature of the interface at ðR 0 − R a ; θ 0 ; φ 0 Þ, and the image of it with respect to the boundary locates at a to-be-determined radial position of ðr ? ; θ 0 ; φ 0 Þ, which will reach ðR 0 þ R a þ 2R b ; θ 0 ; φ 0 Þ as the radius of the cylinder increases toward infinity. medium geometries studied in this work when involving a boundary. In all studied geometries, the straight distance between the physical source and the detector d ¼ jχ −χ 0 j is referred to as the "line-of-sight" source-detector distance, and hereafter specified as source-detector distance. It is with respect to this d, which is to remain identical across different medium geometries considered in this work, that the DPF is to be expressed, facilitating the evaluation of the effect of shape and dimension of the nonsemiinfinite-domains on DPF in comparison to the semi-infinite geometry on DPF.

Infinite-Length Cylindrical Geometry
The homogeneous medium bounded in an infinite-length cylindrical geometry is illustrated in Fig. 1(c). The radius of the cylinder is R 0 , then a detectorχ on the medium-applicator interface locates at ðR 0 ; φ; zÞ, and a directional sourceχ 0 on the mediumapplicator interface locates at ðR 0 ; φ 0 ; z 0 Þ. The geometric symmetry requires that the directional sourceχ 0 be modeled by a "real" isotropic sourceχ 0 real positioned along the radial direction ofχ 0 and inwardly at a radial distance of R a ¼ 1∕μ 0 s , at ðR 0 − R a ; φ 0; z 0Þ. The extrapolated zero-boundary is concentric with and at a radial distance of R b ¼ 2AD outward from the physical boundary. Apparently, as the radius R 0 of the infinite-length cylinder reaches infinity, the photon fluence rate associated with a source-detector pair on the infinite-length cylindrical medium-applicator interface shall reach that on the semi-infinite medium interface for the same source-detector distance. This characteristic has served for "checking" the results developed for infinite-length cylindrical geometry. 25 If one considers reaching the curved geometry from a semiinfinite geometry that extends to infinities at two orthogonal dimensions of the direction, bending on one dimension of the direction will make the effect of an infinite-length cylinder, and concurrent bending on both dimensions of the direction will reach the effect of a sphere. From such an intuitive perspective, we could expect that the effect of the one-dimensional (1-D) bending of the semi-infinite domain on photon diffusion will be seen when evaluating the photon diffusion along the azimuthal direction of the resulted infinite-length cylinder geometry, and that effect on photon diffusion when compared to semiinfinite geometry will be less than the double-dimensional bending of the semi-infinite domain to reaching the sphere for the same source-detector distances. Consequently, we evaluate the infinite-cylinder domain in the azimuthal direction only.

Spherical Geometry
A spherical geometry is illustrated in Fig. 1(d). The radius of the spherical applicator is R 0 , then a detectorχ on the mediumapplicator interface locates at ðR 0 ; θ; φÞ, and a directional sourcẽ χ 0 on the medium-applicator interface locates at ðR 0 ; θ 0 ; φ 0 Þ. The geometric symmetry requires that the directional sourcẽ χ 0 be modeled by a "real" isotropic sourcer 0 real positioned along the radial direction ofχ 0 and inwardly at a radial distance The extrapolated zero boundary is concentric with and at a radial distance of R b ¼ 2AD outward from the physical boundary. Apparently, as the radius R 0 of the spherical geometry reaches infinity, the photon fluence rate associated with a source-detector pair on the spherical medium-applicator interface shall reach that on the semi-infinite medium interface for the same source-detector distance. It shall be noted, however, that the rate at which the photon fluence approaches the value of the semi-infinite geometry will be different in spherical geometry comparing to the infinite-length cylindrical geometry, given the same radius.

Steady-State Photon Diffusion in an Infinite
Homogeneous Domain The equation of steady-state photon diffusion in a homogeneous medium is expressed by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 6 2 7 where Ψ is the photon fluence rate and S is the source term.
is the effective attenuation coefficient. The solution to Eq. (1) also can be written in the form of the following eigenfunction expansion (for derivation, see Appendix A): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 4 5 3 where i l and k l are respectively the l'th order modified spherical Bessel function of the first and the second kinds, r < and r > are respectively the smaller and greater radial coordinates of the source and detector, and Y lm is the spherical harmonics function. The numerical identity between Eq. (3) and the well-known form of Eq. (2) is justified in Appendix B.

Steady-State Photon Diffusion in a Homogeneous Semi-Infinite Medium
We follow the notations summarized by Fantini et al. 34 for the semi-infinite geometry having a directional source and an isotropic detector located on a planar boundary, as illustrated in Fig. 1(b). As stated in Sec. 2.1, the directional source is modeled as an isotropic source placed one reduced scattering distance into the medium. Then, based on the extrapolated boundary and the "image" source approach, the steady-state photon fluence rate reaching the detector located on the physical boundary is determined by the equivalent "real" isotropic source and the "image" of it with respect to the extrapolated boundary. That gives E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 1 5 3 where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 9 8 l r ¼ d

Steady-State Photon Diffusion in a Homogeneous Medium Bounded in a Spherical Geometry
The principles to find the solutions to steady-state photon diffusion associated with a homogeneous infinite-length cylinder 25 will not be duplicated here. The solutions associated with infinite-cylinder geometry will however be presented when needed for the comparison with the solutions obtained for other geometries.

Solution of Steady-State Photon Diffusion Associated with a Homogeneous Medium of Spherical Geometry
In surface imaging of a medium of spherical geometry with a radius of R 0 , the physical source is located at ðR 0 ; θ 0 ; φ 0 Þ and the detector is located at ðR 0 ; θ; φÞ, both on the physical boundary. The equivalent "real" isotropic source locates at ðR 0 − R a ; θ 0 ; φ 0 Þ, and the extrapolated boundary locates at a radial position of ρ r> ¼ R 0 þ 2AD outside the physical boundary. The geometric symmetry determines that the image of the "real" isotropic source with respect to the extrapolated boundary must locate along the radial direction of the "real" isotropic or the physical source. This image source and the "real" isotropic source collectively set zero the photon fluence rate on the extrapolated boundary. Based on Eq. (3), the photon fluence rate associated with the "real" isotropic source and evaluating on the extrapolated boundary, for which the source locates at r < ¼ R 0 − R a and the detector locates at ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; 3 6 3 where the notation "leftjright" indicates evaluating the "left" as a source on the "right" as a boundary. Note that the source term "S" is placed inside the summation. This basically indicates that any "l"'th order of the "real" source has the same intensity of "S,",as compared to different intensity for each order of the "image" source as will be shown in the following equation. Similarly, the photon fluence rate associated with the image source and evaluating on the extrapolated boundary, for which the source locates at radial position of a to-be-determined r > and the detector locates at ; t e m p : i n t r a l i n k -; e 0 0 8 ; 6 3 ; 1 7 0 where the S Ã l terms associated with each "l"'th order are not the same. Based on the essence of "image-source," 35,36 the two unknown terms S Ã l and r > associated with the l'th order "image" source (the k l component) can be expressed by a single unknown term S l associated with the same order "real" source (the i l component), that is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 3 2 6 ; 7 1 2 Applying Eq. (9) to the EBC of Ψ real j extr − Ψ imag j extr ¼ 0, we have E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 3 2 6 ; 6 5 9 Now for the "real" source but evaluating at the physical boundary, the source still locates at r < ¼ R 0 − R a but the detector locates at r > ¼ R 0 . For the "image" source also evaluating at the physical boundary, the detector locates at r < ¼ R 0 and the source terms are now known through Eqs. (9) and (10). Collectively, the photon fluence rate sensed by a detector at the physical boundary becomes E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 3 2 6 ; 5 3 7 Equation (11) is convertible to the following form that may be more convenient for numerical evaluation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 3 2 6 ; 2 7 1 where I l and K l are respectively the l'th order modified Bessel function of the first and the second kinds, P l is the l'th order Legendre polynomial, and γ is the angle formed between the position vector of the source and that of the detector (see also Appendix A). The value determined by Eq. (11) or Eq. (12) is compared against FEM results (configuration of FEM will be detailed in subsequent sections), as shown in Fig. 2, for two radii of 25 mm and 50 mm at three sets of optical properties: the top traces correspond to μ a ¼ 0.01 mm −1 and μ 0 s ¼ 0.25 mm −1 , the middle traces μ a ¼ 0.01 mm −1 and μ 0 s ¼ 0.50 mm −1 , and the bottom traces μ a ¼ 0.01 mm −1 and μ 0 s ¼ 1.0 mm −1 . A total of 51 detector positions are placed at a source-detector distance of 1 to 50 mm at a step of 1 mm from a single source. At the radius of 25 mm, the farthest detector point is positioned at the opposite to the source point, and at the radius of 50 mm, the farthest detector point is positioned at 60 deg with respect to the source point. The results are presented in terms of lnðd 2 Ã ΨÞ versus d, as it leads to a convenient near-linear relationship as the source-detector distance becomes longer than approximately five times of the reduced scattering pathlength. The oscillation observed of the analytical result corresponding to μ 0 s ¼ 1.0 mm −1 is due to arithmetic overflow in computing the modified Bessel functions. It is shown that at source-detector separations less than approximately five times of the reduced scattering pathlength (i.e., in the subdiffusion regime) there is a discrepancy between the analytically evaluated photon fluence rate and its FEM counterpart. As observed in our previous study 37 dedicated to experimental examination of the analytical results of photon fluence rate specific to infinite-length cylinder, 25 the discrepancy between the diffusion-based analytical quantification and the FEM simulation in the subdiffusion regime was largely attributed to the difference between a spatially impulsive source treated in the analytical solution and a spatially distributed Gaussian source profile utilized in the FEM tool available for this work. The spatially impulsive source treated in the analytical solution was physically analogous to the single-point launching of forward directional photons into the medium, whereas the spatially distributed Gaussian source profile implemented in the FEM had the size effect manifested increasingly at shorter distances from the source. Additionally, the treatment to boundary effect (e.g., the probability of a boundary-reaching photon returning to the medium) in the analytical model was based on EBC, whereas a Robin-type boundary condition was programmed in the FEM; the subtle difference in the boundary conditions could have been amplified at the subdiffusion regime. In the diffusion regime, however, the analytical computations marked by the solid lines and the FEM results given in discrete markers clearly match, corroborating that Eq. (11) is the solution to the steady-state photon diffusion in a homogeneous medium bounded in a sphere.

Large-Diameter Approximation of the Solution for a Spherical Domain
As shown in Fig. 1(d), if a plane tangential to the sphere at the physical source position is considered as an imaginary planar boundary to a semi-infinite geometry of the medium to which the medium of the spherical geometry will approach as the radius of the sphere increases without bound, then the "real" isotropic source in this semi-infinite geometry still locates at ðR 0 − R a ; θ 0 ; φ 0 Þ, but the image of the "real" isotropic source with respect to this "hypothetical" semi-infinite boundary will be at ðR 0 þ R a þ 2R b ; θ 0 ; φ 0 Þ. According to Eq. (3), the photon fluence rate sensed by a detector on the spherical boundary due to the image of the "real" isotropic source associated with this "hypothetical" semi-infinite boundary is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 3 2 6 ; 5 7 3 The photon fluence rate sensed by a detector on the spherical boundary due to the image of the "real" isotropic source associated with the spherical boundary, as seen in Eq. (11), can be rewritten to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 3 2 6 ; 4 4 6 where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 3 2 6 ; 3 5 3 Note that R 0 ≫ R a , R b and if the argument of the modified Bessel functions is much greater than 1 (i.e., μ eff R 0 ≥ 10 that corresponds to R 0 ≥ 57.8 mm for μ a ¼ 0.01 mm −1 and   15) can be simplified by their asymptotic expressions, 38 then Eq. (15) becomes E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 6 3 ; 7 0 1 η l ¼ Comparing to its counterpart for infinite cylinder [Eq. (3.1.9) in Ref. 25], Eq. (16) for sphere geometry is different only in the power of the entire fraction, that is, one or one-half-an intuitive analogy of the double-sided bending or single-sided bending from a semi-infinite geometry. Substituting Eq. (16) into Eq. (14) and comparing to Eq. (13) we have E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 6 3 ; 5 9 2 Ψ imag j phys ¼ Ψ semi imag j phys · Hence, for a sphere of large diameter, Eq. (12) approximates to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 6 3 ; 5 2 The Ψ real j phys of Eq. (18) is essentially the same as the Ψ real in Eq. (4). As the radius of the sphere reaches infinity without bound, we have ðR 0 þ R a þ 2R b Þ∕ðR 0 − R a Þ → 1, and the Ψ semi imag j phys of Eq. (18) will become the Ψ imag in Eq. (4) because the detector point locating at ðR 0 ; θ; φÞ will reach the "hypothetical" semi-infinite boundary. By using the expression of the photon fluence rate given in Eq. (2), we can rewrite Eq. (18) as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 6 3 ; 3 7 For the purpose of assessing the general patterns of the steady-state photon fluence rate versus the source-detector distance, the following expressions of the involved distances will be implemented in Eq. (19): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 6 3 ; 2 8 4 l r ¼ d E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 6 3 ; 2 3 8 l i ¼ d The formula for the photon fluence rates associated with the four geometries of consideration, including infinite medium, semi-infinite medium or half-plane, large-radius infinite-cylinder when evaluated along the azimuthal direction, and largeradius sphere is tabulated in Table 1. The numerical results shown in the later sections performed for conditions including μ eff R 0 ¼ 8.66 will indicate that the approximate formula for the infinite cylinder and the sphere is generally valid when the radius satisfies μ eff R 0 ≥ 10, that is, the radius (diameter) is more than 10 (20) times of the effective attenuation pathlength 1∕μ eff .

Differential Pathlength Factors Associated with Homogeneous Semi-Infinite, Infinite-Length Cylinder When Evaluated along the Azimuthal Direction and Spherical Geometries
The DPF was originally introduced in Ref. 6 to relate the changes of the light to the changes of the absorption coefficient, at a given source-detector pair, by the modified Beer-Lambert law of Att ¼ exp½−μ a · DPF · d, where "Att" is the ratio of measured light over the incident light. We use the definition of DPF provided by Fantini et al. 9 or Arridge et al. 18  E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 3 2 6 ; 3 8 9 • Semi-infinite geometry [identical to Eq. (9) of Ref. 9] E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 3 ; 3 2 6 ; 3 2 0 • Infinite-length cylindrical geometry, evaluated along the azimuthal direction E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 4 ; 6 3 ; 7 1 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 ; 6 3 ; 5 7 1 where η azi ¼ • Spherical geometry E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 ; 6 3 ; 5 1 0 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 7 ; 6 3 ; 3 7 5 where η sph ¼ Note that Eq. (23) is the same DPF for semi-infinite medium as was used by Fantini et al., 9 and both Eqs. (24) and (26) will reach Eq. (23) as the radius R 0 becomes infinity.

Numerical Evaluation of the Differential Pathlength Factors Associated with Homogeneous Semi-Infinite, Infinite Cylinder Evaluated along the Azimuthal Direction and Spherical Geometries
The patterns of DPF among the three specific configurations associated respectively with semi-infinite geometry, infinitelength cylinder geometry when evaluated along the azimuthal direction, and sphere were examined by using near-infrared frequency-domain absorption and scattering tomography (NIRFAST), 39,40 an FEM platform for solving photon diffusion problems. A half-plane for semi-infinite geometry and an infinite-length cylinder cannot be produced by FEM, but could be reasonably represented in FEM if the length dimension remained substantially greater than the source-detector spacings. The cylinder and the sphere domains were examined at two relatively small radii of R 0 ¼ 25 mm and R 0 ¼ 50 mm with consistent node density of 1 mm that was comparable to the reduced scattering pathlength ð1∕μ 0 s Þ. The mesh for the semi-infinite geometry was generated for a 100 mm × 100 mm surface area and 50 mm depth, the mesh for the infinite-length cylinder geometry at a 100 mm length with the aforementioned 25 mm or 50 mm radius, and the mesh for the sphere at the same set of radius. Figure 3 displays the meshes (cylindrical and spherical meshes at 25 mm radius only).
In the configuration for the semi-infinite geometry, one source point was set in the middle-sagittal plane and at 30 mm away from one edge of the domain. A total of 50 detector points at an interval of 1 mm were placed in the middle-sagittal plane to cover a source-detector range from 1 mm to 50 mm. In the configurations for the sphere geometry and the cylinder geometry but evaluating along the azimuthal direction, one source point was set at the middle azimuthal plane, and a total of 50 detector points were placed in the middle azimuthal plane to have a range of the source-detector distance ("line-ofsight" distance or the chord distance) from 1 mm to 50 mm at an interval of 1 mm.
The DPFs were calculated using the forward solver of NIRFAST by the following steps: (1) the first set of the photon fluence data Ψ 1 measured at the 50 detector positions were forward calculated at the following setting of the optical properties: The second set of the photon fluence data Ψ 2 measured at the same set of 50 detector positions were forward calculated with the μ a at 0.011 mm −1 and other parameters remaining unchanged from the values of the previous step. (3) The DPF was calculated as f½lnðΨ 1 ∕Ψ 2 Þ∕Δμ a g∕d. Note that NIRFAST by default includes μ a in the calculation of the diffusion coefficient. To make the FEM evaluation as consistent with the analytical approaches as would be available, the μ a has been excluded from the calculation of the diffusion coefficient when using NIRFAST through the workspace open to the users. The DPFs for the source-detector distances ranging from 1 mm to 50 mm in the three specified configurations are shown in the left column of Fig. 4.
The DPFs evaluated according to the set of Eqs. (23), (24), and (26) are given in the right column of Fig. 4. The equations were evaluated with the following parameters: μ a ¼ 0.01 mm −1 , μ 0 s ¼ 1.0 mm −1 , R 0 ¼ 25 mm, and R 0 ¼ 50 mm over the source-detector distance from 3 mm to 50 mm. From both columns of Fig. 4, one may observe the following qualitative patterns of DPF at the same source-detector separation: (1) the DPF evaluated on a semi-infinite medium boundary is greater than the DPF evaluated along the azimuthal direction of the infinite-length cylindrical medium boundary; (2) the latter DPF of (1) is greater than the DPF evaluated on a spherical boundary, given the same radial dimension for the cylinder and the sphere. It is also obvious from Eqs. (23), (24), and (26) that as R 0 → ∞, Eqs. (24) and (26) will all reach Eq. (23), whereas at different rates. The inset at the top row illustrates how the DPF of the sphere calculated using Eq. (26), which is an approximation, fairs with the one calculated with FEM. The analytically approximated DPF for the spherical geometry of 50 mm radius at the given set of optical properties (resulting in μ eff R 0 ¼ 8.66) reached an average error of not more than 5% over the 10 mm to 25 mm source-detector distances when compared with the FEM results performed for the same spherical geometry. If the analytic DPF of the semi-infinite geometry were to be compared to, the DPF approximated for the spherical geometry of 50 mm radius at the set of optical properties as used for Fig. 4 has shown a maximum difference of about 9% at a source-detector distance of 9 mm as shown in Fig. 5. The relative difference of the DPF approximated for the spherical geometry with respect to the analytic semi-infinite one keeps reducing as the radius of the sphere increases-less than 5% at a radius of 110 mm for the given set of optical properties.

Discussions
The dependence of DPF on age has been well documented. In 1996, Duncan et al. showed that the older the subject was the larger the DPF became. 12 At 832 nm, the DPFs measured by using phase-resolved NIRS from the head of neonates were 4.67 AE 0.65, while the DPFs of adults were 5.86 AE 0.98, 7 at the same source-detector separation of 4.3 cm. The increase of DPF from infant to adult may be attributed to either a change of the optical properties of the brain, or the change of the size of the head, or a combination of both. For the appreciation of the effect of the changes of optical properties on DPF, we evaluated the DPFs of a semi-infinite geometry and a spherical geometry of 50 mm radius, using the set of Eqs. (23), (24), and (26), at three different combinations of optical properties. The first case Fig. 3 FEM meshes for the configurations to be compared. In the semi-infinite geometry, one source is set in the middle-sagittal plane and at 3 cm away from one edge, and a total of 50 detector points at an interval of 1 mm are placed in the middle-sagittal plane to cover a source-detector range from 1 mm to 50 mm. In the sphere geometry and the cylinder geometry evaluating along the azimuthal direction, one source is at the middle azimuthal plane, and a total of 50 detector points are placed in the middle azimuthal plane to have a range of the source-detector distance ("line-of-sight" or the chord distance) from 1 mm to 50 mm with an interval of 1 mm. In these two configurations, the detector with 50 mm distance from the source is set at the direct opposite of the source point. is to fix the absorption coefficient at μ a ¼ 0.01 mm −1 , with the μ 0 s changing from 0.5 mm −1 to 1.5 mm −1 . The results are shown in Fig. 6, where it is observed that the three folds of increase of μ 0 s at the specific fixed value of μ a nearly doubled the DPFs for both the semi-infinite geometry (the upper or blue profile) and the spherical geometry (the lower or green profile). It is also observed from Fig. 6 that the absolute difference of the DPFs between the semi-infinite and spherical geometry seems to be insensitive to the μ 0 s . The second case is to fix the reduced scattering coefficient at μ 0 s ¼ 1.0 mm −1 , with the μ a changing from 0.01 mm −1 to 0.1 mm −1 . The results are shown in Fig. 7, where it is observed that the 10 folds of increase of μ a at the given specific fixed μ 0 s has a nearly 50% reduction of the DPFs for both the semi-infinite geometry (the upper or blue profile) and the spherical geometry (the lower or green profile). It is also observed from Fig. 7 that the absolute difference of the DPFs between the semi-infinite and spherical geometry diminishes as the μ a increases. The third case is to fix the ratio of the absorption coefficient to the reduced scattering coefficient to 1∕100, with the μ 0 s changing from 0.5 mm −1 to 1.5 mm −1 . The results are shown in Fig. 8, where it is observed that if the absorption-to-reduced scattering ratio is kept the same, the three folds of increase of μ 0 s has increased the DPFs approximately 20% for both the semi-infinite geometry (the upper or blue profile) and the spherical geometry. In comparison, the same change of μ 0 s at a fixed μ a resulted in nearly 100% changes of DPF as shown in Fig. 6.
The age dependence of the DPFs reported by a number of studies was summarized in Ref. 8. It was shown that DPF increased an amount of more than 1 from neonatal to 25 years of age, and an amount of slightly less than 1 from 25 years of age to 50 years of age. 8 Ijichi et al. 41 41 on the no-or-little change of DPF over the range of postconceptional age of the subjects despite the significant changes of the scattering and absorption properties of the neonatal brain. If the near doubling of the reduced scattering coefficient were nearly counterbalanced by the near doubling of the absorption coefficient as a result of the increase of the cerebral hemoglobin content, the absorption-to-reduced scattering ratio could have remained relatively consistent for the subject. If that were the case, the variation of the DPFs of the subjects over the range of reduced scattering change could be small and not easily distinguishable from the inherent intersubject variations.
This work intends to evaluate how the shape as well as size of the diffusing medium (bounded in a regular geometry like a sphere) will affect the DPF. This initial line of evaluation is based on treating the medium as homogeneous. If one considers reaching the curved geometry from a semi-infinite geometry that extends to infinities at two dimensions of the direction, bending on one dimension of the direction will make the effect of an infinite-length cylinder, and bending on both dimensions of Fig. 5 The effect of the radius (in mm) of the spherical geometry on the DPF evaluated over 3 mm to 50 mm range, when compared to semi-infinite geometry for the same source-detector distances, are shown. The horizontal red dashed line marks the 5% error with respect to the semi-infinite case. Fig. 6 The dependence of DPF on the absorption coefficient from 0.01 to 0.1 mm −1 at a fixed reduced scattering coefficient of 1.0 mm −1 . The upper blue profile corresponds to the DPF of a semi-infinite geometry, and the lower green profile a spherical geometry of 50 mm in radius, for source-detector distances from 1 mm to 50 mm in both geometries. Fig. 7 The dependence of DPF on the reduced scattering coefficient from 0.5 to 1.5 mm −1 at a fixed absorption-to-reduced scattering ratio of 1∕100. The upper blue profile corresponds to the DPF of a semiinfinite geometry, and the lower green profile a spherical geometry of 50 mm in radius, for source-detector distances from 1 mm to 50 mm in both geometries. the direction will reach the effect of a spherical domain. From such an intuitive perspective, we could expect that the effect to photon diffusion due to a spherical curvature (as is the result of double-dimensional bending of the semi-infinite domain) likely will double that due to the azimuthal direction of an infinitelength cylinder geometry (as is the result of 1-D bending of the semi-infinite domain). The expressions of the photon fluence rate as given in the set of Eqs. (23), (24), and (26), which are resulted from the approximations of their respective general solutions, do show a quasidoubling effect of the spherical geometry when compared to the azimuthal direction of the infinite-length cylinder geometry. Such doubling effect of the spherical geometry compared to the azimuthal direction of the cylinder geometry with respect to the semi-infinite geometry is observable in Fig. 4, wherein the trend shown by the FEM becomes more pronounced at the radius of 50 mm compared to the radius of 25 mm. It is also evident that, as the radial dimension of either the spherical or cylindrical geometry increases, the deviation of the DPFs of the two curved geometries from that of the semiinfinite geometry diminishes. Therefore, in the cases of adult head that easily exceeds a local radius of 100 mm, the convention of using the DPFs estimated for a semi-infinite geometry is appropriate for source-detector distance greater than 2.5 cm; whereas for spherical geometries with sizes comparable to infant head, adoption of the DPFs corresponding to semi-infinite geometry should be used with caution.
With the development of denser optodes in fNIRS as well as the attention to neonatal applications, however, the noticeable changes of the DPF at smaller source-detector distances and at a spherical geometry comparing to the semi-infinite geometry may become too significant not to be accounted for when inferring clinical markers for intervention. As shown in Fig. 8, for example, the DPF corresponding to a 2 cm source-detector separation in a sphere of 50 mm radius as is comparable to the size of a neonatal head, can be approximately 5% smaller than that associated with a semi-infinite geometry at the optical properties of the medium set at μ a ¼ 0.01 mm −1 , μ 0 s ¼ 1.0 mm −1 . A 5% reduction of the DPF translates to 5% more changes of absorption coefficient at the same level of fNIRS signal change. If the shape dependency of the DPF is not accounted for, there will be increased probability of underestimating the absorption changes by assuming the DPFs of semi-infinite medium, and an amount of 5% underestimation of the change of absorption coefficient may be a concern if the baseline parameter to be monitored is close to a clinically important threshold. We note that the DPFs reported by Ijichi et al. 41 were experimentally determined by using the phase of the intensity modulated light, where the μ a and μ 0 s were estimated by fitting the absolute DPF to a semi-infinite model. Therefore, when the sets of μ a and μ 0 s were test plotted into a figure similar to Figs. 6-8, unsurprisingly, they always appeared to agree with the semi-infinite model. Were the phase-resolved measurements fitted using a model that had accounted for the smaller head size and shape, the resulted absolute changes of the μ a and μ 0 s of the neonatal brains could have been different.
This study is limited as only the simplest homogeneous domain has been considered. Irrespective of the geometry dependence, the DPF is shown to decrease at short sourcedetector separation down to about 5 mm. This pattern of DPF associated with homogeneous domain is consistent with previous reports developed using MC calculations 42 and the analytical solution of the diffusion equation incorporating the zero-boundary condition. 18 However, there are experimental studies of the adult head 43 suggesting that the DPF increases at short source-detector separation (∼2 cm). This trend of DPF increasing at short source-detector separation obtained from actual tissue was justified as relating to the inhomogeneity of tissue domain with which the use of partial differential pathlength factor (PDPF) 44 was articulated. An immediate improvement of the modeling of the tissue structure to accommodate PDPF is a multilayer geometry accounting for the effect of the blood content of cerebrospinal fluid on the optical properties of brain as indicated by Robertson et al. 45 It, nevertheless, could be expected that, the qualitative pattern of the geometry dependence of PDPF, specifically how it changes over different geometries at fixed optical properties and source-detector separation, may remain indifferent from that of DPF. Future works to extend the presented analytical approach to multilayer geometry are under planning. Furthering the analytical understanding of the geometry effect to the DPF and PDPF may eventually lead to quantitative, rapid, and accurate estimation of the μ a and μ 0 s that would be particularly influential in the clinical evaluation and management of disorders associated with early brain development and acute-care cerebral monitoring.

Conclusion
This work analytically examined the DPF of steady-state photon diffusion in a homogeneous medium associated with a few regularly shaped geometries, including a half-plane, an infinitelength cylinder, and a sphere. With the analytical procedures established for steady-state photon diffusion in an infinite-length cylinder domain that can be extended to deriving DPF, the solution to steady-state photon diffusion in a spherical domain is developed by applying the "image" source method with the EBC and involving the modified spherical Bessel functions of the first and the second kinds. The solution is converted into a format employing the physical source and its image with respect to its associated semi-infinite geometry and a radiusdependent term accounting for the dimension of the sphere. The solutions can be simplified at the general condition of μ eff R 0 ≥ 10, that is, the radius (diameter) is more than 10 (20) times of the effective attenuation pathlength, 1∕μ eff , a form that involve regular functions without lengthy summations. With these analytical preparations, the DPFs of steady-state photon diffusion in homogeneous medium associated with the three bounded geometries, including the semi-infinite geometry, along the azimuthal direction of the infinite-length cylinder geometry, and the spherical geometry, are compared for the same line-of-sight source-detector distance. The patterns of DPF estimated for the two curved geometries with a radius of 25 mm or 50 mm agree in the global patterns with those given by FEM. Quantitatively, the DPF approximated for a source-detector separation range of 10 to 25 mm on a spherical medium of 50 mm radius with μ a ¼ 0.01 mm −1 and μ 0 s ¼ 1.0 mm −1 is on average 5% different compared to the value obtained from FEM. ðR 0 ; θ; φÞ and ðR 0 þ d; θ; φÞ, are presented in Fig. 9. The results are presented in terms of lnðd Ã ΨÞ versus d, as it leads to a convenient linear relationship as evident from Eq. (2). At a radial coordinate of R 0 ¼ 50 mm, it takes the summation of 30 terms with Eq. (41) to reach an error less than 0.1% at source-detector distance greater than 10 mm, when compared to the value given by Eq. (2). At a larger radial coordinate of R 0 ¼ 200 mm, it takes summation of 100 terms with Eq. (41) to reach an error less than 0.1% at source-detector distance greater than 10 mm, when compared to the value given by Eq. (2). In the cases when the source-detector pair is not aligned with the origin (i.e., the source and the detector are separated in either azimuthal or elevational directions), the values given by Eq. (41) as a function of the source-detector distance oscillate with respect to the expected correct values (due to the limit of the numerical arithmetic in the computational software) with respect to the accurate value at the summation up to 200 terms. We note that should Eq. (41) be numerically nonidentical to Eq. (2), the numerical evaluation based on Eq. (41) for any sourcedetector configuration would have systematically deviated from that based on Eq. (2). With the comparison shown in Fig. 9, an upper limit of l ¼ 250 is implemented in the numerical evaluations of the solutions involving Eq. (3) or Eq. (41) when needed in the rest of the study. Fig. 9 Numerical comparison of the values given by Eqs. (41) and (2), at a common set of optical properties of μ a ¼ 0.01 mm −1 and μ 0 s ¼ 1.0 mm −1 , and different upper limit of l for the summation. The source and the detector are respectively at ðR 0 ; θ; φÞ and ðR 0 þ d; θ; φÞ. The results are presented in terms of the linear relationship expected for lnðd Ã ΨÞ versus d . At the radial coordinates of R 0 ¼ 50 mm and R 0 ¼ 200 mm, it takes the summation of 30 terms and 100 terms, respectively, with Eq. (41) to reach an error less than 0.1% at source-detector distance greater than 10 mm when compared to the value given by Eq. (2).