## 2 Theory
With simple mathematical manipulations, the discretized form of Equation (2) can be readily established after integrating the probability over the thickness, twice for each layer:  
\[P(m,n,\theta_{i})=\frac{\left[\exp(-\frac{C_{\text{OD}}(z_{m})}{\cos\theta_{i }})-\exp(-\frac{C_{\text{OD}}(z_{\mu})}{\cos\theta_{i}})\right]\cdot\left[- \exp(\frac{C_{\text{OD}}(z_{m})}{\cos\theta_{i}})+\exp(\frac{C_{\text{OD}}( z_{\mu l})}{\cos\theta_{i}})\right]}{C_{\text{OD}}(z_{mu})-C_{\text{OD}}(z_{ml})} \tag{5}\]  
Furthermore, Equation (4) requires numerical integration instead of analytical solutions. To the authors' knowledge, a general analytical solution for any given phase function of Equation (4) is not readily available for analytical investigation purposes. Therefore, in the Markov chain calculations, \(\varphi\) is discretized (i.e., into 900 equal intervals at a resolution of 0.2\({}^{\circ}\)) to numerically solve for \(P\left(\theta_{i},\theta_{j},n\right)\). In this study, various discretization schemes are utilized to analyze the effect of discretization. The schematic shown in Figure 1 depicts the comparison of uniform discretization (as in [21; 22]) and non-uniform discretization which are proposed by this work. The upper half of Panel (a) represents the uniform discretization and the lower half shows the non-uniform discretization, and in both schemes, 180 intervals were created. For uniform discretizations, every internal is equally divided by a one-degree resolution. For non-uniform discretizations, different discretization schemes are tested, for example, finer meshes/intervals adopted for forward directions (\(\alpha\)-0\({}^{\circ}\)) and coarse meshes/intervals for backward directions (\(\alpha\)-180\({}^{\circ}\)). Panel (a) of Figure 1 is a polar plot demonstration of a Mie scattering phase function with a resolution of 1 degree. In this work, this specific phase function shown in Panel (a) of Figure 1 is defined as phase function 1 in this manuscript. For phase function 1, a uniform particle diameter of 1.2 \(\upmu\)m, incident light wavelength of 532 nm, and ambient refractive index of 1.0 were selected. The refractive index of the particle was set to be 1.33 and unpolarized light was assumed. Different thickness discretization schemes are shown in Panel (b) of Figure 1. The left-hand side of Panel (b) depicts a uniform discretization in which each layer has the same thickness \(\Delta z\). Correspondingly, the right-hand side of Panel (b) shows a non-uniform discretization for which the thickness \(\Delta z\) in each layer may vary. The authors believe that different emphasis should be given for different cases. For example, finer meshes near the incident plane of the turbid are desired if the reflected angular distribution is sought. However, the relevance of such discretization is subject to debate for transmitted photon angular distributions, which the authors will also investigate in this work. Two additional phase functions were considered in this work and the three phase functions are depicted in Figure 2 in Cartesian coordinates. Phase function 2 is defined based on a particle refractive index of 1.57, and everything else is the same as phase function 1. To create a more structurally complex phase function, phase function 3 is defined based on a uniform particle diameter of 10.0 \(\upmu\)m, incident light wavelength of 650 nm, a particle refractive index of 1.57, and an ambient refractive index of 1.0. The phase functions are depicted in Figure 2 and the anisotropy of phase function 1, 2, and 3 are determined to be 0.8314, 0.6653, and 0.8691, respectively. For Markov chain approximations, the transition matrix calculations and the matrix inversions were carried out with MATLAB because of its excellent performance in handling matrix manipulations. An Intel Xeon X5650 CPU (12 threads, 2.66 GHz) was unitized for the MATLAB calculations.