1. Introduction

Diffuse optical tomography (DOT) [1–6] is a biomedical imaging technique that generates a tomographic image of a living body by estimating the distribution of absorption and scattering coefficients of biological tissues in response to near-infrared light using inverse analysis. The behavior of light is dominated by the radiative transfer equation that is an integro-differential equation. Since the radiative transfer equation explicitly handles the angular distribution of photon flux, its numerical solution in three-dimensional geometry requires a higher computational cost. In order to directly solve the radiative transfer equation, the Monte Carlo method [7], the spherical harmonics method [8,9], and the discrete ordinate method [10] can be used. However, these methods require considerable computation time. In DOT, many forward calculations should be carried out to inversely estimate tissue (or material) distribution inside a geometry. Therefore, a diffusion calculation is sometimes used to approximate the angular distribution of photon flux by the 0^th and 1^st Legendre polynomials to reduce the computational cost [11,12].

When considering a system composed of multiple materials with different refractive indices, an internal boundary condition that accounts for reflections and refractions at the boundary is desirable. However, the explicit internal boundary condition has not been commonly used in optical diffusion calculations since the refractive indexes of common biological tissue take similar values (typically in the range of 1.3-1.6) [13]. However, when DOT is applied to a human neck to detect thyroid cancers, the respiratory tract should be taken into account [14–16]. Since the refractive index of air within the respiratory tract is 1.0, the effect of internal boundary condition becomes larger. Therefore, explicit consideration of the internal boundary condition especially between the biological tissue and air in the respiratory tract is desirable.

In the transport calculation that explicitly considers the flight direction of photon flux, the internal boundary condition between two materials having different refractive indexes can be directly treated considering Fresnel's law [14]. However, the consideration of refraction in diffusion theory requires some considerations and there are a small number of references on the study of internal boundary condition for the diffusion theory. Haskel's method [17] and Aronson's method [18] have been proposed as internal boundary conditions for optical diffusion calculations considering reflection and refraction. Haskel's method uses approximations in the derivation of internal boundary condition, while Aronson's method does not, and thus the results obtained by these methods are expected to be different. In this paper, an alternative method is proposed to improve the accuracy by iteratively updating the effective reflection coefficient proposed by Haskel using the scalar fluxes and net currents obtained during the diffusion calculation. The major purpose of the present study is proof of the principle of the present method that estimates an appropriate effective reflection coefficient using the iterative procedure. Therefore, numerical verification calculations are carried out in a simple one-dimensional geometry in this study.

The novelty of this paper is twofold. First, a new method for iteratively updating the effective reflection coefficients is proposed to treat the explicit internal boundary condition for the diffusion theory. Second, systematically investigate the accuracy of Haskel's method, Aronson's method, and the proposed method by comparing them with the reference solutions obtained by a Monte Carlo method. The numerical verification calculations consider the application of DOT to a human neck that includes a respiratory tract in which air with a very different refractive index exists.

The rest of this paper is organized as follows. In Section 2, Haskel's method, the proposed method, and Aronson's method are described, and the computational accuracy of each method is qualitatively discussed. Section 3 presents the numerical results of a benchmark problem for a one-dimensional slab geometry using Haskel's method, the proposed method, and Aronson's method. The computational accuracy of each method is discussed by comparison with the spatial distribution of the photon scalar fluxes obtained by the Monte Carlo method. Section 4 presents the conclusions of this paper.

2. Theories for internal boundary condition

In this section, the theory for Haskel's method, the proposed method, and Aronson’s method are explained. In these methods, the internal boundary condition with different refractive index materials is considered by the effective reflection coefficient for partial photon flux obtained by the diffusion theory. The difference between these methods comes from the assumptions to derive the effective reflection coefficient.

Haskel's method accounts for reflection-refraction phenomena by calculating the reflected partial currents using the effective reflection coefficient, which is the ratio of incoming partial currents to outgoing partial currents at the material boundary. However, it does not take into account incoming partial currents from the adjacent mesh in the derivation of the effective reflection coefficient. Therefore, the derived effective reflection coefficient is appropriate for external boundary surfaces but is expected to be inaccurate for internal boundary.

