## 2 Theory  
The Markov chain model for multiple scattering in mediums has been thoroughly introduced in [PERSON] et al. [21; 22], assuming that the medium consists of parallel layers (the optical properties within a specific layer is uniform) and the light propagation follows the Beer-Lambert law, where the extinction coefficient consists of absorption and scattering. For simplicity, only the Mie scattering was considered in this research. The fidelity of the Markov Chain scheme in dealing with inhomogeneous medium, absorption, and anisotropic phase functions have been demonstrated. Therefore, the authors will only briefly introduce Markov chain theory and its application to photon multiple scattering. Basically, the Markov chain model for anisotropic scattering defines a transition state by the location \(z\) in the 1D turbid slab, where a scattering event takes place, and the direction/angle \(\theta\) (zenith angle) in which this photon will propagate after the scattering event. Essentially, the multiple scattering problem can be considered as a continuous Markov chain problem. To solve the Markov process, the scattering medium is discretized to convert this problem from a continuous challenge into a finite Markov chain problem, so that matrix forms can be implemented. With the discretization approach, a transition matrix \(T\) can be defined with each entry \(T((z_{m},\,\theta_{i}),\,(z_{n},\,\theta_{j}))\), representing the probability of the photon transit from one state \((z_{m},\,\theta_{i})\) to another \((z_{n},\,\theta_{j})\). Absorbing states are incorporated to represent the probability of photons exiting the slab or being absorbed by the medium. More details on the Markov formation can be found in [PERSON] et al. [22]. From photon scattering physics, the following relationships can be inferred for the multiple scattering problem, as in [PERSON] et al. [22]:  
\[T((z_{m},\theta_{i}),(z_{m},\theta_{j}))=P(z_{m},z_{m},\theta_{i})\cdot P(\theta_ {i},\theta_{j},n) \tag{1}\]  
\[P(z_{m},z_{n},\theta_{i})=\exp(-\frac{\int_{z_{m}}^{z_{m}}\Sigma_{e,z}dz}{\cos \theta_{i}}) \tag{2}\]  
\[\alpha=\arccos(\cos\theta_{i}\cos\theta_{j}\cos\varphi+\sin\theta_{i}\sin \theta_{j}) \tag{3}\]  
\[P(\theta_{i},\theta_{j},n)=\int_{0}^{2\pi}\Gamma_{n}[\arccos\alpha]d\varphi \tag{4}\]  
Whereas, in Equation (1), \(P(z_{m},z_{n},\theta_{i})\) is the probability of a photon with two consequent scattering events taking place at \(z_{m}\), \(z_{n}\) and propagates in the direction of \(\theta_{i}\), \(z_{m}\) represents the location of the prior scattering event and \(z_{n}\) represents the location of the latter scattering event. \(P(\theta_{i},\theta_{j},n)\) is the probability that a photon propagates in the direction of \(\theta_{i}\), scattered with a certain phase function in the \(n\)th layer, then propagates in the new direction of \(\theta_{j}\). In Equation (2), \(\Sigma_{e,z}\) is the extinction coefficient at the location \(z\). In Equation (3), \(\varphi\) is the azimuth angle and \(\alpha\) is the included angle between the incident direction and the exiting direction. In Equation (4), \(\Gamma_{n}(\alpha)\) is the phase function (probability density function, PDF) in the \(n\)th layer. It is worth noting that Equation (2) represents the propagating probability from point \(z_{m}\) to point \(z_{n}\). To solve the probability of the photon propagating from layer \(m\) to layer \(n\), proper discretizations and integrations are necessary. To derive such a probability, the cumulative probability density function (CPDF) of the OD can be defined as \(C_{\text{OD}}(z)\). The upper lower and upper bound of layer \(m\) and \(n\) are defined to be \(z_{ml}\), \(z_{m\mu}\), and \(z_{nl}\), \(z_{mu}\), respectively.