A New Nonlinear Photothermal Iterative Theory for Port-Wine Stain Detection

The development of appropriate photothermal detection of skin diseases to meet complex clinical demands is an urgent challenge for the prevention and therapy of skin cancer. An extensive body of literature has ignored all high-order harmonics above the second order and their influences on low-order harmonics. In this paper, a new iterative numerical method is developed for solving the nonlinear thermal diffusion equation to improve nonlinear photothermal detection for the noninvasive assessment of the thickness of port-wine stain (PWS). First, based on the anatomical and structural properties of skin tissue of PWS, a nonlinear theoretical model for photothermal detection is established. Second, a corresponding nonlinear thermal diffusion equation is solved by using the new iterative numerical method and taking into account harmonics above the second-order and their effects on lower-order harmonics. Finally, the thickness and excitation light intensity of PWS samples are numerically simulated. The simulation results show that the numerical solution converges fasterand the physical meaning of the solution is clearerwith the new method than with the traditional perturbation method. The rate of change in each harmonic with the sample thickness for the new method is higher than that for the conventional perturbation method, suggesting that the proposed numerical method may provide greater detection sensitivity. The results of the study provide a theoretical basis for the clinical treatment of PWS.


Introduction
Port-wine stain (PWS), also known as nevus flammeus, is a congenital telangiectasia deformity. It is the most common type of benign vascular malformation and is difficult to cure [1,2]. Wine discoloration often occurs on the head, face, and neck, and severe cases are accompanied by overgrowth of soft tissues and bones in the lesion area, resulting in local enlargement and deformation [3]. These lesions greatly affect the patient's appearance, decrease their quality of life, and cause considerable mental stress [4,5]. Therefore, early and effective intervention is particularly important.
Currently, the evaluation and prediction of PWS treatment consist of invasive and noninvasive approaches. Although biopsy has long been considered the gold standard for treatment, it is invasive and not widely performed. Noninvasive treatments include chromatography [6], dermoscopy [7], high−frequency ultrasound [8,9], and laser scatter imaging [10]. However, none of these commonly used imaging techniques provide adequate imaging depth and contrast to accurately assess PWS. A recent trial has shown that the use of photoacoustic techniques for the clinical evaluation of PWS disorders provides a new method for the quantitative evaluation of PWS [11].
Due to combining the advantages of both deep penetration provided by ultrasound imaging [12,13] and high contrast provided by optical imaging [14], photoacoustic technology [15,16] has become a research frontier and hot spot in the field of biomedical imaging. 2 of 12 Most studies of photoacoustic techniques ignore the effect of the local temperature increase of the medium caused by light absorption on the thermodynamic parameters of the medium (e.g., thermal conductivity, density, and isobaric specific heat capacity) and assume that the thermodynamic parameters are constant. However, a statistically significant increase in the local temperature can change the values of thermodynamic parameters of the medium, and contribute to nonlinear photoacoustic conversion. The nonlinearity describing the thermal conductivity problem of laser irradiated tissue can be caused by various physical reasons, e.g., laser−induced formation of bubbles due to temperature dependence of gas solubility [17,18]; temperature dependence of thermodynamic parameters [19], etc. The nonlinear photoacoustic effect has attracted increasing attention as a possible means of selective detection of contrast agents by heat accumulation and local temperature increase thus enhancing the photoacoustic signal, and it is necessary to consider the nonlinearity of the thermal parameters in the thermal diffusion equation [20,21].
Therefore, this paper investigates the theory of thermal field imaging of PWS using nonlinear photoacoustic effect in the frequency domain by introducing nonlinear thermal conductivity coefficients. Based on previous work [22,23], a new semianalytic numerical iterative method is proposed in this paper. First, the temperature field is expanded in a Fourier series in the frequency domain to separate time variables and spatial coordinates. Then, the nonlinear diffusion equation is solved by selecting the appropriate high−frequency harmonics according to the specific requirements for calculation accuracy. Finally, the thickness and excitation light intensity are numerically simulated in light absorbers of different thicknesses using the new iterative numerical method and the conventional perturbation method. The results show that the solution bythe numerical method has greater sensitivity and bandwidth than that of the perturbation method and can better distinguish between PWS samples of different thicknesses. This work extends the application of nonlinear thermal field theory in clinical medicine and contributes to a better understanding of PWS in lesions during different stages of development.

