## 2 Theoretical Models of Light Propagation and Light Scattering
### Numerical Schemes and Conditions for the RTE  
We numerically calculated the RTE based on the finite difference method. We employed the third order upwind scheme for spatial discretization and discrete ordinate method (DOM) for angular discretization [30; 31; 32]. As discrete angular directions of light intensity, we used the level symmetric even (LSE) quadrature set [33] with the total number of the directions, \(N_{\Omega}=48\). Several researchers have reported the superiority of the LSE quadrature set [34; 35]. The discrete phase function is renormalized by [PERSON]'s method [36] for accurate calculations of the scattering integral (first term in the right-hand side of Equation (1)). [PERSON]'s method works well in the current condition where the anisotropy factor \(g\) is small. For details of the upwind scheme and DOM, see references [30; 31; 32], respectively. Although the references [31; 32] focus on a case of highly forward peaked scattering rather than a case of weakly forward peaked scattering, they summarize the DOM with [PERSON]'s method. We numerically solved the discrete forms of the RTE iteratively by Bi-CGSTAB method [37] based on Eigen 3 [38]. We set the tolerance error to \(10^{-16}\). We numerically treated the RMBC based on [PERSON]'s procedure [29]; we operated the matrix for the RMBC to the results for the light intensity vector after iterations. We preliminarily confirmed that results when operating the RMBC matrix at each iteration were almost the same as the results in the current procedure. However, this procedure requires more computational time. The spatial step size was uniformly given as \(\Delta r=0.025\) cm. Our preliminary study showed that the RTE results with a smaller \(\Delta r\)-value were almost the same as those with \(\Delta r=0.025\) cm. Our in-house source code was written in the C++ programming language, and all the matrices were compressed to vectors in the compressed row storage format. We  
Figure 1: (**a**) Scattering medium, which models the alumina suspension for the RTT. The optical properties are defined in an infinitesimal volume of the medium. The light source (the red arrow) is incident on the origin. The blue dotted line represents positions where the fluence rates are calculated. (**b**) Outward and reflected directions (red and yellow arrows) at \(z=0\) surface (gray color) based on the LSE quadrature set. Black arrows denote outward, and corresponding reflected directions.  
implemented parallel CPU programming with 128 threads work station (AMD EPYC 7452) by OpenMP. The computational time at each volume fraction was approximately 1 h.  
We considered a homogeneous cubic medium with a length of 2.5 cm and an isotropic point source at the origin. In the RTE calculations, we employed the DST results for the scattering coefficient \(\mu_{s}\) and normalized phase function \(p(\theta)\) at each volume fraction. For the refractive index of the alumina suspension, we used the following equation: \(n_{cs}(\eta)=\eta n_{p}+(1-\eta)n_{w}\)[39], where \(\eta\) is the volume fraction; and \(n_{p}\) and \(n_{w}\) represent the refractive indices for alumina particles and water, respectively. We referred values of \(n_{p}\) and \(n_{w}\) as 1.768 and 1.332, respectively [40], with the optical wavelength of 600 nm. As shown in Figure 2a, the \(n_{cs}\)-value increases linearly as the volume fraction increases. According to the change of \(n_{cs}\), the reflectivity coefficient for the RMBC changes with the scattering angle as shown in Figure 2b. For the absorption coefficient, we set it to 0.1 cm\({}^{-1}\) because we focus on the scattering process of light propagation.