# Finite Element Method (FEM) Modeling of Laser-Tissue Interaction during Hair Removal
## Appendix B Heat Diffusion Equation  
Formalism of the heat diffusion process was described for each tissue type and each FEM element as follows [14]:  
\[\rho c_{p}\frac{\partial T}{\partial t}=\
abla(k\
abla T)+Q\]  
where \(T\) is temperature, \(\rho\) is density, \(c_{p}\) is specific heat capacity, \(k\) is thermal conductivity, and \(Q\) is the heat source due to absorption of photons, which was calculated by the mesh-based Monte Carlo method (MMCM). The initial temperature profile was in the form of a linear temperature gradient, namely, the temperature varied continuously from the surface of the skin (i.e., the top of the epidermis) to the bottom of the subcutis from 32 \({}^{\circ}\)C to 37 \({}^{\circ}\)C, after which it remained constant. The model was subject to the convective thermal boundary condition.  
\[k\frac{\partial T}{\partial z}=h(T-T_{air})\]  
where \(k\) is the conductivity of hair or epidermis, \(h\) is a constant heat transfer coefficient, and \(T_{air}\) is the ambient temperature (20 \({}^{\circ}\)C).  
Appendix C Validating a Finite-Element Method Model for Simulating Laser-Tissue Interaction in Hair Removal Systems  
The FEM model for validating our simulations with the results of [PERSON] et al. [14] had a cross section of 2 mm \(\times\) 2 mm, a height of 3.5 mm, and included the epidermis, dermis, hair shaft, and follicle, as shown in Figure 1. The thicknesses of the epidermis and dermis were 0.07 mm and 3.43 mm, respectively. The hair was placed at the center of the model. The hair shaft was modeled as a tiny cylinder with a radius of 0.1 mm and a length of 3 mm, while the follicle had a radius of 0.5 mm and was located at a depth of 2.5 mm. In simulations, we applied laser beams with a radius of no more than 1 mm; therefore, the width of the model was deemed sufficient to account for the multiple scatterings of photons. The FEM model had a total of 800,000 tetrahedral elements.  
In the simulation protocol, we applied the following:  
* The simulated disc-shaped laser beam (1064 nm) with a radius of 1 mm was directed perpendicularly to the cross-section of the FEM model.
* The main outcomes of simulations were temperature distribution maps and the average temperature of the follicle.
* Simulations were undertaken by varying the fluence and pulse duration: we closely followed the published study and accordingly used the sequence of three pulses of fluence and pulse duration of 4.8 J/cm\({}^{2}\) and 2.5 ms, respectively, with a pause of 200 ms between each pulse when analyzing temperature distribution maps; when analyzing average temperature of the follicle, the fluence was between 5.0 J/cm\({}^{2}\) and 7.0 J/cm\({}^{2}\) with 0.5 J/cm\({}^{2}\) steps, each pulse duration of 2.5 ms, with a 200 ms pause between each pulse.  
Figure A5 shows the temperature distributions at \(y=0\) mm immediately after the first, second, and third pulses, respectively. Comparison with Figure 6 in the publication of [PERSON] et al. [14] shows excellent qualitative correspondence; the only difference is due to the slightly different shape of follicles; in our simulations, the follicle was of a spherical shape, while in the publication it was of a truncated sphere.  
The results of the analysis of the average temperature of the hair follicle are shown in Figure A6. Comparison with Figure 8 in the article by [PERSON] et al. [14] shows agreement within 5 percent. There are a couple of potential reasons for the relatively minor deviation: one is due to the different shapes of the follicle, and the other is due to different ways of calculating the average temperatures when using FEM and finite-difference models [14].