Semianalytical solution for the transient temperature in a scattering and absorbing slab consisting of three layers heated by a light source

We derived a semianalytical solution for the time-dependent temperature distribution in a three-layered laterally infinite scattering and absorbing slab illuminated by an obliquely incident collimated beam of light. The light propagation was modeled by the low-order \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_1$$\end{document}P1 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_3$$\end{document}P3 approximations to the radiative transfer equation with closed form expressions for eigenvalues and eigenvectors, yielding a quickly computable solution, while the heat conduction was modeled by the Fourier equation. The solution was compared to a numerical solution using a Monte Carlo simulation for the light propagation and an FEM method for the heat conduction. The results showed that using the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_3$$\end{document}P3 solution for the light propagation offers a large advantage in accuracy with only a moderate increase in calculation time compared to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_1$$\end{document}P1 solution. Also, while the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_3$$\end{document}P3 solution is not a very good approximation for the spatially resolved absorbance itself, its application as a source term for the heat conduction equation does yield a very good approximation for the time-dependent temperature.

Non-contact thermal property measurement methods are widely used because they avoid contact resistances and require little to no sample preparation 1 . While there is a large variety of detection principles, most techniques use optical heating in the visual and near-infrared range. Specifically in photothermal radiometry (PTR), the detection is based on recording the thermal emission using an IR sensor. In the recent past, PTR has e.g. been used to retrieve microhardness profiles in case hardened steel 2 , detect partial curing in dental resins 3 or study drug diffusion in human skin 4 . For semi-transparent, non-scattering media, Ravi et al. studied the reconstruction of absorption profiles using modulated PTR 5 , while Salazar et al. used the same method to simultaneously retrieve the absorption and thermal diffusivity profiles in semi-infinite media 6 and multi-layered slabs 7 , utilizing analytical solutions of the heat conduction equation. Recently, Ren et al. extended this to temperature-dependent medium parameters for a 1D planar symmetric medium with coupled radiation and conduction using numerical solutions 8 , but also specifically excluded scattering. Apart from parameter retrieval with PTR, the required forward solutions can also be used to predict temperature distributions and, especially in biomedical applications, thermal damage after laser irradiation 9 . Here, as in many other media, light scattering is a large effect and therefore must be accounted for in a model for optical heating. The radiative transfer equation (RTE) is generally considered a very accurate model for light propagation on mesoscopic and macroscopic scales for scattering and absorbing media 10 . But due to its complexity, the RTE is often replaced by the diffusion equation (DE) 11,12 , which may be derived as an approximation to the RTE. Solutions to the DE are usually simple and can be computed quickly, but they are known to be inaccurate in many cases 13 . The P N approximations to the RTE 14 on the other hand offer a much higher accuracy, but the computation time increases rapidly with the approximation order N. For low approximation orders N however, quickly computable closed form expressions can still be found 15 . While these low approximation order solutions can still exhibit large errors for the spatially resolved fluence or absorbance rate, their low spatial frequency components are already quite accurate. Using them as a source term for the heat conduction equation should therefore already produce highly accurate solutions for the medium temperature. The aim of this work therefore is to derive a semianalytical solution for the time-dependent temperature distribution in a three-layered scattering and absorbing medium heated by an incident beam of light, where the www.nature.com/scientificreports/ light propagation is modeled using low order P N approximation to the RTE. Specifically, we use the expansion orders N = 1 , which leads to a diffusion-like approximation, and N = 3 , for which closed form expressions are available for the homogeneous solution, resulting in short computation times.
In "Solution of the RTE" section, we present our light propagation model using the P N approximation to the RTE, including the used source and the optical boundary conditions for the three-layered medium. From this, we derive the source term for the heat conduction model, which is described in "Solution of the heat equation" section. Finally, in "Results and discussion" section, we compare the results of our solution to a purely numerical solution for a simple model of human skin.