In the proposed method, the effective reflection coefficient is calculated from the scalar fluxes and net currents through iterative calculations. The proposed method takes into account incoming partial currents from the adjacent mesh, so the accuracy at the internal boundary is expected to be improved over Haskel's method.

Aronson's method assumes a hypothetical boundary having a small thickness and then calculates the net currents inside the boundary from the difference of the partial currents transmitted to the internal boundary. The accuracy at the internal boundary is expected to be better than Haskel's method because it can account for incoming partial currents from the adjacent mesh in the calculation of the net photon currents at the internal boundary.

In the following subsections, some details on these three methods are provided.

2.1 Haskel’s method

In Haskel's method [17], the reflection of partial currents $J_{}^{refle}$ is calculated by Eqs. (1) and (2) using the effective reflection coefficient.

It is noted that notations are summarized in Fig. 1. As shown in Fig. 1, region i is on the left, and region $i + 1$ is on the right, with the direction to the right being positive. The subscripts L and R indicate that they are physical quantities related to the left and right sides of the internal boundary. The effective reflection coefficient is the ratio of incoming partial currents to outgoing partial currents at the internal boundary and is defined by Eqs. (3) and (4).

 

Fig. 1. Physical quantities at an internal boundary

Download Full Size | PDF

Based on the following notations, the derivation of the effective reflection coefficient is described below.

$J_{i + 1/2,L}^n = \; $ net currents on the left side of the boundary $i + 1/2$

$J_{i + 1/2,R}^n = \; $ net currents on the right side of the boundary $i + 1/2$

$J_{i + 1/2,L}^{out} = \; $ outgoing partial currents from the region i to the left of the boundary $i + 1/2$

$J_{i + 1/2,R}^{out} = \; $ outgoing partial currents from the region $i + 1$ to the right of the boundary $i + 1/2$

$J_{i + 1/2,L}^{in}\; = \; $ incoming partial currents from the left side of the boundary $i + 1/2$ into the region $i$

$J_{i + 1/2,R}^{in} = \; $ incoming partial currents from the right side of the boundary $i + 1/2$ into the region $i + 1$

${\phi _{i + 1/2,L}} = \; $ Scalar fluxes on the left side of the boundary $i + 1/2$

${\phi _{i + 1/2,R}} = $ Scalar fluxes on the right side of the boundary $i + 1/2$

${R_{eff,i + 1/2,L}} = $ Effective reflection coefficient of the boundary $i + 1/2$ from the left side

${R_{eff,i + 1/2,R}} = $ Effective reflection coefficient of the boundary $i + 1/2$ from the right side

${T_{eff,i + 1/2,L}}$ = Effective transmission coefficient of the boundary $i + 1/2$ from the left side and ${T_{eff,i + 1/2,L}} = {1–R_{eff,i + 1/2,L}}$

${T_{eff,i + 1/2,R}} = $ Effective transmission coefficient of the boundary $i + 1/2$ from the right side and ${T_{eff,i + 1/2,R}} = {1–R_{eff,i + 1/2,R}}$

$\mathrm{\Delta }{x_i} = $ width of a spatial mesh i

${D_i} = $ Diffusion coefficient of region $i$

The reflected partial current of photons entering the internal boundary from the left side is given by Eq. (5) and that entering from the right side is given by Eq. (6).

It is noted that $\mu = 1$ in the x-axis positive ($+ $) direction and ${\mu _ - } = 1$ in the x-axis negative ($- $) direction. There is a relationship of $\mu ={-} {\mu _ - }$. The incident angle-dependent reflection coefficient is expressed by Eq. (7) from Fresnel's law [17].

The angle of refraction and the critical angle are calculated by Eqs. (8) and (9) [17].

In the diffusion calculation for a one-dimensional slab geometry, the angular distribution of the flux is treated with P0 and P1 components of the Legendre polynomials as shown in Eq. (10).

