Finite Element Method (FEM) Modeling of Laser-Tissue Interaction during Hair Removal

: In this study, a comprehensive and realistic model of laser light interaction with skin and hair was constructed. The model was applied to study the characteristics of laser-tissue interaction for the deeply penetrating Nd:YAG laser. Three types of ﬁnite element method (FEM) models were developed. In the ﬁrst model, the hair shaft grew straight out of the follicle; in the second model, it grew at a variable angle; and in the third model, an array of hair was considered. The transport equation and heat diffusion equation were solved with the mesh-based Monte Carlo method and partial differential equations, respectively. The results of the simulations indicated that the area of necrosis increased with increasing ﬂuence; cooling had a limited effect on the extent of necrosis, particularly at a ﬂuence of 80 J/cm 2 . The thermal damage to hair follicles on the periphery of an irradiated array of hair may be insufﬁcient for achieving necrosis. The pulse itself and the short cooling-down period after the pulse contributed the most to the ﬁnal thermal damage to the hair follicle. The FEM modeling of laser-tissue interaction has proven to be a useful tool for studying the inﬂuence of different therapeutic parameters on the resulting hair and skin damage.


Introduction
In recent decades, laser technology has become an increasingly prevalent tool for hair removal, both for aesthetic and health-related purposes [1][2][3][4]. When compared to conventional practices such as shaving, waxing, chemical depilatories, and electrolysis, laser hair removal procedures appear to be fast, with long-lasting results and minimal side effects [5,6]. The mechanism of laser hair removal has been well described and is based on irreversible photothermal damage to hair follicles through the use of a light source at the red or infrared end of the spectrum [1]. Among the critical factors influencing the effectiveness of laser hair removal are the laser wavelength, duration of the pulse, fluence, and cooling method [1]. The main advantage of laser hair removal is selective photothermolysis, while most complications can be prevented by correctly selecting the laser type tailored to each patient and adjusting parameters and cooling [7]. It is imperative to have well-trained personnel with in-depth knowledge of the procedures to perform laser treatments effectively and safely [7].
The most frequently used laser hair removal systems are ruby laser (694 nm), alexandrite laser (755 nm), pulsed diode lasers (800, 810 nm), and Nd:YAG laser (1064 nm); in addition, intense pulsed light (IPL) sources with a broad spectrum (within 500-1200 nm) are also being used for hair removal [7]. In this article, we exclusively focus on the long-pulse Nd:YAG laser (1064 nm), which is characterized by the largest penetration within skin due to its low absorption and scattering, while at the same time the absorption in the hair follicle is sufficiently high to achieve irreversible damage to this area.
Ideally, the duration of the laser pulse should be shorter than the thermal relaxation time of the hair and longer than that of the epidermis in order to achieve effective hair removal while avoiding unnecessary skin damage [1]. However, since the epidermal and hair relaxation times both vary in the range of 1-100 ms, with thinner hair having a shorter relaxation time, this approach may not always be viable [8]. For example, one of the challenges involved with commercial diode laser hair removal systems is their very long pulse duration [1].
Most commonly, induced skin cooling is used to prevent epidermal damage. The fastest cooling method, albeit with certain risks involved, is a cryogenic cooling spray that evaporates immediately and very quickly removes heat from the epidermis without substantially changing the temperature in the deeper layers of the skin [9,10]. Safer methods include contact cooling with cold metal and induced cooling with cold air [11,12]. Most recently, a DMC (Dry Molecular Cooling) technology has been introduced that is based on a controlled fine water spray, combining the positive features of the cryogenic spray and cold air cooling [13].
Numerical modeling of laser-tissue interaction is useful as it can aid in optimizing the design and validation of laser hair removal protocols and systems by methodically assessing the influence of various parameters before submitting patients to unnecessary clinical trials. It is therefore somewhat surprising that, according to our knowledge, there has been only one published study [14] that applied an adequately realistic numerical modeling technique to simulate light and heat transport during laser hair removal. Specifically, Sun et al. [14] systematically simulated hair removal with an IPL device for wavelengths between 400 and 1200 nm, with fluences from 5.0-7.0 J/cm 2 and pulse lengths of 2.5-25 ms. Their approach was based on the three-dimensional modular adaptable grid model of Pfefer et al. [15] and contained epidermis and dermis, while the hair was approximated as a cylindrical sleeve with a truncated sphere.
The primary aim of the present study was to extend the formalism of Sun et al. [14] by adding an additional tissue type to the finite element method (FEM) model in order to be able to study hair removal characteristics for the deeply penetrating (up to~1 cm) Nd:YAG laser. The developed model was then used to study laser-tissue interaction depending on the exit angle of the hair relative to the skin surface and for a more realistic situation where multiple hairs are concurrently irradiated by a large spot-size laser beam. Additionally, the influence of externally applied cooling on the resulting skin and hair damage was evaluated.