Solution of the RTE
General solution. The steady-state radiative transport equation (RTE) governing the radiance I(r,ŝ) at position r and direction ŝ for a point source at r = 0 radiating in direction ŝ 0 is given by with the phase function where µ a is the absorption coefficient, µ s is the scattering coefficient, µ t = µ a + µ s is the attenuation coefficient, f l are the phase function moments and P l the Legendre polynomials. In order to improve the accuracy of low order P N approximations, we first approximate the phase function using the delta-M method 16,17 of order N as Next, we split a ballistic part off the radiance as I = I 0 + I d . The ballistic part I 0 is then given by 18 where the modified extinction coefficient becomes μ t = µ a + (1 − f N+1 )µ s . The diffuse radiance I d then satisfies the modified RTE with the new source term This new source corresponds to a line source in direction ŝ 0 with exponentially decreasing fluence, radiating with the modified phase function rotated in beam direction. In the P N method, the diffuse radiance is expanded in spherical harmonics Y lm with expansion coefficients lm as Using the transformed quantity the P N equations are then given by 15 for 0 ≤ m ≤ N and m ≤ l ≤ N with the condition where the right hand side E lm must be determined from the source term (see Eq. (15)). The coefficients are given by 15 (1) s · ∇I(r,ŝ) + µ t I(r,ŝ) − µ s www.nature.com/scientificreports/ and Note that the rotation e −imφ q in Eq. (8) makes the left hand side of the system (9) independent of φ q , allowing a much more efficient numerical evaluation. The σ l in Eq. (11) are the coefficients that follow directly from applying the P N method to Eq. (1) without the delta-M method and the ballistic part separation 14 . It is therefore evident that the homogeneous solutions of Eqs. (1) and (5) are identical. These solutions can for example be found using the method of rotated reference frames, which was studied intensively over the last few years. The solution for the diffuse radiance then has the form 14 Here, C i and C ′ i are unknown coefficients determined by the boundary conditions, q 2 + 1 ξ 2 i =: ζ i (q) are the positive eigenvalues and � i lm (q) the corresponding eigenvector components. We refer the reader to previous publications 14,15,[19][20][21] for details on this method.
For the particular solution � (p) lm (q, φ q , z) with the source term (6), we first calculate its spherical harmonics decomposition Applying a 2D Fourier transform, setting ŝ 0 = (θ 0 , 0) in spherical coordinates and making use of condition (10), we obtain the desired moments with µ 0 = cos(θ 0 ) . Introducing the abbreviation 22 the similarity ansatz � (p) lm = χ lm e −µ c z can be used to solve the system (9) with the source moments (15), resulting in the system of equations for z ≥ 0 and 1 ≥ µ 0 > 0 In many cases, only the fluence rate �(q, φ q , z) = S 2 I(q, φ q , z,ŝ)d 2 s is required instead of the radiance from Eqs. (4) and (13). Using the same set of polar coordinates for the ballistic part, the fluence rate is found to be RTE boundary conditions for a three-layered slab. Our goal here is to construct a solution for a stack of three laterally infinite layers, where for simplicity all layers share the same refractive index. The refractive index outside the stack may, however, be different from the one inside so that Fresnel reflection at the top and bottom surfaces must be accounted for in the boundary conditions. Figure 1 illustrates this problem geometry schematically. www.nature.com/scientificreports/ To that end, we first rescale the coefficients C i and C ′ i from equation (13) for each layer n ∈ {1, 2, 3} with thickness L (n) and take the extinction and shift D(n) of the incident light by the layers above into account, resulting in where we defined and Using µ = cos(θ) and µ > 0 , the boundary conditions are given by 13 where R(µ) is the Fresnel reflection coefficient for the top and bottom surfaces 23 . Multiplying all conditions with Y * l ′ m ′ (µ, ϕ) and integrating over the half-space µ > 0 then yields the generalized Marshak boundary conditions 14,22,24   i , completing the solution for the P N approximation to the RTE. In all of the above, a single source with unit strength just below the upper boundary was assumed. To model an incident collimated beam from outside the medium, we have to take the Fresnel reflection outside the medium as well as internal reflections of the unscattered light into account. If R d (µ 1 ) is the reflection coefficient outside the medium and R 0 = R(µ 0 ) the reflection coefficient inside the medium, where µ 1 and µ 0 are connected by Snell's law, the sources due to multiple internal reflections can be reduced to two sources with strengths S o and S u for the top and bottom surface, respectively, with In many cases, R 0 D(4) is very small and the lower source may be neglected. If however the medium is optically thin enough for the lower source to matter, the coefficients C (n) i and C ′(n) i must be recalculated with reversed layer order for a complete solution.