From Eqs. (12) and (13), the reflectance for the 0^th-order isotropic components ${\phi _{0,L}}$ and ${\phi _{0,R}}$ can be calculated by $\mathop \int \nolimits_0^1 R(\mu )\mu d\mu $ and $\mathop \int \nolimits_0^1 R({{\mu_ - }} ){\mu _ - }d{\mu _ - }$, respectively. Similarly, for the 1^st -order anisotropic components ${\phi _{1,L}}$ and ${\phi _{1,R}}$, reflectance is given by $\mathop \int \nolimits_0^1 R(\mu ){\mu ^2}d\mu $ and $\mathop \int \nolimits_0^1 R({{\mu_ - }} )\mu _ - ^2d{\mu _ - }$, respectively. This relation is illustrated in Fig. 2. Finally, $J_{i + 1/2,L}^{refle}$ and $J_{i + 1/2,R}^{refle}$ can be derived as Eqs. (14) and (15).

 

Fig. 2. Conceptual diagram of reflection for 0^th and 1^st order components

Download Full Size | PDF

$J_{i + 1/2,L}^{out}$ and $J_{i + 1/2,R}^{out}$ can be expressed as Eqs. (20) and (21) using scalar fluxes (0^th component) and net currents (1^st component).

By substituting Eqs. (14) and (20) into Eq. (3), and Eqs. (15) and (21) into Eq. (4), we have:

$J_{i + 1/2,L}^{in}$ and $J_{i + 1/2,R}^{in}$ can be expressed as Eqs. (24) and (25) using scalar fluxes and net currents.

By assuming that incoming partial currents from adjacent regions are negligible (i.e., $J_{i + 1/2,L}^{refle} = J_{i + 1/2,L}^{in},\; \; J_{i + 1/2,R}^{refle} = J_{i + 1/2,R}^{in}$), the following equations are obtained from Eqs. (14), (24), and Eqs. (15), (25).

Equations (26) and (27) can be rearranged to Eqs. (28) and (29).

However, to derive Eqs. (26) and (27), incoming partial currents from adjacent regions are neglected. Therefore, the effective reflection coefficient obtained from Eqs. (30) and (31) correctly satisfies Eq. (1) at the external boundary where no incoming partial current exists. However, it does not exactly satisfy Eq. (1) at the internal boundary.

The finite-difference diffusion equation considering the effective reflection is derived as follows. Based on Fick's law and the finite-difference approximation, $J_{i + 1/2,L}^n$ and $J_{i + 1/2,R}^n$ can be calculated using Eqs. (32) and (33), respectively.

Expressing the net currents by the subtraction of the partial currents, Eqs. (34) and (35) are obtained. Here, the incoming partial currents are given by the addition of the reflected and the transmitted partial currents.

Substitution of Eqs. (20) and (21) into Eqs. (34) and (35), and utilization of the relation for the reflection and transmission probabilities (${T_{eff}} = 1 - {R_{eff}}$) yield Eqs. (36) and (37), respectively.

Using Eqs. (32) and (33), Eqs. (38) and (39) can be rearranged for scalar fluxes.

By substituting Eqs. (31) and (32) into Eqs. (29) and (30), we have Eqs. (33) and (34).

By eliminating $J_{i + 1/2,R}^n$ from Eqs. (40) and (41), utilizing the relation for the reflection and transmission probabilities (${T_{eff}} = 1 - {R_{eff}}$), and rearranging, $J_{i + 1/2,L}^n$ is given by Eq. (42).

Similarly, we have the following relation.

The right-hand sides of Eqs. (42) and (43) are identical. It confirms that the net currents at the left and right sides of the boundary are identical. It is noted that the equivalence between $J_{i + 1/2,L}^n$ and $J_{i + 1/2,R}^n$ can be derived by subtracting both sides of Eqs. (40) and (41).

Next, let us consider the external boundary at the left end of the system. The net current on the right side of the external boundary can be calculated using Eq. (44).

Using the partial current, the net current is given by Eq. (45).

$J_{1/2,R}^{out}$ can be expressed as Eq. (46) using scalar fluxes and net currents.