Theoretical Model
The skin is the largest and most important tissue in the human body. Anatomically, the skin can be divided into the epidermis, dermis, hypodermis, and muscle layers, as shown in Figure 1a. The epidermis consists of the high-fat, low-water stratum corneum and the melanin-containing living epidermis. Similarly, the dermis has two sublayers: the papillary dermis and the reticular dermis, which contain two vascular plexuses; the upper and deep blood plexuses are located in the upper and lower reticular layers of the dermis. The subcutaneous tissue consists mainly of fat cells. technology [15,16] has become a research frontier and hot spot in the field of biomedical imaging. Most studies of photoacoustic techniques ignore the effect of the local temperature increase of the medium caused by light absorption on the thermodynamic parameters of the medium (e.g., thermal conductivity, density, and isobaric specific heat capacity) and assume that the thermodynamic parameters are constant. However, a statistically significant increase in the local temperature can change the values of thermodynamic parameters of the medium, and contribute to nonlinear photoacoustic conversion. The nonlinearity describing the thermal conductivity problem of laser irradiated tissue can be caused by various physical reasons, e.g., laser−induced formation of bubbles due to temperature dependence of gas solubility [17,18]; temperature dependence of thermodynamic parameters [19], etc. The nonlinear photoacoustic effect has attracted increasing attention as a possible means of selective detection of contrast agents by heat accumulation and local temperature increase thus enhancing the photoacoustic signal, and it is necessary to consider the nonlinearity of the thermal parameters in the thermal diffusion equation [20,21]. Therefore, this paper investigates the theory of thermal field imaging of PWS using nonlinear photoacoustic effect in the frequency domain by introducing nonlinear thermal conductivity coefficients. Based on previous work [22,23], a new semianalytic numerical iterative method is proposed in this paper. First, the temperature field is expanded in a Fourier series in the frequency domain to separate time variables and spatial coordinates. Then, the nonlinear diffusion equation is solved by selecting the appropriate high−frequency harmonics according to the specific requirements for calculation accuracy. Finally, the thickness and excitation light intensity are numerically simulated in light absorbers of different thicknesses using the new iterative numerical method and the conventional perturbation method. The results show that the solution bythe numerical method has greater sensitivity and bandwidth than that of the perturbation method and can better distinguish between PWS samples of different thicknesses. This work extends the application of nonlinear thermal field theory in clinical medicine and contributes to a better understanding of PWS in lesions during different stages of development.

