## 1 Introduction  
Understanding light transportation in skin tissues can help improve clinical efficacy in the laser treatment of dermatosis, such as port-wine stains (PWS) [1]. As an effective and easy-to-parallelize technique [2], the Monte Carlo (MC) method has been successfully used in modeling light transportation in biological tissues. In this method, a \"random walk\" with a sufficient number of photon packets is traced to build reliable trends of light transportation on the basis of their statistical scattering and absorption. Non-uniform energy distribution via selective photothermolysis [3] is caused by various optical properties in different skin tissues, leading to preferential absorption in blood vessels instead of in other skin tissues, including epidermis, dermis, and hair follicles [4]. Moreover, patient-specific PWS vessels are typically structurally complicated in terms of computation. Multilayered PWS vessels, particularly cross-bridge vessels, considerably influence laser propagation due to shading effects [5, 6] (Figure 1). Therefore, an unknown deposition of light energy always appears in PWS vessels [7]. The effect of the adaptability of complicated vessel structures has elicited considerable attention from researchers in MC simulation. [PERSON] et al. [9] developed the MC method in multilayered tissues (MCML) and applied it to a succession of homogeneous layers with different optical properties that are parallel to one another. MCML has been extended to adapt to the cylindrical surface of blood vessels by using Snell's law in the curved interface to trace the transport of light. This algorithm is called the geometric Monte Carlo (GMC) [10; 11] because all tissue structures and their interfaces require geometric descriptions. However, GMC is limited in the burdensome mathematical definitions of complicated tissue structures. In terms of light propagation in arbitrarily complex 3D heterogeneous tissue morphology, the voxel-based Monte Carlo (VMC) algorithm has attracted considerable attention among scholars [12; 13; 14; 15]. In VMC, the computational domain is divided into voxels by using a cuboid mesh with a simple connection relationship. Photon radiation is investigated by assigning certain optical properties to each voxel [13]. However, the round or oblique boundaries between neighboring tissues are approximated via serrated surfaces due to the fixed mesh connection. Consequently, large errors (80-120%) of photon intensity occur in boundary voxels [7]. Moreover, mesh refinement does not change the shape and arrangement of grids, and thus is unable to improve the description of light reflection and refraction at curved boundaries. Therefore, errors exist even with fine discretization [16]. Compared with voxel-based medium discretization, a tetrahedron-based approach is ideally adaptable to modeling targets with curved boundaries or locally refined structures. [PERSON] and [PERSON] [17] proposed a tetrahedron-based inhomogeneous MC optical simulator (TIM-OS) using tetrahedral meshes for MC simulations. Although TIM-OS was validated by MCML, its computational domains are mostly regular regions with less complexity. [PERSON] et al. [18] developed a tetrahedron-based MC (TMC) solver that is suitable for computing photon energy deposition in complicated domains with a body-fitted tetrahedron mesh. The error caused by confused optical parameters between neighboring cells when photons are incident along the cell boundary can be avoided in TMC by introducing the distance threshold [18; 19]. The influence of discretization on light propagation was discussed to compare the relative errors of VMC and TMC with the benchmark of GMC. A comparative study showed that energy deposition error via TMC in the interfacial region is one-tenth to one-fourth that via VMC; thus, the former yields a more accurate computation of photon transportation [18]. In TMC, the escape boundary condition (BC) is applied to computational surfaces.