Solution of the heat equation
We assume that the heat transport inside the medium is governed by the differential equation for the temperature T(r, t) , where ρ is the density, c p the heat capacity and k the isotropic thermal conductivity. n is the layer number and m the layer containing a source sheet at depth z 0 with spatial profile Q s (x, y) and time dependence Q t (t)�(t) . For the heat equation (35) to be valid, the spectral components of Q t (t) must be restricted to frequencies well below 1 GHz 25 . This allows using the steady-state RTE absorbance for the spatial source profile. Also, the layers should be at least several micrometers thick. Applying the same 2D-Fourier transform we used for the RTE and a Laplace transform with respect to time to (35) using T = 0 everywhere as initial condition, we get with the Laplace variable s www.nature.com/scientificreports/ where the tilde marks transformed quantities. We note that in our previous work 26 , we derived a more general solution of the heat conduction equation for N layers with anisotropic thermal conductivities. The solution we require here, therefore, arises as a special case and is given by The coefficients A 1 for a three-layered system and sources in the first layer. Since in our case there are sources in all layers, we also need the expressions for the remaining 12 coefficients. For sources in the third layer, they follow immediately from the coefficients for sources in the first layer with reversed layer properties. The coefficients A m (s, q, φ q , z, z 0 ) =A (n) m e α (n) (z− (n) −L (n) ) + B (n) m e −α (n) (z− (n) ) + δ nmP (m) (s, q, φ q , z, z 0 ),

Results and discussion
In this section, the obtained solution is tested and validated against a purely numerical solution. As an exemplary system, we choose a three-layered generic model of human tissue in the near-infrared range 27,28 . The three layers represent skin, subcutaneous fat and muscle tissue, respectively. The thermal parameters of the model are listed in Table 1 and the optical parameters in Table 2. For all layers, we assume a refractive index of n i = 1.4 with n o = 1.0 outside the medium and the Henyey-Greenstein phase function with g = 0.8 using the moments f l = g l . For the upper boundary, we use a heat transfer coefficient of h 1 = 100 W m −2 K −1 with an ambient temperature of T a = 0 equal to the initial medium temperature and assume an adiabatic lower boundary with h 2 = 0.
As an input beam, we use a Gaussian beam with beam radius r w = 1 mm hitting the top surface with an angle outside the medium of µ 1 = cos(θ 1 ) , centered at the origin. We use θ 1 = 0 • and θ 1 = 60 • as incidence angles. Taking the distortion of the beam profile due to the oblique incidence into account, the input beam profile is then given by 22 Applying a 2D Fourier transform to (49) again yields .   Table 2. Optical parameters of the considered medium. For the phase function, the Henyey-Greenstein function is used. The refractive index n i must be identical for all layers, while the refractive index outside the medium is set to 1.0. www.nature.com/scientificreports/ The lower source strength from Eq. (34) depends on the approximation order and on the angle of incidence, but even for normal incidence with the P 1 approximation, it is nine orders of magnitude weaker than the upper source. Therefore, we can safely neglect the lower source in all calculations. As approximation orders, we consider only the P 1 and P 3 approximations, because for these, there are closed form expressions available for the eigenvalues and eigenvectors. For the P 3 , we use the expressions given by Liemert et al. 15 , while the P 1 expressions are given in Appendix B. Although theoretically possible up to P 7 , no higher order closed form expressions are available to the authors' knowledge. Computing the eigenvalues and eigenvectors numerically is possible for higher orders, but computationally quite expensive. Also, the solution would become more and more unstable for high values of q, which would render the numerical inverse Fourier transform extremely tricky. To obtain a numerical solution to compare to, we use a two step process. First, the light propagation is calculated using a custom GPU-accelerated Monte Carlo solver, recording the spatially resolved fluence rate inside the medium. Then, we use the COMSOL Multiphysics software 29 to interpolate the resulting absorbance from the Monte Carlo simulation and calculate the time-dependent temperature rise. Since we expect the P N approximation to be the main source of error in the final results, we first look at the predicted fluence rate. Figure 2 contains the depth-resolved total fluence rate calculated in P 1 and P 3 approximation, integrated over x and y. Here, we do not yet have to resort to Monte Carlo methods to obtain a correct reference, since for planar symmetry or q = 0 , the P N equations can be solved efficiently for high orders N. Integrating the depth resolved total fluence rate over a layer and multiplying by the respective µ a yields the total absorbance inside that layer. As the dependence on z of the fluence rate is very simple (compare Eqs. (18) and (19)), this integration can be done analytically and these total absorbances per layer can, therefore, be calculated with little computational effort. Table 3 contains the errors of the total absorbances per layer for the P 1 and P 3 approximations relative to the practically exact P 201 . It can be seen that the P 3 is clearly superior to the simpler P 1 approximation. Especially in the vicinity of the light source near the upper boundary, the P 1 introduces large errors. However, the P 3 still underestimates the total fluence and therefore the total heat source strength per layer by 0.7 − 0.9% . Being able to quickly quantify these errors, in addition, enables us to correct them. This results in a weighting factor for the heat source strength in each layer, making sure that the total source strength in each layer, and consequently also in the whole medium, is correct. The corresponding fluence rates are also shown in Fig. 2.
We implemented our solution in Python. The implementation includes the planar symmetric P N approximation solution, the 3D P 1 and P 3 approximation solutions, the solution of the heat equation with the P 1 and  Table 3. Total absorbance error for each layer of the P 1 and P 3 approximations relative to a P 201 approximation that may be regarded as an exact solution. www.nature.com/scientificreports/ P 3 absorbance sources and the required numerical transform algorithms and is freely available together with the reference results of the numerical calculations 30 . The transform algorithms are basically the same that were used in a previous work 26 . For the comparison, we calculate the results along the line y = 0 mm at two depths z = 0 mm and z = 4 mm , t = 1 s and t = 15 s after the source was switched on. Figures 3 and 4 show the results for normal incidence, while the results for an incident angle of θ 1 = 60 • are shown in Figs. 5 and 6. For each curve, 400 points were calculated. The performance depends, of course, on the efficiency of the numerical transforms and on the required accuracy. For the present calculations, we used 39 points for the inverse Laplace transform and 240 × 80 points for the inverse 2D Fourier transform, which resulted in a computation time of 3.4 s for the P 3 approximation and 2.7 s for the P 1 approximation per curve on a standard desktop PC. For the numerical solution, the Monte Carlo simulation took 16 min and the subsequent finite element solution took 44 min to complete with reasonable accuracy. As can be seen, the P 3 solution agrees quite well with the numerical result, while the P 1 predictably shows larger errors up to 7% . For short times, the errors are slightly larger for all cases. This is to be expected, since the errors in the P N absorbance are mostly contained in the high spatial frequency components. These are most prominent in the temperature solution for short times and then decrease due to the effects of heat conduction. To confirm that the remaining error of the P 3 approximation is indeed due to a slight misprediction of the heat source geometry instead of the interpolation or the finite element calculation, we repeated the numerical simulation with absorbance values computed from the P 3 approximation instead of the Monte Carlo simulation. For this comparison, we observed a practically exact agreement of the solutions.