By substituting Eq. (39) into Eq. (38), we have Eq. (40).

By substituting Eq. (47) into Eq. (44) and by rearranging, $J_{1/2,R}^n$ is given by Eq. (48).

Finally, let us consider the external boundary at the right end of the system. Similar to the Eqs. (44), (45), and (46), we have Eqs. (49), (50), and (51).

By substituting Eq. (45) into Eq. (42) and by rearranging, $J_{{N_x} + 1/2,L}^n$ is given by Eq. (46).

The finite-difference diffusion equation considering the internal and external boundary conditions is given by Eq. (54) because of Eqs. (42), (43), (48), and (53).

Therefore, to obtain the numerical solution of the diffusion calculation, the three-point differentiation equation shown in Eq. (56) should be solved.

2.2 Proposed method

As explained in Section 2.1, the Haskel's method neglects the incoming partial currents from adjacent regions to derive the effective reflection coefficient. Thus, the accuracy of Haskel's method would be less accurate at an internal boundary. To address this issue, the proposed method the proposed method uses iterative calculations.

The effective reflection coefficients can be calculated using not only Eqs. (30) and (31), but also Eqs. (22) and (23). Using the effective reflection coefficients obtained by Eqs. (30) and (31) as initial values, the scalar fluxes and net currents are calculated by performing a diffusion calculation. Substituting the obtained scalar fluxes ($\phi _{i + 1/2,L}^{Diff},\phi _{i + 1/2,R}^{Diff}$) and net currents ($J_{i + 1/2,L}^{Diff},J_{i + 1/2,R}^{Diff}$) into Eqs. (22) and (23), the effective reflection coefficients can be updated as in Eqs. (60) and (61).

Equations (60) and (61) use net currents to calculate effective reflection coefficients, so they can account for incoming partial currents from adjacent regions. The calculation procedure of the proposed method is summarized as follows.

2.3 Aronson's method

Aronson's method considers [18] a virtual internal boundary having a very thin thickness, as shown in Fig. 3.

 

Fig. 3. Internal boundary condition of Aronson's method

Download Full Size | PDF

$J_L^{out}$ and $J_R^{out}$ are expressed by Eqs. (62) and (63), respectively.

Since the fraction of partial currents transmitting the boundary is $1 - R$, the currents $J_L^ + $ flying rightward inside the boundary and the currents $J_R^ - $ flying leftward inside the boundary are given by Eqs. (64) and (65) from Eqs. (62) and (63), respectively.

Net currents on the left side of the boundary, inside the boundary, and right side of the boundary should be identical.

Therefore, $J_L^ + $ and $J_R^ - $ can be rearranged as in Eqs. (67) and (68).

${T_{\phi ,L}}$, ${T_{J,L}}$, ${T_{\phi ,R}}$, and ${T_{J,R}}$ are given by:

In the integration of Eqs. (69)–(72), the critical angle must be taken into account. ${T_{\phi ,L}}$ and ${T_{\phi ,R}}$ correspond to the transmittance for the 0^th-order isotropic component, and ${T_{J,L}}$ and ${T_{J,R}}$ correspond to the transmittance for the 1^st-order anisotropic component. It is noted that in the Haskel method, these coefficients (${R_\phi }$, ${R_J}$) are related as ${R_\phi } + {T_\phi } = 1$ and ${R_J} + {T_J} = 1$. The net currents are given by Eq. (73).

By substituting Eqs. (59) and (60) into Eq. (65) and using the relationship of ${R_J} + {T_J} = 1$, we have Eq. (66).

${J^n}$ includes the reflection of the 0^th-order and 1^st-order components of the flux as in Haskel's method, but the approximation in Haskel's method ($J_{i + 1/2,L}^{refle} \approx J_{i + 1/2,L}^{in}$) is not used. Therefore, Eq. (66) is expected to be more accurate than Haskel's method. From Eq. (74),

From Fick's Law, we have:

By substituting Eq. (70) into Eqs. (67) and (68), we have Eqs. (71) and (72).

