## 3 Computational Results  
### Uniform Discretizations  
In this section, Markov chain solutions with uniform discretizations were solved and compared with Monte Carlo simulations. It is worth noting that in the zenith angle discretizations (to be discussed in Section 3.3), both the incident angle \(\theta_{1}\) and new propagation angle \(\theta_{2}\) need to be discretized and the same discretization scheme must be applied to both angles. Such a requirement is posed by the nature of the Markov chain approximation, in that the transition states cannot be re-defined during the Markov process. Similarly, the resulting transmitted and reflected angular distributions tend to follow the same discretization that one sets at the beginning. Therefore, direct comparison among different discretization schemes will not offer sizably meaningful insights. Instead, the Markov calculations will be compared against one common benchmark with the results from the Monte Carlo simulations. Transmitted probabilities were normalized with a unit of \(1^{\circ-1}\) in this work to make sure that the summation of the probability over the range from \(0^{\circ}\) to \(180^{\circ}\) was unity minus the percentage of ballistic photons (transmitted photons without scattering, not discussed in this work). The same method was applied for all results presented in this work. Figure 3 depicts Markov chain predictions of the transmitted photon angular distributions compared with those from the Monte Carlo simulations. In the simulations, uniform OD distribution was assumed, and the total OD was set to either 2 or 5. On the \(x\) axis, the transmitted angle from \(0^{\circ}\) to \(90^{\circ}\) indicates that the photons leave the turbid slab from the bottom exit plane, whilst the transmitted angle from \(90^{\circ}\) to \(180^{\circ}\) represents the photons that are reflected and leave the slab from the entrance plane. Since the range of the OD in the simulations was within the intermediate regime, the transmitted angular distribution exhibited both features from the phase function and multiple scattering averaging effect, as reflected by the complex angular distribution structures, as reflected in Figure 3. From the results, it is evident that [PERSON] chain approximations agree well with Monte Carlo simulation results with satisfactory fidelity, for all the six cases considered in this work. Figure 3: Markov chain angular distribution predictions compared with Monte Carlo. Panel (**a**). Phase function 1, OD = 2; Panel (**b**). Phase function 1, OD = 5; Panel (**c**). Phase function 2, OD = 2; Panel (**d**). Phase function 2, OD = 5; Panel (**e**). Phase function 3, OD = 2; Panel (**f**). Phase function 3, OD = 5. Figure 4 shows the absolute relative error of the Markov chain approximations compared to the Monte Carlo simulations for the results shown in Figure 3, to further study the performance of Markov chain approximations. Note that other algorithms, such as the random walk [7], adding doubling [34], and the imperialist competitive algorithm [35] are also frequently used to evaluate such computational tasks. Since we have incorporated such comparisons in our previous work [22], we will not elaborate on such comparisons in this work. The absolute relative error in this paper is defined as:  
\[\text{error}=\left|\frac{\text{Markov }-\text{ Monte Carlo}}{\text{Monte Carlo}}\right|\times 100\% \tag{6}\]  
The absolute relative error as a function of the transmitted angle is plotted in logarithm scale in Figure 4. We used Red dash lines to indicate an error threshold of 1.0%, and blue dash lines for a threshold of 0.1%. It is worth noting that we did not use the arbitrarily chosen thresholds to define any regime, but rather to evaluate the performance of each discretization scheme as a reference. From the data, several observations could be made. Firstly, the absolute relative error for all the cases was generally quite small (mostly below 0.1% for phase function 1 and phase function 2, below 1% for phase function 3).