Conclusions
In conclusion, we derived a semianalytical solution for the heating of a scattering and absorbing three-layered medium by an incident beam of light, where the light transport is modeled using the P 1 and P 3 approximations to the RTE. We compared our solution to a purely numerical one and observed a good agreement in the case of the P 3 approximation for the light transport. The simpler P 1 approximation showed significantly larger errors and should not be used, given the modest advantage in calculation time compared to the P 3 . For the comparison above, the semianalytical solution was several orders of magnitude faster than the numerical solution. An entirely fair comparison is of course difficult to perform, since the numerical simulation unavoidably also yields the absorbance for all positions inside the medium and the temperature for all positions and times. But in most situations, only a small subset of this data is actually required and this is where the semianalytical solution offers large improvements. Finally, some simple generalizations are possible for the heat conduction, since we only require the form of (36). For example, the equation for a moving heat source has the same form with a slightly different α where u x and u y are the velocity components of the source in the x-y-plane.

A Coefficients A (n) 2 and B (n) 2 for sources in the second layer
Here, we give closed form expressions for the coefficients A 3 follow immediately by reversing the layer order. Adopting the same notation, the required coefficients follow from the system of equations (51) α (n) u = ρ (n) c (n) p k (n) s + iq(u x cos φ q + u y sin φ q ) + q 2 , Figure 5. Temperature after t = 15 s with an incidence angle of θ 1 = 60 • along the line y = 0 for two different depths. While the P 3 is still clearly superior to the P 1 approximation, neither of the approximations deteriorates much with oblique incidence.