From Eqs. (79) and (80), the net currents are given by:

The diffusion equation in $i \ne 1,\; {N_x}$ becomes Eq. (83) from Eqs. (81) and (82).

From the above, the three-point differential equation should be solved to obtain a numerical solution to the diffusion calculation. In Eq. (83), the coefficients ${A_{i - 1/2}}$, ${A_i}$, and ${A_{i + 1/2}}$ are given by:

The coefficients shown in Eqs. (84), (85), and (86) are similar when compared to those derived from Haskel's method (Eqs. (57), (58), and (59)). However, reflection and transmission coefficients in the equation are different. Haskel's method uses effective reflection coefficients for all terms. In Aronson's method, the third term in the denominator uses the transmission coefficient for the anisotropic component, while the other terms use the transmission coefficient for the isotropic component.

2.4 Comparison of the three methods

Haskel's method and the proposed method take into account reflection and refraction phenomena by calculating the reflected partial currents at the internal boundary using the effective reflection coefficient ${R_{eff}}$. It is important to accurately calculate the value of the effective reflection coefficient because the reflected partial currents are determined by the value of the effective reflection coefficient.

Haskel's method defines the reflectance ${R_\phi }$ and ${R_J}$ for the 0^th-order isotropic and 1^st-order anisotropic components, respectively. In Haskel's method, ${R_\phi }$ and ${R_J}$ are used to calculate the effective reflection coefficient for outgoing partial currents. However, when currents are incoming from the adjacent region, the ratio of isotropic to anisotropic components is not uniquely determined, and therefore the effective reflection coefficient is also not uniquely determined. Haskel's method eliminates the above issue by approximating the incoming partial currents from adjacent regions as zero. This approximation is expected to obtain an accurate result at the external boundary, but less accurate at the internal boundary where the incoming partial current is not zero.

In contrast, the proposed method uses an iterative method to calculate the effective reflection coefficients. This method can take into account incoming partial currents from adjacent regions. By performing a diffusion calculation using the effective reflection coefficient of Haskel's method as the initial value and determining the ratio of isotropic to anisotropic components from the results, the updated effective reflection coefficient can be obtained, which is expected to be closer to the true value. The converged effective reflection coefficients will improve the accuracy at the internal boundary.

Unlike the other two methods, Aronson's method does not use the effective reflection coefficient. Considering a hypothetical boundary with a very thin thickness, the net currents inside the boundary are calculated from the difference of the transmitted partial currents into the boundary interior. The accuracy of the calculations at the internal boundary is expected to be higher because incoming partial currents from adjacent regions are taken into account in the process of calculating the net currents inside the boundary.

From the above, it is expected that the accuracy of the three methods is equivalent for the external boundary, while for the internal boundary, the proposed method and Aronson's method are expected to be more accurate than Haskel's method.

3. Results and discussions

To compare the accuracy of Haskel's method, the proposed method, and Aronson's method, a light diffusion calculation code based on these methods was developed to compute the spatial distribution of photon scalar fluxes. These results were compared to a reference solution obtained by the Monte Carlo code MCML [19] which was modified to treat a spatially uniform and isotropic photon source in a tissue region.

The calculation geometry is a one-dimensional slab geometry as shown in Fig. 4. The calculation geometry was set up assuming a simplified human neck, which is composed of human tissue and respiratory tract. Photons are generated uniformly and isotropically only in the region of material 1 $(0.0 < x < 0.1\; \textrm{cm})$ and transported through the region of material 2$(0.1 < x < 0.2\; \textrm{cm})$, which is assumed to be air, to the region of material 3 $(0.2 < x < 0.3\; \textrm{cm})$, which is assumed to be tissue. In this process, at the internal boundary at $x = 0.1\; \textrm{cm}$, the incidence is mainly occurred from a high refractive index material to a low refractive index material. At the boundary at $x = 0.2\; \textrm{cm}$, the incidence is occurred from a low refractive index material to a high refractive index material. Note that the outside of the system is assumed to be a vacuum, and once a photon leaves the system, it will not re-enter the system.

 