FEM (Finite Element Method) Modelling
We developed our skin and hair modeling for the Nd:YAG laser using the finite element method (FEM), building upon previous research [14,15]. To ensure accurate results, we adapted the model's geometric shape to a cylinder instead of the conventional rectangular shape with sharp corners. This modification reduced the need for a larger number of volumetric elements while maintaining comparable accuracy. SALOME (v7) software was used to create mesh models. The process began by generating the geometry of the mesh, with careful consideration given to achieving a seamless connection between the surfaces to eliminate any gaps. Subsequently, the NETGEN [16] algorithm was employed to establish a surface mesh, which served as a foundation for generating the volumetric mesh.
The FEM Model #1 had a radius of 10 mm, a height of 8.5 mm, and included the epidermis [14,17], dermis [14,18], subcutis [19], hair shaft [20], and follicle [21], as shown in Figure 1. The thicknesses of the epidermis, dermis, and subcutis were 0.07 mm, 3.43 mm, and 5 mm, respectively. The hair was placed at the center of the model cylinder. The hair shaft was modelled as a tiny cylinder with a radius of 0.03 mm and a length of 2.5 mm, oriented along the long axis of the larger cylinder, while the follicle had a radius of 0.1 mm and was located at a depth of 2.5 mm [19]. In simulations, we applied laser beams with a radius of no more than 4 mm; therefore, the width of the model was deemed sufficient to account for the multiple scatterings of photons [14].
The FEM Model #1 had a radius of 10 mm, a height of 8.5 mm, and included the epidermis [14,17], dermis [14,18], subcutis [19], hair shaft [20], and follicle [21], as shown in Figure 1. The thicknesses of the epidermis, dermis, and subcutis were 0.07 mm, 3.43 mm, and 5 mm, respectively. The hair was placed at the center of the model cylinder. The hair shaft was modelled as a tiny cylinder with a radius of 0.03 mm and a length of 2.5 mm, oriented along the long axis of the larger cylinder, while the follicle had a radius of 0.1 mm and was located at a depth of 2.5 mm [19]. In simulations, we applied laser beams with a radius of no more than 4 mm; therefore, the width of the model was deemed sufficient to account for the multiple scatterings of photons [14]. Hair typically exits the skin at an angle, and for that purpose, FEM Model #2 was built, identical to FEM Model #1, with the only difference being that the hair shaft grew out of the follicle at different angles. In Figure 2, an example is depicted for the angle of 45°. In that specific case, the hair was 3.53 mm long, exiting the skin at the coordinates x = 2.5 mm and y = 0 mm (instead of x = 0 mm and y = 0 mm when considering straight hair). The exit point of the hair from the skin was thus still covered by a laser beam with a radius of 4 mm. The surface density of hair on human skin has on average 14-32 hair follicles per square centimeter [20], which we have taken into account with the even more realistic FEM Model #3, with an array of hair arranged in the shape of a cross ( Figure 3); this arrangement corresponded to a density of 18 follicles per square centimeter. Hair typically exits the skin at an angle, and for that purpose, FEM Model #2 was built, identical to FEM Model #1, with the only difference being that the hair shaft grew out of the follicle at different angles. In Figure 2, an example is depicted for the angle of 45 • . In that specific case, the hair was 3.53 mm long, exiting the skin at the coordinates x = 2.5 mm and y = 0 mm (instead of x = 0 mm and y = 0 mm when considering straight hair). The exit point of the hair from the skin was thus still covered by a laser beam with a radius of 4 mm.
The FEM Model #1 had a radius of 10 mm, a height of 8.5 mm, and included the epidermis [14,17], dermis [14,18], subcutis [19], hair shaft [20], and follicle [21], as shown in Figure 1. The thicknesses of the epidermis, dermis, and subcutis were 0.07 mm, 3.43 mm, and 5 mm, respectively. The hair was placed at the center of the model cylinder. The hair shaft was modelled as a tiny cylinder with a radius of 0.03 mm and a length of 2.5 mm, oriented along the long axis of the larger cylinder, while the follicle had a radius of 0.1 mm and was located at a depth of 2.5 mm [19]. In simulations, we applied laser beams with a radius of no more than 4 mm; therefore, the width of the model was deemed sufficient to account for the multiple scatterings of photons [14]. Hair typically exits the skin at an angle, and for that purpose, FEM Model #2 was built, identical to FEM Model #1, with the only difference being that the hair shaft grew out of the follicle at different angles. In Figure 2, an example is depicted for the angle of 45°. In that specific case, the hair was 3.53 mm long, exiting the skin at the coordinates x = 2.5 mm and y = 0 mm (instead of x = 0 mm and y = 0 mm when considering straight hair). The exit point of the hair from the skin was thus still covered by a laser beam with a radius of 4 mm. The surface density of hair on human skin has on average 14-32 hair follicles per square centimeter [20], which we have taken into account with the even more realistic FEM Model #3, with an array of hair arranged in the shape of a cross ( Figure 3); this arrangement corresponded to a density of 18 follicles per square centimeter. The surface density of hair on human skin has on average 14-32 hair follicles per square centimeter [20], which we have taken into account with the even more realistic FEM Model #3, with an array of hair arranged in the shape of a cross ( Figure 3); this arrangement corresponded to a density of 18 follicles per square centimeter. Appl