Theoretical Model
The skin is the largest and most important tissue in the human body. Anatomically, the skin can be divided into the epidermis, dermis, hypodermis, and muscle layers, as shown in Figure 1a. The epidermis consists of the high-fat, low-water stratum corneum and the melanin-containing living epidermis. Similarly, the dermis has two sublayers: the papillary dermis and the reticular dermis, which contain two vascular plexuses; the upper and deep blood plexuses are located in the upper and lower reticular layers of the dermis. The subcutaneous tissue consists mainly of fat cells. In human dermatology, PWS is one of the most common benign tumors involving vascular malformations, with an incidence of between approximately 0.3 and 0.5% in the general population. PWS are mainly located in the papillary layer of the dermis and the upper layer of the reticular layer, with a diameter of approximately 0.01~0.15 mm and a thickness of approximately 0.001~1.5 mm [24]. Lesions can expand in size to cover a larger dermal area over time, as shown in Figure 1b.
Based on the anatomical and structural properties of skin tissue lesions, this paper constructs a nonlinear photothermal detection model for PWS, as shown in Figure 1c. Based on the optical attenuation, imaging depth, and other information reported by Chen et al., a laser pulse with a wavelength of 840 nm was chosen for PWS detection in this study. The parameters of each skin layer at a wavelength of 840 nm are shown in Table 1 [11,[25][26][27]. Table 1. Thickness, optical and thermal parameters of the skin model at 840 nm [11,[25][26][27].

Nonlinear Thermal Diffusion Equation
A wave pulsed laser I(t) = 2 ln (2) π ω P τ exp −4 ln(2)( t−t 0 τ ) 2 irradiates the sample, and light is absorbed by the sample and converted into heat, generating a photoacoustic signal through thermal expansion, where τ is the pulse width, t 0 is the pulse center, and w P is the luminous flux (mJ/cm 2 ). The fundamental equation of the photoacoustic imaging theory is based on the thermal diffusion equation; therefore, the thermal diffusion equation and boundary conditions in the sample are [28,29]: where T s (z, t), T g (z, t) are the temperature rise and fall of the sample and air, respectively, α = K s /(ρC P ) is the thermal diffusion coefficient, ρ and C p are the density and isobaric specific heat capacity of the sample, respectively, H = h + 4εσT 3 , h is the convection coefficient and 4εσT 3 is the heat radiation exchange term. According to 4εσT 3 h [30], the heat radiation exchange term can be neglected. Subscripts behind the comma mean the calculation of partial derivatives for the corresponding subscripts. The endothermic source in the sample is Q(z, t), which can be expressed as [28]: where β is the bulk absorption coefficient of the sample, β is the surface absorption rate of the sample, ω is the modulation angular frequency, and K s , K g is the thermal conductivity of the sample and air. In linear theory, K s is generally assumed to be constant; however, when the local temperature rise is significantly higher than the average temperature, K s this assumption is not valid. Therefore, this paper introduces nonlinear thermal conductivity coefficients into the nonlinear thermal diffusion equation to investigate the theory of nonlinear photothermal imaging.
In nonlinear thermal diffusion theory, the dependence of thermal conductivity on the temperature is usually considered to be linear to simplify the problem [31]: where K 0 is the thermal conductivity of the sample at steady−state temperature, b is the temperature coefficient of thermal conductivity, and in general, b 1 [32]. Substituting Equation (4) into Equation (1), the one−dimensional nonlinear thermal diffusion equation is obtained as follows: where α b = K 0 /(ρC P ). Subscripts z and t behind the comma mean the calculation of partial derivatives for the corresponding subscripts.

Iterative Numerical Method for Solving the Nonlinear Heat Diffusion Equation
Because it is difficult to obtain a general solution of the one−dimensional nonlinear thermal diffusion equation, a new numerical method is proposed in this paper as follows. In many cases, it is advantageous to separate the variables t and z in T(z, t). In general, the Fourier series expansion T(z, t) in the frequency domain can be expressed as [33]: where j is an imaginary unit, A 0 is a real field variable, and A 0 /2 represents the DC component of the sound wave. When A n (n ≥ 1) is the complex field variable (i.e., complex amplitude) of the nth−order harmonic, the real part of A n exp(jnωt) is the real displacement of the nth-order harmonic, and A * n (n ≥ 1) is the complex conjugate field variable of A n . Note that A n and A * n no longer contain the time variable t, which is a function of the spatial coordinate z.
Generally, higher−order harmonics are weaker, and some higher−order harmonics can be ignored according to the requirements for calculation accuracy. To simplify the theoretical description, only harmonics of the order less than or equal to N(N ≤ 5) are considered in this paper, and other higher−order harmonics are ignored; this is referred to as an N − order approximation.
Substituting Equation (6) into Equation (5), the following equation can be obtained by orthogonality: where, Notably, Equation (7) gives only an equation of field variables and not an equation of conjugate field variables. The equation of conjugate field variables can be obtained by taking the complex conjugate of Equation (7); therefore, Equation (7) is complete. Substituting i = 0, ±1, ±2, ±3, ±4, ±5 into Equation (7) yields a set of nonlinear equations. However, it is difficult to solve these equations directly. Therefore, a simple iterative method for solving these equation is proposed, where A * (m) and A (m) denotes the field quantities obtained from the m − th iterative calculation. In the iterative calculation, the following method is used. First, the field quantities on the left side of Equation (7) are taken as A . Therefore, the following equation is used in the m − th iterative calculation: when i = 0, ±1, ±2, ±3, ±4, ±5, an uncoupled set of equations can be obtained from (14).
can be calculated separately and independently from the other iterations, which means that the computational effort does not increase dramatically when higherorder harmonics are involved.
Outputting results A i (i = 1) of the 2 iteration calculation are same with those predicted by the perturbation theory. The depletion of pump waves has already been taken into account in the 2 interaction calculation. However, it is not taken into account in the perturbation theory [31].
In the iterative calculation, this study let A . At both endpoints z = 0 and z = d, this work considers that the photothermal radiation signal is mainly due to the alternating temperature field. To simplify the calculation, this paper ignores the DC term and the effect of other layers and convective radiation on the temperature

Numerical Results and Discussion
In this section, we numerically calculate the solution to Equation (5) using the new numerical iterative method and the conventional perturbation method and discuss the effect of different sample parameters on the amplitude of the posterior surface of the sample. Table 1 lists the physical parameters used in the calculation of the PWS samples. The thermal conductivity is K 0 = 121.5 W/m and the thermal diffusion coefficient is α b = 5.9419 × 10 −5 m 2 /s for the 2219 aluminum alloy sample.
The amplitude of the signal from the posterior surface of PWS samples of different thicknesses obtained by the two numerical methods decreases with increasing frequency in the low−frequency and high−frequency ranges as shown in Figures 2 and 3, respectively. In the low−frequency range, Figure 2 shows that the amplitude of each order harmonic of the signal from the posterior surface of the PWS samples obtained by both numerical methods decreases with increasing frequency when other parameters are constant and the sample thicknesses are d = 0.5 mm and d = 0.8 mm. In the high−frequency range,   As seen in Figures 2 and 3, the results obtained by the new numeri method are more sensitive to thickness than those obtained by the convention tion method. Figure 2 shows the results in the low−frequency range PWS sa A1(normalize) A2(normalize) As seen in Figures 2 and 3, the results obtained by the new numerical iterative method are more sensitive to thickness than those obtained by the conventional perturbation method. Figure 2 shows the results in the low−frequency range PWS samples with thicknesses d = 0.5 mm and d = 0.8 mm. Both the fundamental and second harmonics on the posterior surface of the samples obtained by the two numerical methods decrease with increasing thickness, which is in good agreement with the results predicted in the literature [28,29]. For the fundamental and second harmonics obtained by the conventional perturbation method, there is no significant change in the effect of the sample thickness on the results, and there is no significant difference between the fundamental and second harmonics. However, the fundamental and second harmonics obtained by the new numerical method have a statistically significant effect on the results due to the change in the sample thickness, and the change is more pronounced in the second harmonic, so that it is necessary to consider the second harmonic. Figure 3 shows the results in the high-frequency range for the PWS samples with the thicknesses of d = 0.01 mm and d = 0.02 mm. Both fundamental frequency waves and second harmonics on the posterior surface of the samples obtained by the new method decrease with increasing thickness. However, both fundamental frequency waves and second harmonics on the posterior surface of the samples obtained by the conventional perturbation method increase with thickness, contradicting the literature predictions. Therefore, it is possible that the conventional perturbation method is not applicable at high frequencies. Additionally, the new numerical method shows a more pronounced change in the second harmonic than in the fundamental frequency wave when a weak change in the sample thickness occurs. Therefore, the second harmonic has stronger sensitivity to thickness, and the new method has the potential for important applications in the noninvasive assessment of PWS thickness.
Reference [31] applied the perturbation method to the solution of the nonlinear thermal diffusion equation to theoretically study the nonlinear photoacoustic effect; while the physical meaning of its simple method and solution is clear, Reference [31] considers only the effect of lower−order harmonics on higher−order harmonics and ignores the inverse effect. Figures 2 and 3 show that the sensitivity of the perturbation method to the sample thickness and the applicable frequency range is slightly lower effective than that of the new method.
In Figure 4, the solid line shows the thickness of the PWS sample inv d = 0.01 mm, and the marked dashed line shows the thickness of the PWS sample in d = 0.02 mm. The numbers 1−5 indicate the fundamental frequency wave and second through fifth harmonics, respectively. Figure 4 shows the variation in the amplitude of each order of harmonics the frequency on the posterior surface of PWS samples at different thicknesses obtained by the new method. From Figure 4, the following conclusion can be drawnwhile other parameters held constant. First, each order of harmonic decreases with increasing frequency. Second, each order of harmonic decreases with increasing thickness. Third, the rate of change with thickness is larger for higher−order harmonics than lower−order harmonics, indicating that the effects of higher−order harmonics are more sensitive to the change in the sample thickness than those of the lower−order harmonics. Therefore, it is necessary to consider higher−order harmonics. In addition, the difference between the rates of change of the fourth and fifth harmonics with thickness is not very obvious. Therefore, the choice N ≤ 5 is appropriate in the theoretical derivation. Figure 5 shows the variation of each order of harmonics with light energy on the posterior surface of the PWS samples at different thicknesses obtained with the new method. Figure 5 shows that the amplitude of each order of harmonics increases with increasing light energy when other parameters are constant. Additionally, when other parameters are constant, the amplitude of each order of harmonic decreases with increasing thickness. As observed from Figure 5a, the fundamental frequency amplitude is proportional to the optical energy w p , which is consistent with the results of linear theory. As observed from    Figure 4, the following conclusion can be drawnwhile ot ters held constant. First, each order of harmonic decreases with increasing freq ond, each order of harmonic decreases with increasing thickness. Third, the ra with thickness is larger for higher−order harmonics than lower−order harmon ing that the effects of higher−order harmonics are more sensitive to the change ple thickness than those of the lower−order harmonics. Therefore, it is neces sider higher−order harmonics. In addition, the difference between the rates o the fourth and fifth harmonics with thickness is not very obvious. Therefore 5 N ≤ is appropriate in the theoretical derivation. Figure 5 shows the variation of each order of harmonics with light en posterior surface of the PWS samples at different thicknesses obtained w method. Figure 5 shows that the amplitude of each order of harmonics inc increasing light energy when other parameters are constant. Additionally, wh rameters are constant, the amplitude of each order of harmonic decreases wit thickness. As observed from Figure 5a, the fundamental frequency amplitud tional to the optical energy p w , which is consistent with the results of linea observed from Figure 5b,c, the amplitudes of the higher−order harmonics are p to the square of the light energy, which is consistent with the theoretical deriv To demonstrate the validity of the new theory, two numerical methods were used to analyze and compare the results for the variation in the amplitude with sample thickness on the posterior surface of 2219 aluminum alloy samples. The two methods are the new iterative numerical method and the conventional perturbation method [28,29].
The effectiveness of the new iterative numerical method for solving the nonlinear heat diffusion equation is demonstrated in Figures 6 and 7. Numerical simulations of the nonlinear heat diffusion equation for the 2219 aluminum alloy sample using the new numerical iterative method show two effects. First, the amplitude of each order of harmonics decreases with increasing frequency. Second, the second−order harmonic amplitude decreases with increasing thickness. The above conclusions are consistent with those obtained by the conventional perturbation method, and Figure 6 shows that the fundamental frequency wave results obtained by the two methods are very consistent. Therefore, the effectiveness of the method is demonstrated in Figures 6 and 7. (a) Fundamental frequency wave; (b) second harmonic; (c) third harmonic; (d) fourth harm To demonstrate the validity of the new theory, two numerical methods were analyze and compare the results for the variation in the amplitude with sample th on the posterior surface of 2219 aluminum alloy samples. The two methods are t iterative numerical method and the conventional perturbation method [28,29].
The effectiveness of the new iterative numerical method for solving the no heat diffusion equation is demonstrated in Figures 6 and 7. Numerical simulation nonlinear heat diffusion equation for the 2219 aluminum alloy sample using the n merical iterative method show two effects. First, the amplitude of each order of har decreases with increasing frequency. Second, the second−order harmonic amplitu creases with increasing thickness. The above conclusions are consistent with th tained by the conventional perturbation method, and Figure 6 shows that the funda frequency wave results obtained by the two methods are very consistent. Therefo effectiveness of the method is demonstrated in Figures 6 and 7. An examination of Figures 6 and 7 shows the superiority of the new numerical method. For the 2219 aluminum alloy sample, the rate of change with the frequency of the second harmonic is larger than that of the fundamental frequency wave, and the rate of change of the results obtained by the new numerical method is larger than that of the second harmonic results obtained by the traditional perturbation method. Additionally, Figure 7 shows that the rate of change of the second harmonic amplitude with increasing thickness with other parameters unchanged is greater with the new method than with the traditional method. Therefore, the new numerical method is better than the traditional perturbation method. This comparison reflects the greater sensitivity of the results to the sample thickness obtained by the proposed method compared tothe traditional method. Figure 8 describes the relationship between the second harmonic amplitude on the posterior surface of the 2219 aluminum alloy sample and the iteration frequency m when the light energy w p differs. Figure 8 shows that when the light energy is 5 mJ/cm 2 and the number of iterations is greater than or equal to 6, the fundamental frequency amplitude converges to 7.6485 K. In addition, when the light energy is 8 mJ/cm 2 and the number of iterations is greater than or equal to 10, the fundamental frequency amplitude converges to 19.8287 K. Moreover, when the light energy is 10 mJ/cm 2 and the number of iterations is greater than or equal to 12, the fundamental frequency amplitude converges to 31.3493 K. These results imply that the proposed method is converges well.

Conclusions
In this paper, a theoretical model for thermal nonlinear photoacoustic detec lated to port-wine stain samples is constructed. A new numerical iteration app developed and compared to the perturbation method for analyzing the thickness sorption coefficient of PWS samples. The main conclusions are as follows: (1) The rates of change with frequency, thickness, and optical energy intensity ar for higher−order harmonics than lower-order harmonics; higher−order ha are more sensitive to sample detection than lower-order harmonics. (2) For the same parameter values, the proposed new numerical iterative met greater sensitivity and a wider frequency band than the perturbation meth thermore, the calculation time of our proposed method will not drastically when additional high−order harmonics are included.

Conclusions
In this paper, a theoretical model for thermal nonlinear photoacoustic detection related to port-wine stain samples is constructed. A new numerical iteration approach is developed and compared to the perturbation method for analyzing the thickness and absorption coefficient of PWS samples. The main conclusions are as follows: (1) The rates of change with frequency, thickness, and optical energy intensity are larger for higher−order harmonics than lower-order harmonics; higher−order harmonics are more sensitive to sample detection than lower-order harmonics. (2) For the same parameter values, the proposed new numerical iterative method has greater sensitivity and a wider frequency band than the perturbation method. Furthermore, the calculation time of our proposed method will not drastically increase when additional high−order harmonics are included.
Author Contributions: N.C. initiated the idea; H.C. advised on the study design and the analysis plan; N.C. facilitated the study design, specifically the integration of disciplines, and wrote the first draft of the paper, which was reviewed and edited by H.L., R.Z. and Y.L. All authors have read and agreed to the published version of the manuscript.

Funding:
The study was supported by the National Natural Science Foundation of China, grant number 11074159, and the Fundamental Research Funds for the Central Universities, grant number 2021TS088.
Institutional Review Board Statement: The study received ethical approval from the Shaanxi Normal University Institutional Review Board, and all procedures performed in studies involving human participants were following the 1964 Helsinki declaration and its later amendments.
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 (H.C.), upon reasonable request.