Fig. 4. Calculation geometry.

Download Full Size | PDF

The optical characteristic values of each material are shown in Table 1. The scattering coefficient ${\mu _s}$ and the absorption coefficient ${\mu _a}$ of the system were set based on the typical values of human tissues [13,14]. The refractive indices of material 2 and the outside of the system were set to 1.0, and the refractive indices of material 1 and material 3 were set to five conditions: 1.0, 1.2, 1.4, 1.6, and 1.8. The same values were set for material 1 and material 3. In this system, the reflection-refraction phenomenon occurs at all internal and external boundaries. When the refractive indexes given to materials 1 and 3 become larger, the effect of the reflection-refraction phenomenon is larger. The number of mesh divisions is 100 for each material (i.e. 300 meshes in total), and the convergence criterion for the effective reflection coefficient in the proposed method is ${10^{ - 15}}$. The number of photon histories used in the MCML calculation is 10¹⁰.

Table 1. Optical characteristic values of each material

View Table

The scalar fluxes calculated by Haskel's method, the proposed method, and Aronson's method, as well as the reference solution of MCML, are shown in Figs. 5 to 9. The relative error of each method to the MCML reference solution is also shown. Relative errors were calculated by Eq. (87).

 

Fig. 5. Spatial distribution of scalar fluxes and errors relative to MCML ($n = 1.0$) (Data File 1)

Download Full Size | PDF

 

Fig. 6. Spatial distribution of scalar fluxes and errors relative to MCML ($n = 1.2$). (Data File 2)

Download Full Size | PDF

 

Fig. 7. Spatial distribution of scalar fluxes and errors relative to MCML ($n = 1.4$). (Data File 3)

Download Full Size | PDF

 

Fig. 8. Spatial distribution of scalar fluxes and errors relative to MCML ($n = 1.6$). (Data File 4)

Download Full Size | PDF

 

Fig. 9. Spatial distribution of scalar fluxes and errors relative to MCML ($n = 1.8$). (Data File 5)

Download Full Size | PDF

The following trends can be observed from Fig. 5 to Fig. 9.

The reason for the last observation is given as follows. As the refractive index increases, the reflection at the internal boundary increases, and thus the transmission decreases. In other words, incoming partial current from adjacent region decreases. Thus, the results of Haskel's method become accurate because Haskel's method assumes no incoming partial current from adjacent region ($J_{i + 1/2,L}^{refle} = J_{i + 1/2,L}^{in},\; J_{i + 1/2,R}^{refle} = J_{i + 1/2,R}^{in}$).

Focusing on the accuracy of the proposed method, the error of the proposed method is approximately 4% to 7% lower than that of Haskel's method near the internal boundary, which is identical to that of Aronson's method. Note that the error at the external boundary exceeds 9%. The angular distribution of the flux becomes anisotropic due to the outflow of photons from the system near the external boundary. As a result, the diffusion approximation is no longer valid and the calculation accuracy is reduced. Reducing this error at an external boundary is beyond the scope of this paper and will not be further discussed. There are several works on the treatment of external boundary, e.g., Ref. [18].

The present and Aronson’s methods show higher accuracy when the difference of refractive indexes becomes larger. Since the typical refractive index for a biological tissue is in the range of 1.3-1.6, consideration of the internal boundary condition is important when the air (the refractive index $n = 1.0$) within the tissue is considered. A typical situation is the detection of thyroid cancers in a human neck. For such a situation, the present and Aronson’s methods will be more suitable for DOT.

Next, the convergence of the effective reflection coefficient in the proposed method is discussed. The residuals of the effective reflection coefficient at $n = 1.4$ during iteration is shown in Fig. 10. The residual of the effective reflection coefficient was calculated by Eq. (88).

 

Fig. 10. The residuals of the effective reflection coefficients during iteration

Download Full Size | PDF

From Fig. 10, it is observed that the effective reflection coefficient satisfies the convergence criterion after 9 iterations. This result suggests that the convergence of the proposed method is good.