Laser-Tissue Interaction
The transport of photons in the tissue is the first of the two main physical processes occurring during laser hair removal. In the first step, we simulated the irradiation of each of our FEM models and, to that end, used the mesh-based Monte Carlo method (MMCM) [22] as an implementation algorithm for solving the transport equation. The MMCM algorithm assumed homogeneity of the optical properties of the tissue within each tetrahedral element and allowed for parallelization via several processor cores or a graphics card, which significantly accelerated the simulations [23]. Our software enabled varying the type, dimension, and direction of the laser beam. The optical properties of different tissue types used in simulations are summarized in Table 1. While a substantial amount of work has been carried out in determining the optical properties of human skin and hair [24,25], it has been noted [26] that the tabulated values of optical properties may have a range of values and may be determined only for specific wavelengths; additionally, the absorption coefficients of epidermis and hair have proven to be particularly difficult to assess because they contain melanin, which exhibits pronounced inter-subject variability [27][28][29][30][31][32]. For that reason, we have developed, verified, and piloted a protocol for in vivo measurement of absorption coefficients of epidermis and hair in human subjects, which is described in detail in Appendix A. Table 1. Optical properties of different tissue types used in simulations along with references; µa and µs are absorption and scattering coefficients, respectively; n is the refractive index; and g is an anisotropy factor. The second modeled physical process is heat diffusion within the tissue. The absorbance calculated by MMCM for each tissue type and each FEM element was converted into a heat source (W/mm 3 ). The calculated heat source was then fed into the heat diffusion equation, subject to the initial conditions, the boundary conditions between the tissue surfaces, and the convective thermal boundary condition with a constant heat transfer coefficient h, as specified in Appendix B [14]. Again, homogeneity of the thermal properties of the tissue within each tetrahedral FEM element was assumed. We solved the heat diffusion equation in the MATLAB software environment (R2019b) using the PDE toolbox. The "createpde" function was employed with the parameters "thermal" and "transient" to simulate the time-varying temperature distribution in a model. To simulate different

Laser-Tissue Interaction
The transport of photons in the tissue is the first of the two main physical processes occurring during laser hair removal. In the first step, we simulated the irradiation of each of our FEM models and, to that end, used the mesh-based Monte Carlo method (MMCM) [22] as an implementation algorithm for solving the transport equation. The MMCM algorithm assumed homogeneity of the optical properties of the tissue within each tetrahedral element and allowed for parallelization via several processor cores or a graphics card, which significantly accelerated the simulations [23]. Our software enabled varying the type, dimension, and direction of the laser beam. The optical properties of different tissue types used in simulations are summarized in Table 1. While a substantial amount of work has been carried out in determining the optical properties of human skin and hair [24,25], it has been noted [26] that the tabulated values of optical properties may have a range of values and may be determined only for specific wavelengths; additionally, the absorption coefficients of epidermis and hair have proven to be particularly difficult to assess because they contain melanin, which exhibits pronounced inter-subject variability [27][28][29][30][31][32]. For that reason, we have developed, verified, and piloted a protocol for in vivo measurement of absorption coefficients of epidermis and hair in human subjects, which is described in detail in Appendix A. Table 1. Optical properties of different tissue types used in simulations along with references; µ a and µ s are absorption and scattering coefficients, respectively; n is the refractive index; and g is an anisotropy factor. The second modeled physical process is heat diffusion within the tissue. The absorbance calculated by MMCM for each tissue type and each FEM element was converted into a heat source (W/mm 3 ). The calculated heat source was then fed into the heat diffusion equation, subject to the initial conditions, the boundary conditions between the tissue surfaces, and the convective thermal boundary condition with a constant heat transfer coefficient h, as specified in Appendix B [14]. Again, homogeneity of the thermal properties of the tissue within each tetrahedral FEM element was assumed. We solved the heat diffusion equation in the MATLAB software environment (R2019b) using the PDE toolbox. The "createpde" function was employed with the parameters "thermal" and "transient" to simulate the time-varying temperature distribution in a model. To simulate different therapeutic protocols, we varied the laser fluence and pulse duration, or forced cooling  Table 2. Table 2. Thermal properties of the tissues used in simulations, along with references: c p is specific heat capacity; k is thermal conductivity; and ρ is density.

Thermal Damage Estimation in the Tissues
Although thermal damage to the tissue depends on many complex mechanisms, it is accepted that the time course can be well approximated by a single process of protein denaturation. It is described by a simple kinetic first-order expression known as a damage integral Ω (x, y, z, t), based on the Arrhenius equation [33][34][35], The values for the Arrhenius constants A and E a were assumed to be the same for all tissue types and were set at 2.9 × 10 37 s −1 and 240 J/mol [36]. When the value of Ω (x, y, z, t) equaled at least 1, then at least 63% of cells within the given tetrahedral FEM element were deemed to be dead. This meant that the skin tissue within that element was irreversibly damaged.

Simulation Protocols of Thermal Damage in the Tissue
Before applying simulation protocols, we thoroughly tested our FEM modeling approach by comparing the results of simulations with the publication of Sun et al. [14], as outlined in Appendix C.
In the first simulation protocol (Protocol #1), we applied the following: -The simulated laser beam, with a wavelength of 1064 nm and a disc-shaped radius of 4 mm, was directed perpendicularly to the short-axis center of FEM Model #1, FEM Model #2, and FEM Model #3. - The parametric study was undertaken by varying the fluence and pulse duration: we closely followed the clinical study [37] and accordingly used combinations of fluence and pulse duration of (i) 50 J/cm 2 and 25 ms, (ii) 60 J/cm 2 and 50 ms, and (iii) 80 J/cm 2 and 50 ms. -We explicitly considered the (i) mode without forced cooling, with a value of the free convection coefficient h of 10 W/m 2 K [14,38] and ambient temperature T = 22 • C for all three FEM models, and the (ii) mode with forced cooling by means of a cryogenic spray with a convection coefficient h of 5000 W/m 2 K [39] and spray temperature T = −50 • C [39] for FEM Model #1 and FEM Model #2 only. The forced cooling by cryogenic spray consisted of 10 ms of preliminary cooling, with a 5-ms pause immediately after the pulse and 15 ms of cooling thereafter, as outlined in the clinical study [37]. - The main outcomes of the simulations were temperature and thermal damage distributions after the completion of the irradiation.
In the second simulation protocol (Protocol #2), we applied a wider laser beam with a diameter of 6 mm and a combination of fluence of 80 J/cm 2 and a pulse duration of 50 ms; in this protocol, we considered only the mode without forced cooling. We used the more realistic models, FEM Model #2 and FEM Model #3. The main outcome was a time-dependent curve of cumulative hair-follicle damage for each of the two FEM models.

Optimization of the Laser Beam
The Monte Carlo methodology is a stochastic technique with random sampling, and thus its accuracy inherently depends on the number of simulated photons. To address this attribute, we calculated the absorbance in the hair follicle of the FEM Model #1 for different numbers of photons, ranging from 10 4 to 10 9 . The simulated laser beam was disc-shaped with a radius of 4 mm and directed perpendicularly to the short-axis center of the FEM Model #1. The average absorbance value of the entire hair follicle was determined by interpolating and spatially integrating values at its vertices. The mean value over the follicle was evaluated ten times (from ten independent simulations), and then the average and standard deviations were calculated, which are presented in Figure 4. The average absorbance was normalized, meaning that the source power was set at 1 W.

Optimization of the Laser Beam
The Monte Carlo methodology is a stochastic technique with random sampling, and thus its accuracy inherently depends on the number of simulated photons. To address this attribute, we calculated the absorbance in the hair follicle of the FEM Model #1 for different numbers of photons, ranging from 10 4 to 10 9 . The simulated laser beam was discshaped with a radius of 4 mm and directed perpendicularly to the short-axis center of the FEM Model #1. The average absorbance value of the entire hair follicle was determined by interpolating and spatially integrating values at its vertices. The mean value over the follicle was evaluated ten times (from ten independent simulations), and then the average and standard deviations were calculated, which are presented in Figure 4. The average absorbance was normalized, meaning that the source power was set at 1 W. As can be clearly seen in Figure 4, average hair-follicle absorbance converged with an increasing number of photons, while the number of photons beyond 10 7 had virtually no impact on the results. Based on these results, we decided to use 10 8 photons as the optimal quantity in all simulations.

Optimization of the FEM Model
The FEM models were constructed so that the density of the tetrahedral elements close to the hair shaft and follicle was much higher than the density on the periphery. A high density of tetrahedral elements was also present in the whole epidermis, while in the dermis, far from the hair, the density of tetrahedral elements was lower, with the lowest in the subcutis. The reason for this is that the description of a homogeneous volume, for example the subcutis far from the hair, does not require as many elements as, e.g., a thin epidermis or a tiny hair shaft and follicle. We found that the following ratio in the number of tetrahedral elements between the hair shaft/follicle (Nhair), epidermis (Nepidermis), dermis (Ndermis), and subcutis (Nsubcutis) would be appropriate: Nhair:Nepidermis:Ndermis:Nsubcutis = 10:10:100:1.
To estimate the minimal number of tetrahedral elements, we simulated the average temperature within the hair follicle using a fluence of 5 J/cm 2 and a sequence of three As can be clearly seen in Figure 4, average hair-follicle absorbance converged with an increasing number of photons, while the number of photons beyond 10 7 had virtually no impact on the results. Based on these results, we decided to use 10 8 photons as the optimal quantity in all simulations.

Optimization of the FEM Model
The FEM models were constructed so that the density of the tetrahedral elements close to the hair shaft and follicle was much higher than the density on the periphery. A high density of tetrahedral elements was also present in the whole epidermis, while in the dermis, far from the hair, the density of tetrahedral elements was lower, with the lowest in the subcutis. The reason for this is that the description of a homogeneous volume, for example the subcutis far from the hair, does not require as many elements as, e.g., a thin epidermis or a tiny hair shaft and follicle. We found that the following ratio in the number of tetrahedral elements between the hair shaft/follicle (Nhair), epidermis (Nepidermis), dermis (Ndermis), and subcutis (Nsubcutis) would be appropriate: Nhair:Nepidermis:Ndermis:Nsubcutis = 10:10:100:1.
To estimate the minimal number of tetrahedral elements, we simulated the average temperature within the hair follicle using a fluence of 5 J/cm 2 and a sequence of three consecutive pulses with a duration of 2.5 ms and a pause of 200 ms between each pulse. The average temperature was calculated in the same manner as described above. Figure 5 shows the dependence of the average temperature of the hair follicle immediately after the first, second, and third pulses, respectively, on the number of tetrahedral FEM elements. We can conclude that the average temperature of the hair follicle starts converging at roughly 500,000 elements, which was taken as the minimal number of elements in our simulations. FEM  consecutive pulses with a duration of 2.5 ms and a pause of 200 ms between each pulse. The average temperature was calculated in the same manner as described above. Figure 5 shows the dependence of the average temperature of the hair follicle immediately after the first, second, and third pulses, respectively, on the number of tetrahedral FEM elements. We can conclude that the average temperature of the hair follicle starts converging at roughly 500,000 elements, which was taken as the minimal number of elements in our simulations. FEM Model #1 thereby had 422,742 elements, FEM Model #2 had 453,775 elements, and FEM Model #3 had 586,633 elements. For the selected three models, the simulation speed of the laser-tissue interaction on the Intel i5-8250u CPU was approximately 50,000 photons per second, while the single heat diffusion simulation took approximately 3 h to complete. Figure 6 shows temperature distribution maps for FEM Model #1 immediately after irradiation for all three combinations of fluence and pulse duration, both without and with forced cooling. As expected, the maximal temperature (Table 3) occurred at the very tip of the hair shaft due to the hair having a larger absorption coefficient than the epidermis and dermis. The largest gradient was observed along the hair shaft. Interestingly, the temperature in the hair follicle was higher than in the hair shaft just above it, which can be attributed to the hair follicle absorbing a much larger number of scattered photons than the tiny hair shaft. All these observations are congruent with those of Sun et al. [14]. The combination of fluence and pulse duration had a non-intuitive effect, as the maximum temperature increased with increasing fluence but decreased with increasing pulse duration. The maximum temperature at the tip of the hair shaft was thus higher for a fluence of 50 J/cm 2 and a pulse duration of 25 ms than for a fluence of 60 J/cm 2 and a pulse duration of 50 ms. This happened because the difference in irradiation power normalized to the surface area of the laser beam (2000 W/cm 2 in the first case and 1200 W/cm 2 in the second case) was not compensated by delivering 17% more energy in twice the amount of time in the second case. The temperature was highest when irradiating the model with a fluence For the selected three models, the simulation speed of the laser-tissue interaction on the Intel i5-8250u CPU was approximately 50,000 photons per second, while the single heat diffusion simulation took approximately 3 h to complete. Figure 6 shows temperature distribution maps for FEM Model #1 immediately after irradiation for all three combinations of fluence and pulse duration, both without and with forced cooling. As expected, the maximal temperature (Table 3) occurred at the very tip of the hair shaft due to the hair having a larger absorption coefficient than the epidermis and dermis. The largest gradient was observed along the hair shaft. Interestingly, the temperature in the hair follicle was higher than in the hair shaft just above it, which can be attributed to the hair follicle absorbing a much larger number of scattered photons than the tiny hair shaft. All these observations are congruent with those of Sun et al. [14]. The combination of fluence and pulse duration had a non-intuitive effect, as the maximum temperature increased with increasing fluence but decreased with increasing pulse duration. The maximum temperature at the tip of the hair shaft was thus higher for a fluence of 50 J/cm 2 and a pulse duration of 25 ms than for a fluence of 60 J/cm 2 and a pulse duration of 50 ms. This happened because the difference in irradiation power normalized to the surface area of the laser beam (2000 W/cm 2 in the first case and 1200 W/cm 2 in the second case) was not compensated by delivering 17% more energy in twice the amount of time in the second case. The temperature was highest when irradiating the model with a fluence of 80 J/cm 2 and a pulse duration of 50 ms. In this case, the irradiation power, normalized to the surface of the laser beam, equaled 1600 W/cm 2 , but the higher energy and longer pulse were sufficient to deliver the highest temperature among the three combinations. These observations were corroborated by the maximum temperatures calculated within the follicle (Table 3). We observed that cooling had no impact on the follicle temperature. of 80 J/cm 2 and a pulse duration of 50 ms. In this case, the irradiation power, normalized to the surface of the laser beam, equaled 1600 W/cm 2 , but the higher energy and longer pulse were sufficient to deliver the highest temperature among the three combinations. These observations were corroborated by the maximum temperatures calculated within the follicle (Table 3). We observed that cooling had no impact on the follicle temperature.   Figure 7 shows distribution maps of thermal damage Ω for FEM Model #1 immediately after irradiation for three combinations of fluence and pulse duration, for modes without forced cooling and with forced cooling. As expected, the area of thermal necrosis (where Ω exceeds 1) increased with increasing fluence, while cooling led to the area of necrosis in the epidermis being substantially smaller compared to that without cooling, particularly for the fluence of 50 J/cm 2 with a pulse duration of 25 ms and the fluence of 60 J/cm 2 with a pulse duration of 50 ms. It is notable that the area of thermal necrosis along the hair shaft toward the hair follicle was not connected, reflecting the temperature gradient along the hair shaft and the fact that the temperature was higher in the follicle than in the hair shaft just above it ( Figure 6).   Figure 7 shows distribution maps of thermal damage Ω for FEM Model #1 immediately after irradiation for three combinations of fluence and pulse duration, for modes without forced cooling and with forced cooling. As expected, the area of thermal necrosis (where Ω exceeds 1) increased with increasing fluence, while cooling led to the area of necrosis in the epidermis being substantially smaller compared to that without cooling, particularly for the fluence of 50 J/cm 2 with a pulse duration of 25 ms and the fluence of 60 J/cm 2 with a pulse duration of 50 ms. It is notable that the area of thermal necrosis along the hair shaft toward the hair follicle was not connected, reflecting the temperature gradient along the hair shaft and the fact that the temperature was higher in the follicle than in the hair shaft just above it ( Figure 6). Appl Figure 8 shows distribution maps of thermal damage for FEM Model #2 immediately after irradiation for three combinations of fluence and pulse duration, for modes without forced cooling and with forced cooling. The area of necrosis in the hair follicle is smaller than in FEM Model #1, and the discontinuity of thermal necrosis along the hair shaft is more extensive. The main reason for this was the lower temperature achieved in the hair follicle due to the hair shaft being longer and at an angle relative to the incoming laser beam. Table 4 reveals that the maximum follicle temperature in FEM Model #2 was consistently lower than in FEM Model #1.  Figure 8 shows distribution maps of thermal damage for FEM Model #2 immediately after irradiation for three combinations of fluence and pulse duration, for modes without forced cooling and with forced cooling. The area of necrosis in the hair follicle is smaller than in FEM Model #1, and the discontinuity of thermal necrosis along the hair shaft is more extensive. The main reason for this was the lower temperature achieved in the hair follicle due to the hair shaft being longer and at an angle relative to the incoming laser beam. Table 4 reveals that the maximum follicle temperature in FEM Model #2 was consistently lower than in FEM Model #1. Figure 9 shows temperature and thermal damage distribution maps for FEM Model #3 immediately after irradiation for three combinations of fluence and pulse duration, in a mode without forced cooling. The non-central hair shafts and follicles heated up substantially less than the central hair shaft and follicle; as a result, for a fluence of 50 J/cm 2 and a pulse duration of 25 ms, no necrosis occurred in the non-central hair follicles. Table 5 indicates that the maximum temperature of the central follicle was lower than in FEM Model #1 and FEM Model #2 and that a non-central follicle had a 10 • C lower maximum temperature than the central follicle.   Figure 9 shows temperature and thermal damage distribution maps for FEM Model #3 immediately after irradiation for three combinations of fluence and pulse duration, in a mode without forced cooling. The non-central hair shafts and follicles heated up substantially less than the central hair shaft and follicle; as a result, for a fluence of 50 J/cm 2 and a pulse duration of 25 ms, no necrosis occurred in the non-central hair follicles. Table  5 indicates that the maximum temperature of the central follicle was lower than in FEM Model #1 and FEM Model #2 and that a non-central follicle had a 10 °C lower maximum temperature than the central follicle.

Time-Dependent Curve of Cumulative Hair-Follicle Damage-Protocol #2
When inspecting the results of cumulative heat damage to the hair follicle for FEM Model #2 and FEM Model #3, shown in Figure 10, we notice that the pulse itself and a short cooling-down period after the pulse contributed the most to the final thermal damage to the hair follicle; in our case, this amounted to 50 ms of the pulse duration and about 20 ms after the pulse, i.e., 70 ms in total. We can hypothesize that the heat damage to an array of hair may have been insufficient to reach necrosis of the hair follicle, which is in line with observations in FEM Model #3 under Protocol #1 (Figure 9).

Discussion
This is the first report of an approach to realistic modeling of skin and hair in which laser beam-tissue interaction was studied within an FEM framework. A modeling approach allowed us to examine temperature and thermal damage distributions following

Discussion
This is the first report of an approach to realistic modeling of skin and hair in which laser beam-tissue interaction was studied within an FEM framework. A modeling approach allowed us to examine temperature and thermal damage distributions following irradiation at a wavelength of 1064 nm for three combinations of fluence and pulse duration, with and without forced cooling. Our approach also accounted for the presence of angled hair and hair arrays and is comprehensive enough that it could be applied to further clinical research studies.
In this paper, we have undertaken a parametric study by varying light fluence and pulse duration, closely following the clinical study [37]. We have found that for all combinations of fluence and pulse duration, at least a partial necrosis of the hair follicle at a depth of 2.5 mm occurred, which corroborated the outcomes of the clinical study [37]. Increasing fluence resulted in an increased area of necrosis, particularly in the epidermis and dermis, which again agreed with the clinical study where an increase in the likelihood of damage to the epidermis and dermis was demonstrated, typically manifesting as erythema and perifocal edema [37]. The use of forced cooling made the area of necrosis less extensive, especially in the epidermis, which again agreed with the findings of the clinical study [37].
We have shown that the hair follicle of a hair exiting the skin surface at an angle undergoes less thermal damage than the follicle of a straight hair. In addition, we observed that the thermal exposure of hair follicles on the periphery of an irradiated array of hair may be insufficient for reaching necrosis. This is possibly a partial explanation for why laser hair removal does not remove all of the hair, since during the hair removal procedure, the laser beam is not centered on each hair separately but on a given area. Finally, our simulation study indicated that the pulse itself and a short relaxation time after the pulse contributed the most to the final thermal damage to the hair follicle.
We have also demonstrated that the FEM model of light-tissue interaction in combination with measurements of surface temperature has proven useful in assessing the absorption coefficients of human epidermis and hair. We attempted to develop, verify, and pilot a protocol for in vivo assessment of absorption coefficients of epidermis and hair in humans that could efficiently account for inter-subject variability. As the protocol of our study is simple and comprehensive enough, it may be potentially helpful in a clinical setting.
Because pulses significantly shorter than the hair's thermal relaxation time are expected to provide the highest temperatures, the initial hair removal lasers were Q-switched Nd:YAG lasers with nanosecond-long pulse durations [40]. However, since the achievable pulse energies of Q-switched lasers limited the fluences to below 2-4 J/cm 2 for spot radiuses of 6-4 mm, and typically required using a topical carbon-based solution to enhance hair absorption, most of today's hair removal procedures are performed with lasers operating in the millisecond (1-50 ms) pulse duration regime [1,3,41,42]. The longer Nd:YAG laser pulse durations, in the range of 25-50 ms, as used in our simulations, are often a preferred choice in order to avoid overheating of the epidermis with a shorter thermal relaxation time. However, our model can be easily adapted to incorporate shorter pulses. Another potential extension of our FEM modeling approach would be to simulate laser treatment of keratinocyte carcinoma [43].
An obvious extension of our study would be to consider thinner hair at different depths and different dimensions of the epidermis, dermis, and subcutis. Another possibility would be to vary skin and hair types. While our study specifically focused on the Nd:YAG laser, it could be easily applied to other lasers operating at shorter wavelengths and then compared to the results of corresponding clinical studies. In our study, we used the well-established Arrhenius damage integral [44], however, we could also consider other models of thermal damage [45]. Further, the model could be applied to study hair removal while taking into account the recently reported "avalanche-type" enhancement of hair temperature increase under repetitive pulse exposures [46]. Our ultimate goal is that the developed FEM model of laser beam-tissue interaction would also prove useful in a clinical setting by aiding the optimal selection of treatment parameters during hair-removal procedures. Using the novel hyperspectral imaging modality [47] could aid in, e.g., collecting information about the quantity of melanin and estimating epidermal thickness, which would in turn be helpful for selecting optimal laser treatment parameters.
In this study, we focused on the Nd:YAG laser, which is known for its ability to treat darker skin tones since the longer wavelength bypasses the melanin in the skin. However, our approach could also be applied to evaluate the effects of other hair removal laser wavelengths and the effects of multiple wavelengths present during IPL treatments. The IPL (Intense Pulsed Light) uses broad-spectrum light to target and remove unwanted hair and is particularly efficient in removing dark hair on fair to medium skin tones. The method is fast and inexpensive, but the major drawback is that it may lead to pigmentation of the darker skin [48]. In spite of our focusing on the Nd:YAG wavelength only, we believe that the major conclusions of our study related to the dependence of the efficacy of the hair removal on the hair density and the angle of the hair shaft apply as well to other light-based hair removal methods.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request. (protocol number 0120-365/2017-9); the subject signed an informed consent form. Temperature rises on the skin surface were recorded in different areas of freshly shaved forearm skin, hand-cleaned with cotton wool and ethanol, and finally air-dried. This anatomical location was chosen because of its easy accessibility for the experimental setup and ease of handling. Experiments were also undertaken using a PlatSil SiliGlass type C tissuemimicking phantom with known optical and thermal properties [49]. All personnel present during the laboratory experiment wore appropriate eye protection.

Appendix A.3. Simulations
To simulate the transport of photons and heat diffusion in the tissue, we used the finite element method framework (FEM) in combination with the mesh-based Monte Carlo method (MMCM), as described in the main body of this article and Appendix B. The FEM model of the PlatSil SiliGlass type C tissue-mimicking phantom had the following optical and thermal properties. The absorption and scattering coefficients µ a and µ s were 0.035 mm −1 and 5.9 mm −1 , respectively; the anisotropy factor g was set to 0.7; the refractive index n was 1.44 [49], while the specific heat c p , thermal conductivity k, and density ρ were 703 J/kgK, 1.38 W/mK, and 2203 kg/m 3 , respectively [49]. The FEM model of the phantom was homogeneous, i.e., had only one tissue type and had 8,697 tetrahedral elements. To assess the absorption coefficients of the epidermis and hair, we used the realistic FEM Model #1 with the optical and thermal properties of different tissue types specified in Tables 1 and 2, respectively, with the µ a for the epidermis and hair obviously missing.
Verification of the protocol using the PlatSil SiliGlass type C tissue-mimicking phantom. Figure A1 shows the temperature distribution map on the surface of the PlatSil Sili-Glass type C phantom immediately following irradiation with a fluence of 50 J/cm 2 and a pulse duration of 30 ms. In the center of the irradiated area, the temperature increase amounted to more than 30 K. The time course of the temperature change in the center of irradiation following the pulse is shown in Figure A1, where we can see that experimental measurements agree well with simulation results. Appl. Sci. 2023, 13, x FOR PEER REVIEW 15 of 20 Figure A1 shows the temperature distribution map on the surface of the PlatSil SiliGlass type C phantom immediately following irradiation with a fluence of 50 J/cm 2 and a pulse duration of 30 ms. In the center of the irradiated area, the temperature increase amounted to more than 30 K. The time course of the temperature change in the center of irradiation following the pulse is shown in Figure A1, where we can see that experimental measurements agree well with simulation results.   Figure A2 shows the time course of the temperature change in the human subject following irradiation of (i) hair on the skin surface and (ii) skin at two different locations distant from visible hair. For the case with a fluence of 20 J/cm 2 and a duration of 30 ms, the temperature increase was 16 °C, while in the case with a fluence of 40 J/cm 2 and a duration of 30 ms, the corresponding increase was 41 °C.  Appendix A. 5

. Simulations
As the absorption coefficient mainly affects the temperature increase on the surface of the skin during irradiation, we decided to assess the absorption coefficient of skin and hair in such a way that the experimental measurements and simulated results were matched at the moment of maximum temperature increase. Figure A3 shows the measured and simulated time course of the temperature change in both subjects following irradiation of hair and skin; we exhibit a fluence of 40 J/cm 2 and a duration of 30 ms, while a fluence of 20 J/cm 2 and a duration of 30 ms gave similar results. As we can see, the measured and simulated results for the epidermis matched very well, and the absorption coefficient of the epidermis µa was 0.069 mm −1 . For hair, the correspondence between the measured and simulated time course of the temperature change was worse; we attribute this to the fact that the absorption coefficient regulates only the temperature increase but has no effect on relaxation, which may very much depend on the thickness of the hair and the proximity of adjacent hairs. The absorption coefficient of the hair µa was estimated at 0.56 mm −1 .

. Simulations
As the absorption coefficient mainly affects the temperature increase on the sur of the skin during irradiation, we decided to assess the absorption coefficient of skin hair in such a way that the experimental measurements and simulated results w matched at the moment of maximum temperature increase. Figure A3 shows the m ured and simulated time course of the temperature change in both subjects following radiation of hair and skin; we exhibit a fluence of 40 J/cm 2 and a duration of 30 ms, w a fluence of 20 J/cm 2 and a duration of 30 ms gave similar results. As we can see, the m ured and simulated results for the epidermis matched very well, and the absorption c ficient of the epidermis µa was 0.069 mm −1 . For hair, the correspondence between measured and simulated time course of the temperature change was worse; we attrib this to the fact that the absorption coefficient regulates only the temperature increase has no effect on relaxation, which may very much depend on the thickness of the hair the proximity of adjacent hairs. The absorption coefficient of the hair µa was estimate 0.56 mm −1 .

Appendix B. Heat Diffusion Equation
Formalism of the heat diffusion process was described for each tissue type and e FEM element as follows [14]:   [14]. The model includes epidermis, dermis, hair shaft, and follicle and has a radius of 10 mm and a height of 8.5 mm; the thickness of epidermis and dermis is 0.07 mm and 3.43 mm, respectively; the cylindrical hair shaft has a radius of 0.1 mm and a length of 3 mm; the follicle has a radius of 0.5 mm and is located at a depth of 2.5 mm.
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. Figure A4. A three-dimensional representation of the Finite-Element Method (FEM) model of skin and hair for validating our simulations with the results of Sun et al. [14]. The model includes epidermis, dermis, hair shaft, and follicle and has a radius of 10 mm and a height of 8.5 mm; the thickness of epidermis and dermis is 0.07 mm and 3.43 mm, respectively; the cylindrical hair shaft has a radius of 0.1 mm and a length of 3 mm; the follicle has a radius of 0.5 mm and is located at a depth of 2.5 mm.
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 Sun 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. Compare with Figure 6 in the article published by Sun et al. [14].
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 Sun 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]. Figure A5. Temperature distribution maps shown in the y = 0 plane for the FEM model immediately (left) after the first pulse, (middle) after the second pulse, and (right) after the third pulse (fluence per pulse is 4.8 J/cm 2 and the duration of pulse is 2.5 ms, with a 200 ms pause between each pulse). Compare with Figure 6 in the article published by Sun et al. [14].
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 Sun 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].
The results of the analysis of the average temperature of the hair follicle are in Figure A6. Comparison with Figure 8 in the article by Sun et al. [14] shows agr within 5 percent. There are a couple of potential reasons for the relatively minor de one is due to the different shapes of the follicle, and the other is due to different calculating the average temperatures when using FEM and finite-difference mode Figure A6. Effect of fluence (x axis) on average follicle temperature (y axis) immediately first pulse (purple line), second pulse (blue line), and third pulse (red line). Compare with in the article by Sun et al. [14].