The proposed method would be considered a candidate for practical DOT calculations since the complicated nature of reflection and transmission at the material interface can be treated with one effective reflection coefficient. In the framework of diffusion approximation, Aronson’s method will be faster than the present method since no iteration is necessary. The present method will be advantageous for accurate calculations considering higher-order anisotropy of photon flux considering transport effect. For example, let us consider a local photon transport problem considering a small region near an internal boundary. The effective reflection coefficient of the internal boundary can be calculated using Eqs. (3) and (4) instead of Eqs. (60) and (61) using the result of photon angular flux distribution obtained in the local photon transport problem. Once the effective reflection coefficient is obtained, it can be used in a diffusion calculation for a whole geometry. In this case, though the diffusion theory is used to calculate a whole geometry, the internal boundary condition (the effective reflection coefficient) captures the transport effect at the material interface. When the photon flux distribution is obtained in a whole geometry, the local photon transport problem near the internal boundary can be recalculated using the updated boundary condition (e.g., more appropriate net photon current obtained in a whole geometry calculation). By repeating this calculation procedure, most of the transport effect at the internal boundary would be captured.

On the contrary, the treatment of Aronson’s method will be more complicated since the independent reflection and transmission coefficients are necessary for higher-order moments of photon angular flux. However, the present method still can be used with the identical framework when the partial current photon fluxes are calculated by the higher-order transport method as described above. Since the objective of the present study is proof of the principle of the present method using the iterative approach, incorporation of higher-order transport calculation and quantitative comparison of computational performance with other methods will be future tasks.

4. Conclusions

In this paper, a new method for the treatment of internal boundary condition in the diffusion calculation using the effective reflection coefficient is proposed.

An existing internal boundary condition using the effective reflection coefficient is Haskel's method. Haskel's method defines the reflection coefficients ${R_\phi }$ and ${R_J}$ for the 0^th-order isotropic and 1^st-order anisotropic components of the photon angular flux, respectively, and uses ${R_\phi }$ and ${R_J}$ to calculate the effective reflection coefficient ${R_{eff}}$ for outgoing partial currents. However, when entering partial currents from adjacent regions exists, ${R_{eff}}$ is not uniquely determined because the ratio of isotropic and anisotropic components is not uniquely determined. Haskel's method eliminates the above issue by approximating the incoming partial current from the adjacent mesh is zero and then calculating ${R_{eff}}$. This approximation will be correct for the external boundary, but the accuracy for the internal boundary will decrease.

The proposed method uses a self-consistent iterative method to compute the effective reflection coefficient, which takes into account the incoming partial current from adjacent regions. By performing a diffusion calculation using the effective reflection coefficient from Haskel's method as the initial value and determining the ratio of the isotropic and anisotropic components of the photon angular flux from the results, the effective reflection coefficient can be obtained that is closer to the true value. Namely, a more accurate solution can be numerically solved by updating the effective reflection coefficients using the above procedure and using the converged effective reflection coefficients. The proposed method would have an advantage for a more complicated calculation considering higher order angular distribution of photons since only one effective reflection coefficient is considered to capture the cumbersome nature of reflection/transmission at an internal boundary with different material.

The alternative method for handling internal boundary condition that does not use an effective reflection coefficient is Aronson's method. In the Aronson method, the net currents inside the internal boundary are calculated from the difference of the partial currents transmitted inside the boundary having a hypothetical small width. Since incoming partial currents from adjacent regions can be taken into account in the net currents evaluation inside the boundary, the diffusion calculation can be performed accurately at the internal boundary.

To compare the computational accuracy of Haskel's method, the proposed method, and Aronson's method, the spatial distribution of photon scalar fluxes in a one-dimensional slab geometry was calculated using the three methods. In addition, they were compared to a reference solution obtained by the Monte Carlo code MCML. As a result, it was confirmed that the results of the proposed method and Aronson's method are consistent and the proposed method and Aronson's method are more accurate than Haskel's method near the internal boundary. The convergence of the effective reflection coefficient in the proposed method is good.