Numerical Simulation of Enhancement of Superﬁcial Tumor Laser Hyperthermia with Silicon Nanoparticles

: Biodegradable and low-toxic silicon nanoparticles (SiNPs) have potential in different biomedical applications. Previous experimental studies revealed the efﬁciency of some types of SiNPs in tumor hyperthermia. To analyse the feasibility of employing SiNPs produced by the laser ablation of silicon nanowire arrays in water and ethanol as agents for laser tumor hyperthermia, we numerically simulated effects of heating a millimeter-size nodal basal-cell carcinoma with embedded nanoparticles by continuous-wave laser radiation at 633 nm. Based on scanning electron microscopy data for the synthesized SiNPs size distributions, we used Mie theory to calculate their optical properties and carried out Monte Carlo simulations of light absorption inside the tumor, with and without the embedded nanoparticles, followed by an evaluation of local temperature increase based on the bioheat transfer equation. Given the same mass concentration, SiNPs obtained by the laser ablation of silicon nanowires in ethanol (eSiNPs) are characterized by smaller absorption and scattering coefﬁcients compared to those synthesized in water (wSiNPs). In contrast, wSiNPs embedded in the tumor provide a lower overall temperature increase than eSiNPs due to the effect of shielding the laser irradiation by the highly absorbing wSiNPs-containing region at the top of the tumor. Effective tumor hyperthermia (temperature increase above 42 ◦ C) can be performed with eSiNPs at nanoparticle mass concentrations of 3 mg/mL and higher, provided that the neighboring healthy tissues remain underheated at the applied irradiation power. The use of a laser beam with the diameter ﬁtting the size of the tumor allows to obtain a higher temperature contrast between the tumor and surrounding normal tissues compared to the case when the beam diameter exceeds the tumor size at the comparable power.


Introduction
Hyperthermia is a treatment modality based on the local heating of pathological tissue volumes to temperatures exceeding 42 • C. The heat effect produced by these temperatures can damage or kill pathological cells or make them more vulnerable to other therapy modalities, such as radiation therapy or chemotherapy, with minimal damage to surrounding normal tissues [1,2]. In practice, hyperthermia protocols aim at concentrating the heating primarily in the tumor area [3,4]. A promising direction to improve hyperthermia 2-5 nm agglomerated into big particles of sizes ranging from 40 to 300 nm [33]. SiNPs of diameter less than 220 nm were reported to be formed with the help of fracturing porous silicon by ultrasonication in water [37] and ethanol [38] followed by membrane filtration.
The laser ablation of silicon or silicon nanostructures in liquids opens up a promising way to manufacture almost ready-to-use suspensions of SiNPs of the desired size free of any chemical impurities [40][41][42]. For example, SiNPs of mean size of 20-30 nm formed by the laser ablation of porous silicon were employed as agents for photohyperthermia: irradiation of the SiNP aqueous suspension (0.4 mg/mL) containing unicellular Paramecium caudatum organisms by continuous-wave laser radiation with the intensity of 25 kW/cm 2 at a wavelength of 808 nm was found to result in temperature rise up to 20 • C and death of the living cells [43].
Employing arrays of silicon nanowires (SiNWs) as ablation targets seems to have various advantages for producing the SiNPs for biophotonics applications, including satisfactory stability of water-and ethanol-based suspensions with an absolute value of the zeta potential up to 30 mV [44]. This value indicates the typical time scale of spontaneous agglomeration as tens of hours, and the sonification of the SiNPs suspension prior to application allows to exclude the effect of agglomeration in the course of the procedure. Stability can be potentially increased via coating with biocompatible and biodegradable polymers as it has been done for porous silicon [45].
SiNW arrays fabricated via a metal-assisted chemical etching technique [46] contain SiNWs of about 100 nm in diameter and up to 100 µm in length which demonstrate strong scattering and high absorption of visible and near-infrared light [47,48]. Ablation threshold of SiNWs is up to eight times lower than in crystalline silicon [34,35], which increases product yield and reduces the estimated production cost. However, to our knowledge, the laser-ablated SiNWs have not been reported to be used in photohyperthermia.
In this work, we aimed at estimating the potential of SiNPs fabricated by the laser ablation of the SiNW arrays in water and ethanol as thermal coupling agents for hyperthermia. We performed a numerical simulation of heating a millimeter-sized nodular basal cell carcinoma (nBCC), which is the most common form of human skin cancer, located in human facial skin tissue [49] using continuous-wave (cw) laser irradiation. The choice of cw laser radiation for hyperthermia is determined by limits of laser radiation intensity for safe laser treatment of biotissues (below 200 mW/cm 2 [50]), especially for healthy tissues surrounding the tumor. Employing pulsed irradiation has no advantages due to the decrease of safety energy density limits with shortening pulse duration [50,51]. Although silicon possesses significant cubic nonlinear susceptibility responsible for the light self-action [52], simple estimations demonstrate that for a laser pulse energy of the order of safety threshold [51], the relative variation of absorption for bulk silicon is in the order of 10 −4 .
We considered two configurations of the beam: (1) beam width corresponds to tumor transversal size to avoid excessive surrounding tumor heating, and (2) beam width exceeds it to ensure tumor heating at transversal boundaries. The considered beam power range was 60-200 mW. A systematic analysis of the effect of SiNPs concentration and beam configuration of hyperthermia performance was conducted. Optimal irradiation parameters providing hyperthermia of the entire tumor without significant overheating of surrounding healthy tissues were found.

Materials and Methods
The developed numerical model of laser-induced tumor hyperthermia with the embedded SiNPs consists of three parts. In the first step, we calculated the partial scattering and absorption coefficients of the nanoparticles embedded in biotissue using Mie theory based on the measured size distributions of the two types of employed nanoparticles. In the second step, we simulated light absorption distribution in a biotissue model mimicking skin with a superficial nBCC upon laser irradiation with a particular beam configuration using optical parameters accounting for the SiNPs presence in the tumor. In the third step, a bioheat transfer equation was employed to evaluate the temperature distribution in the tumor and the surrounding tissues. Several configurations of laser beam and power, as well as various SiNPs concentrations, have been considered in order to analyse hyperthermia efficiency and ensure the heating of healthy tissues below the damaging threshold. Figure 1 shows histology images of normal skin ( Figure 1a) and typical nBCC (Figure 1b). The histology images were obtained from the institutional database and were fully anonymized. In the case of nBCC, stratum corneum and upper epidermis layers have normal morphology, while the tumor node has a shape close to the ellipsoidal one. To reproduce the observed morphological structure in the simulations, we modeled human facial skin by a three-layered medium consisting of epidermis, dermis, and subcutaneous fat. The nBCC is represented by a spheroid located at a depth of 0.02 mm below the skin surface, with semiaxes of 2.5 mm in the two directions parallel to the skin surface and 1 mm in the direction normal to skin surface, which is in agreement with the morphological observations.  Figure 2 demonstrates the layout of the nBCC hyperthermia scheme with the position of the laser beam at the skin surface. We considered cw laser irradiation with the circular homogeneous intensity distribution. The radius of the circle is r, and the laser beam axis coincides with the axial symmetry axis of the tumor. Laser radiation is assumed to be non-polarized. This beam shape is typical for hyperthermia carried out by light-emitting diodes [53] or lasers generating in multimode regime [54]. In these cases, circular homogeneous intensity distribution can be achieved with the help of laser radiation homogenizers, enabling a spatially flat-top laser spot of non-polarized laser radiation (see, e.g., [55]). Applying a homogeneous beam seems to be more practical due to the more homogeneous transversal distribution of the intensity as compared to alternatively employed Gaussian beams.

Geometry and Optical Parameters of the Model
To ensure deep penetration of the laser radiation into the treated biotissue region, a wavelength from the optical transparency window (600-1300 nm) should be used. At the same time, in this wavelength range, the absorption of light by silicon decreases with the increase in wavelength [56]. Our previous study [39] compared the application of the wavelengths of 633 and 800 nm for nanoparticle-mediated tumor hyperthermia. It was shown that at equal intensities the wavelength of 633 nm provides higher local heating, which determined the choice of this wavelength for this study. Typical thicknesses of the simulated biotissue layers and their optical parameters at 633 nm are shown in Table 1 in accordance with Refs. [57][58][59].

Optical Characteristics of SiNPs Employed in the Simulations
In simulations, we considered the SiNPS obtained via the ablation of SiNW arrays by picosecond laser at wavelength of 1064 nm in water and ethanol (wSiNPs and eSiNPs, correspondingly) and reported in our earlier work [44]. According to previous studies of similar laser-ablated SiNPs, they may have a thin oxide shell [36,60] as well as a small fraction of amorphous phase [61]. Besides, eSiNPs undergo surface carbonization due to high power laser action to ethanol [36]. However, these effects may be excluded from the following consideration in our calculations due to their negligible contribution to SiNPs optical properties, and pure crystalline silicon particles are considered in the study.
Based on scanning electron microscopy data, a normalized size distribution function f i = n i N was obtained for wSiNPs ( Figure 3a) and eSiNPs ( Figure 3b) where index value i corresponds to the i-th size fraction with step of 4 nm, n i is the number of particles in a unit of volume (partial volume concentration) for i-th size fraction, and N = ∑ i n i is the number of all SiNPs in the unit of volume (total volume concentration). In practice, the number of nanoparticles in the tumor is characterized by the mass concentration C m (in mg/mL), given by the relation: where ρ is the silicon density (ρ = 2.33 g/cm 3 ) [62] and V i is the average volume of a single particle with the diameter corresponding to i-th size fraction. Therefore, the concentration n i of particles with diameter d i can be expressed via the known size distribution function f i and mass concentration C m as follows: Assuming the spherical shapes and high crystallinity of both synthesized wSiNPs and eSiNPs, Mie theory [63] can be employed to calculate scattering cross-section σ s,i , absorption cross-section σ a,i (Figure 3c), and anisotropy factor g i for i-th size fraction of silicon nanocrystals. Due to the small SiNP volume fraction, which does not exceed 0.003 for the maximal SiNP concentration, scattering (µ tum+NP s ) and absorption (µ tum+NP a ) coefficients of the tumor with the embedded SiNPs can be found by adding the scattering and absorption coefficient of all nanoparticles (∑ i σ s,i n i and ∑ i σ a,i n i , respectively) to scattering (µ tum s ) and absorption (µ tum a ) coefficients of the tumor without SiNPs, respectively [19]: The relation between scattering anisotropy factor g tum+NP for the tumor with the SiNPs, scattering anisotropy factor g tum of the tumor without SiNPs, and scattering anisotropy factor g i of i-th SiNP fraction is given by the expression [19]: and the effective attenuation coefficient of diffusive light (µ tum+NP eff ) for the tumor with the embedded SiNPs [64] is: Given that the considered SiNP size exceeds 8 nm, whereas quantum-confinement effects in silicon resulting in variations of band gap and optical parameters could be detected for SiNP diameter below 5 nm [65], to calculate the SiNP optical parameters we used refractive index n NP = 3.8823 and extinction coefficient k NP = 0.0196 for crystalline silicon [66]. Relative variations of these values with temperature increase up to 10 • C for bulk silicon are 10 −3 and 10 −5 , respectively [67,68]. Hence, the thermooptical effect on the SiNPs embedded into the tumor is expected to be negligible. Thermooptical effects in the biotissue are also extremely small. According to Ref. [69], the value of dµ a /dT for biotissues does not exceed 0.02 cm −1 K −1 . Considering the typical temperature increase in the treated biotissues and their absorption coefficients (see Tables 1 and 2), the possible relative variation in absorption is estimated to not exceed 7% for tumor with embedded 1 mg/mL of eSiNPs and an order of magnitude smaller for tumor with 7 mg/mL wSiNPs. The refractive index for tumor tissue at the wavelength of 633 nm is assumed to be equal to n tum = 1.4 [57,70,71]. Optical properties of skin layers and tumor tissue are presented in Table 1. It is worth nothing that all optical parameters values were measured for a non-polarized light source [57].
We expect that the mass concentration of SiNPs required for hyperthermia should be rather high due to the comparatively low absorption of red light. This excludes any possibility of intravenous injection due to low accumulation efficiency [72,73]. Meanwhile, intratumoral injection allows up to 15% of all injected NPs to be concentrated inside the tumor, and NPs can be contained in a tumor in a concentration of up to 9 mg/mL [74]. To reveal the effect of SiNP concentration on the result of hyperthermia, the mass concentration of SiNPs employed in the simulations varied from 1 to 7 mg/mL within the tumor. We expect that during laser hyperthermia the tissue temperature does not significantly exceed 42 • C, and therefore SiNP agglomeration would not occur.
Optical parameters for a tumor containing SiNPs of various types and mass concentrations calculated using formulas (1)-(6) are presented in Table 2. In the calculations, it was assumed that after injection the SiNPs are contained only in the tumor area. The errors for µ a , µ s , and g of the tumor with SiNPs originate from the error of size distribution measurement that were estimated as 6.0% and 10.4% for wSiNPs and eSiNPs, respectively. This yields an estimated error of 2% in scattering and absorption characteristics evaluated by Mie formulas, which caused no relevant change in the obtained results.
According to Tables 1 and 2, light absorption in nBCC without SiNPs is lower than that in healthy skin and human subcutaneous fat. The injection of SiNPs into the tumor at a concentration of 7 mg/mL increases the absorption coefficient of the tumor site 15 times for wSiNPs and 5.5 times for eSiNPs at the wavelength of 633 nm. However, an increase in scattering occurs upon SiNP injection, which yields a decrease of radiation penetration depth in biotissue. The effect of scattering increase is more pronounced for wSiNPs characterized by the maximum diameter of 140 nm compared to eSiNPs characterized by the maximum diameter of 80 nm.

Simulations of Light Absorption Distribution
The Monte Carlo method for light transport in tissues [75][76][77] was employed to calculate the absorption map, i.e., 3D spatial distribution of the power density absorbed by biological tissue upon its irradiation by a continuous-wave laser beam at the wavelength of 633 nm, as shown in Figure 2. Simulations were performed for N = 10 7 photons uniformly distributed within the circle of radius r at the biotissue surface. All photons enter the tissue normal to its surface. The computational grid for the absorption map consisted of 200 × 200 × 500 voxels (the voxel size was ∆x × ∆y × ∆z = 0.1 mm × 0.1 mm × 0.02 mm).
To study the effect of the irradiation area size on tumor heating efficiency, the calculations were performed for two laser beam sizes: r = 2.5 mm and r = 5 mm (Figure 2). In the first case, the laser beam size coincides with the width of the tumor, whereas in the latter case the laser beam covers a larger area, which may simulate hyperthermia on account of the safety margins.
At each step of Monte Carlo trajectory tracing, the absorbed part of the photon packet energy was added to the current value of the absorption map element B(x,y,z) corresponding to the coordinates (x,y,z) of the center of the voxel in which the current trajectory node was located. The resulting absorption map was normalized for the volume of the voxel and for the intensity of probing radiation expressed in photon weight per area unit to calculate the normalized volumetric density of absorbed energy: To derive the discrete dimensional volumetric density of heat generation rate Q ext (x,y,z) within the medium, the normalized volumetric density of absorbed energy should be multiplied by the beam intensity: where P is the laser beam power. The laser power P and the beam radius r were varied in the simulations in order to evaluate the range of optimal parameters to perform SiNPs mediated hyperthermia of nBCC. The range of the irradiation power was limited to ensure healthy biotissue is not heated to temperatures exceeding 42 • C to avoid damaging effects, such as tissue thermal destruction [78,79].

Calculation of Volumetric Temperature Distribution
At the next stage of simulations, the obtained discrete functions of heat sources Q ext (x, y, z) were used to calculate temperature distribution within the biotissue during the laser treatment by means of solving a bioheat transfer equation. The standard way to obtain a temperature distribution within the tissues during a hyperthermia procedure is solving the bioheat transfer equation, which was initially proposed by Pennes in the following form [80,81]: where T is the biotissue temperature, ρ is the biotissue density, C p is the biotissue heat capacity at a constant pressure; k is the thermal conductivity, Q met is the speed of metabolic heat generation the per unit volume, Q ext is the external source power density, which in our case is a result of the laser radiation absorption, and Q perf is power density of the heat loss caused by perfusion (blood transfer through capillaries and extracellular spaces in tissue [82]). Despite great progress achieved during more than 70 years since its first proposal and various improvements and refinements of the equation (see, e.g., [83][84][85]), the Pennes equation is, nevertheless, applicable for the stationary case of continuous-wave irradiation of the biotissue [83]: We solved Equation (10) using a finite-difference approach, taking the absorption map simulated by the Monte Carlo approach as the external source power density. In its stationary form, the Pennes Equation (9) does not depend on biotissue density and heat capacity: Here, k is the thermal conductivity coefficient, Q met = 420 W/cm 3 is the volumetric rate of metabolic heat generation [86], and the last term of Equation (11) describes the perfusion heat loss power density, with ρ bl = 1060 kg/m 3 , C bl = 3770 J/(kg·К) [87,88], T bl = 37.2 • C [89,90], ω bl (Table 3) being the blood density, heat capacity, temperature, and perfusion coefficient (transfer of blood through capillaries and extracellular spaces [82]), respectively. The values of thermal conductivity coefficient k for different biotissues are presented in Table 3. Generally, the thermal parameters of nanoparticles can differ from those for bulk material due to the variation of the phonon spectrum and well-developed surface [91]. However, previously, we have shown that Raman spectra of the employed SiNPs coincide with the spectrum for crystalline silicon [34], while in the current study the volume fraction of the SiNPs does not exceed 0.003 and SiNP heating is expected to occur at several degrees. Based on this, we consider no variations of thermal parameters both for the SiNPs and for the tumor with embedded SiNPs.
In our preliminary studies, we demonstrated that a stationary temperature distribution is achieved in less than 400 s of cw irradiation [39]. It is worth noting that, according to the requirements of the hyperthermia procedure, the irradiation time in cw regime may take up to 60 min [92]. Thus, for therapeutic procedures longer than 400 s, employing a stationary bioheat transfer equation is reasonable.
Equation (11) was solved using COMSOL Multiphysics ® software package using the finite element method [93][94][95]. The computational grid was tetrahedral with the minimum and maximum element sizes of 0.15 mm and 1 mm, respectively. Tetrahedral mesh allows accounting for the complicated shape of the ellipsoidal tumor boundaries when solving the equation. It also significantly reduces the calculation time by varying the grid voxel size depending on the typical scale of the simulated structure (larger voxels for subcutaneous fat and smaller voxels for skin and tumor) compared to the grid with rectangular voxels.
To translate previously calculated distributed heat sources to COMSOL Multiphysics ® environment, we used linear interpolation of the discrete function Q ext (x,y,z) on the tetrahedral grid nodes. Constant temperature of 37.0 • C was set at the rear boundary of the medium located within the biological tissue (z = 10 mm), while at the skin-air interface convection was taken into account [96]: where convective heat transfer coefficient h = 18 W/(m 2 ·K) and air temperature T air = 25.0 • C. Temperature distribution at the lateral boundaries T sides (z) was evaluated from the solution of the 1D stationary bioheat transfer equation in the infinitely wide medium containing no tumor and in the absence of heat sources in biotissue: Boundary conditions for Equation (11) were T| z=10mm = 37.0 • C at the lower boundary (z = 10 mm) and the convection condition at the skin-air interface (z = 0 mm): The values of bioheat transfer equation coefficients employed in the simulation are presented in Table 3. It is worth noting that blood perfusion in a tumor is greater than that in healthy tissue owing to a more developed network of small blood vessels [97,98]. As a result of the performed modeling, thermal fields in tissues containing nBCC tumor without and with embedded SiNPs were calculated. We aim to determine the conditions (nanoparticles concentration and the irradiation parameters) at which the employment of SiNPs allows to heat the tumor over 42 • C completely with the area of overheated healthy tissue being insignificant, which cannot be achieved in the absence of SiNPs in the treated tissue. Figure 4 shows the in-depth profiles of normalized volumetric density of absorbed energy at the beam axis q(x = 0,y = 0,z) in the biotissue with and without different types of SiNPs. The profile for the case without SiNPs shows that the maximum in the absorbed energy is located in the superficial tissue layer above the tumor. This is caused by a lower absorption coefficient of the tumor compared to that of the remaining biotissue together with the in-depth attenuation of the irradiation intensity. The embedding of SiNPs into the tumor leads to a shift of the absorption profile maximum to the tumor region. On the other hand, all profiles corresponding to presence of SiNPs demonstrate lower energy absorption in normal tissue layers underlying the tumor. Since wSiNPs are characterized by a higher absorption coefficient as compared to eSiNPs for the same mass concentration, for them the absorption maximum is located in the upper layer of the tumor, and they shield the penetration of the radiation into deeper tumor regions. For the eSiNPs, the maximum is shifted to the tumor center, and the density of absorbed energy in deeper tumor regions (z ≥ 1 mm) is higher than that for wSiNPs due to weaker shielding. Note that with the increase of the SiNPs concentration the maximum shifts towards tissue surface due to the higher extinction of radiation coming from the top of the tissue. Note that embedding SiNPs into the tumor reduces light absorption in epidermis, since the light is partially absorbed by SiNPs in the tumor, which results in a decrease of backscattering from this area.

The Effect of eSiNPs and wSiNPs on Radiation Absorption Efficiency
Assessing the dependence of function q(x = 0, y = 0, z) on the scattering and absorption parameters of the medium may assist in designing SiNPs with the optimal optical properties to ensure sufficient absorption over the whole tumor region. Note that the in-depth decay rate of q(x = 0, y = 0, z) within the tumor up to the depths of z = 1.5 mm is strongly related to the value of the effective attenuation coefficient µ eff (see Table 2). Moreover, the dependencies q(x = 0, y = 0, z) for eSiNPs with the mass concentration C m = 5 mg/mL and for wSiNPs with the mass concentration C m = 1 mg/mL are characterized by quite similar behavior, while the values of µ eff are very close for these two cases. We found out that the function q fit (z) ∝ exp(−µ eff z), which describes attenuation of light in diffusive regime [75], provides a good fit to the simulated dependencies q(x = 0, y = 0, z) within the tumor region (Figure 4b). The exponential fit is shown for one curve in Figure 4b using the corresponding value of µ eff taken from Table 2. Note that the slope of the function q(x = 0, y = 0, z) increases with the concentration of SiNPs as well as its amplitude defined by the absorption coefficient µ a . Therefore, the effect of the concentration increase on the absorption efficiency is not obvious at first sight: the higher the concentration, the more pronounced are both the absorption and the signal drop. Assuming both the amplitude and the longitudinal scale of the function q(x = 0, y = 0, z), we propose the parameter A for the estimation of the light absorption efficiency within the whole tumor node with the embedded SiNPs that is defined as the product of the absorption coefficient in the tumor over the thickness of the "skin layer" of radiation penetration into the tumor: The dependencies of the absorption efficiency A on the concentration for the two types of SiNPs are presented in Figure 4c. One can see that in the given range of mass concentrations C m , the absorption efficiency almost does not depend on C m for wSiNPs, while for eSiNPs the increase in concentration provides a noticeable rise of the parameter A. However, at very large eSiNPs concentrations (C m > 10 mg/mL) the growth of A(C m ) slows down and obviously tends to saturate. In general, the increase in absorption efficiency obtained by embedding wSiNPs does not exceed 25% at their highest concentration compared to that for the tumor solely, and it is about 40% and 60% for eSiNPs at concentrations of 3 mg/mL and 7 mg/mL, respectively.
The cross-sections q(x, y = 0, z) of the normalized volumetric density of absorbed energy for a beam radius of 5 mm in the absence of SiNPs in nBCC and with SiNPs at concentrations of 3 and 7 mg/mL are shown in Figure 5. As already demonstrated in Figure 4, radiation absorption in nBCC without SiNPs is significantly lower than that of surrounding healthy tissue, which in particular results in discontinuity at the rear border of the tumor (Figure 4a,b). Embedding nanoparticles in the nBCC node significantly enhances absorption efficiency in it, resulting in a higher density of absorbed energy as compared to healthy tissue in superficial tumor region (Figure 5b Figure 4a,b). However, the shielding effect is significantly weaker, resulting in a more uniform in-depth distribution of additional absorption in the tumor, thus confirming higher potential of these particles for hyperthermia.  (Figures 4b and 5c,e), manifested by the deviation from exponential decay for z > 1.5 mm, is caused by the inflow of backscattered diffusive light from the surrounding normal tissues where the fluence level is considerably higher with respect to the tumor due to smaller absorption and scattering. Such an effect is hardly noticed with eSiNPs because of smaller differences in the absorption and scattering coefficients between the tumor with embedded nanoparticles and the surrounding tissues (Figure 5b,d).

The Effect of eSiNPs and wSiNPs on Tumor Heating
The obtained absorption map values were used to calculate Q ext , the volumetric heating rate from inner heat sources, in the bioheat transfer equation. For therapeutic applications, it is important to initiate tissue thermal destruction as a result of irradiation only in the tumor, not in surrounding normal tissues. In this context, the laser power P was varied to obtain temperatures above the threshold of T 42 = 42 • C in the entire tumor volume while trying to minimize the temperature increase in surrounding healthy tissue above T 42 . Such a situation is referred to as optimal heating. The analysis of the choice of P ensuring tumor hyperthermia with minimal overheated normal tissue volume was performed for two beam sizes: one corresponding to the transversal tumor size and one exceeding it for two times aiming to ensure sufficient heating of the tumor boundary region. Figure 6 shows axial and transversal profiles of the temperature distribution for r = 2.5 mm and power 110 mW (corresponding intensity I = 560 mW/cm 2 [50]) in biotissue with and without SiNPs. This combination of r and I parameters ensures heating of the entire tumor volume to temperatures above T 42 in the presence of eSiNPs at concentrations higher than 5 mg/mL, although it is worth mentioning that, besides quite a high value of intensity itself, it requires longer than 400 s to reach the stationary heating mode considered in simulations [5,39]. For the tumor without SiNPs, the employed beam intensity allows reaching T 42 at the beam axis through the whole tumor length (Figure 6a). However, the side regions of the tumor are left underheated, which can be seen from the transversal temperature profiles at depth z = 1 mm (Figure 6b, black solid line). The profiles have the bell shape owing to high tumor thermal conductivity, which smooths the temperature distribution despite the size of the beam exactly fitting the tumor width, and the intensity distribution is uniform.

Beam Size Equals the Tumor Transversal Size
Embedding eSiNPs in nBCC provides additional temperature increase up to 4.6 • C at the laser beam axis for a maximum mass concentration of 7 mg/mL (Figure 6a, solid red line) and broadens the transversal and longitudinal sizes of the region heated above T 42 , thus providing sufficient heating of the entire tumor for all considered concentrations (Figure 6b). The application of wSiNPs results in less efficient additional heating of the tumor: though the whole volume is heated over T 42 , the maximal temperature is reached only within a few hundred microns at the tumor top, and the back half of the tumor is sufficiently less heated compared to the front part. Note that the maximal temperature value obtained with the largest concentration of wSiNPs (7 mg/mL) is almost the same as for the eSiNPs at a concentration of 5 mg/mL. These results can be explained by the increase of absorption efficiency of eSiNPs with their concentration while the increase of C m has almost no effect on the absorption efficiency in the presence of wSiNPs (Figure 4c).
It is worth noting that it is impossible to avoid overheating of surrounding normal tissues. However, it primarily occurs in the superficial layer above the tissue. Despite the decrease in absorption in this region owing to the presence of SiNPs in the tumor (see Figure 4a), a thermal flow from the tumor region results in its overheating ( Figure 6). Nevertheless, the overheating of the normal tissues underlying the tumor does not exceed 2 degrees and quickly decays in-depth, while in the transversal direction the overheating of normal tissue is negligible (Figure 6b).

Beam Size Exceeds the Tumor Transversal Size
The previous considered case for r = 2.5 mm demonstrated the potential of employing SiNPs for the enhancement of laser hyperthermia. However, the required intensity of 560 mW/cm 2 seems to be extremely high from the point of view of the safety requirements which demands intensities of about 200 mW/cm 2 or less [50]. Obviously, the increase in beam radius allows to deliver the same energy to the tissue with lower intensity. Figure 7 depicts the temperature profiles for beam with radius of 5 mm and laser power of 165 mW. The corresponding intensity is 210 mW/cm 2 , which is significantly lower than the maximum permissible radiation intensity in the order of 500 mW/cm 2 for only heating effects to occur.
In this case, the maximal additional temperature increase due to embedding SiNPs is 1.9 • C achieved when adding eSiNPs at concentration of 7 mg/mL. Figure 7a shows that the presence of SiNPs in tumor volume allows to reach temperatures higher than T 42 in the entire tumor volume given that without SiNPs this threshold temperature is not reached. However, it can be expected that due to the complex structure of biotissue the temperatures in real tissue may vary, and the obtained increase is not so significant as compared to the value of 4.5 • C obtained in the case of a narrow beam. A comparison with Figure 6 shows that employing the beam with the diameter exceeding the tumor width allows to arrange more uniform transversal temperature distribution in the tumor region. Figure 8 shows the dependency of temperature in the center of the tumor (x = 0, y = 0, z = 1 mm) and in healthy tissue at the same depth neighboring the tumor (x = 0, y = 3 mm, z = 1 mm) for both considered beam sizes on laser beam power. A smaller beam size allows to reach a more pronounced temperature contrast between tumor tissue and healthy tissue of up to 5 • C compared to 1.5 • C for the larger beam at comparable beam power. We stress that in the case of anti-tumor treatment the ANSI requirements for maximum permissible intensity (200 mW/cm 2 ) [50] can be violated in the tumor region provided that the surrounding healthy tissues are not over-exposed.  Therefore, higher mismatch between temperatures in tumor and normal tissue in the case of smaller beams provides wider opportunities for choosing a beam power for fulfilling optimal heating conditions. Note that the smaller temperature mismatch in the wide beam case is conditioned by the choice of the location in healthy skin within the area of direct laser irradiation by the 5-mm radius beam, while in the case of a narrow beam the temperature increase in healthy skin is caused by thermal conductivity.

Comparison of Flat Beam and Gaussian Beam Irradiation
The Gaussian beam is a widely used alternative to a flat intensity distribution in tissue irradiation configuration. The inhomogeneity of intensity distribution is considered as a drawback in laser hyperthermia due to possible hot spots in the target area [92]. A comparison of these two irradiation regimes evidences the longer time necessary for cancer cell necrosis in the latter case [99]. Nevertheless, Gaussian laser beams are often employed for cancer theranostics (e.g., [100]). To compare effects for both cases, we carried out simulations for the Gaussian laser beam with the same parameters as were employed for the flat beam.
In these simulations, the laser beam axis coincides with the axis of the tumor. Incident laser radiation power depends on x and y coordinates as: where r is the beam radius at 1/e level and P 0 is the total laser power. In Monte Carlo simulations, to obtain Gaussian beam intensity distribution, the Box-Muller transform [101] was employed. The obtained cross-sections q(x, y = 0, z) of the normalized volumetric density of absorbed energy is presented in Figure 9. As one can see, qualitatively, the absorption distribution for the Gaussian beam does not differ significantly (cf. Figures 5a,d,e and 9a-c) from the one for the flat beam, except for the region outside the tumor. Embedding in tumor eSiNPs with smaller absorption results in heating the entire tumor, whereas in the case of wSiNP the rear tumor part is underheated. Detailed calculations of the temperature distribution inside the tissue ( Figure 10) evidence that at equal power the Gaussian laser beam enables heating up to lower temperature than the flat one: the tumor is not completely heated up to 42 • C at P 0 = 110 mW and r = 2.5 mm (precisely of the tumor diameter) (cf. Figures 6 and 10a,b) and precisely up to 42 • C at P 0 = 165 mW and r = 4 mm (cf. Figures 7 and 10c,d, note that the beam radius for uniform circular beam is larger than that for the Gaussian beam, which indicates even less tumor heating efficiency in the latter case). The reasons for this effect are obvious: (i) more than 37% of the laser radiation is absorbed outside the tumor if the beam radius exceeds the tumor radius, and (ii) due to a strongly non-uniform heat rate the heat flux is more prominent than in the case of a uniform intensity distribution. Thus, the simulation of the hyperthermia by Gaussian beam evidences that the circular flat laser beam has certain advantages over the Gaussian one. Since no experimental results are available to assess the validity of the calculations performed, we analysed similar simulation studies. Refs. [37,38,43] report on the numerical simulation of heating water suspensions of SiNPs. However, such studies do not account for the effect of surrounding biotissues both in light transport and heat transfer aspects.
Recently, various inorganic metallic compounds were tested for their effectiveness as thermal coupling agents. In [102], CuFeSe 2 nanosheets with a mean hydrodynamic diameter of 200 nm were experimentally shown to be able to heat water to 30 • C within 600 s. In [103], palladium hydride nanoparticles suspended in water gained additional heating at 10 • C during 180 s at NIR laser power of 154 mW. However, they require a sophisticated fabrication technique, while the pulsed laser ablation of silicon nanostructures in liquids seems to be a promising, scalable technology. As for inert material studies, gold nanoparticles irradiated with NIR light permit selective tumor hyperthermia at 42−43 • C without side effects to surrounding healthy tissues due to the spatial and temporal control of heat generation (selective accumulation and irradiation of heat nanosources) [104]. Nevertheless, silicon demonstrates more prominent biodegradation compared to gold, since silicon presents in living organisms as a natural trace compound.

Conclusions
This paper reports on the results of numerical simulations of human skin tumor hyperthermia mediated by the targeted delivery of SiNPs. Two types of SiNPs formed by the laser ablation of silicon nanowire arrays in water (wSiNPs) and ethanol (eSiNPs) are considered as potential enhancement agents. The simulations were performed for two radii of irradiation beam (one fitting the tumor size and second exceeding it) in order to reveal the optimal conditions of heating the tumor above the threshold temperature while minimizing excessive heating of surrounding normal tissues. A parameter for the evaluation of volume absorption efficiency depending on the optical properties of biotissue and SiNPs was proposed as a measure of NPs efficiency for hyperthermia performance, and the influence of SiNPs concentration on the light absorbance by the tumor was analyzed for both wSiNPs and eSiNPs. The effect of shielding deeper tumor regions from radiation by embedded SiNPs with high absorption coefficient was demonstrated, indicating the requirement for a careful selection of SiNPs concentration in tumor. It was demonstrated that irradiation by the beam with a smaller radius allows to obtain a higher temperature contrast between the tumor and surrounding normal tissues (about 5 • C) than in case of a larger beam radius (1.5 • C) for the comparable beam power at which healthy tissue neighboring the tumor remains underheated, while eSiNPs provide a higher additional heating effect as compared to wSiNPs for the same concentration of nanoparticles in the tumor.
As a result, the laser-ablated SiNPs seem to be promising for the enhancement of photohyperthermia. In comparison to gold nanoparticles, widely discussed as hyperthermia agents, the SiNPs are less expensive in terms of raw material and are shown to be biodegradable. However, laser ablation technology remains relatively expensive for use in mass production. This problem may be solved by using silicon targets with a reduced ablation threshold, e.g., the SiNWs discussed in this article, as well as by the optimization of laser source parameters. In the latter case, relatively cheap, high-power nanosecond lasers have high potential to produce a high yield of SiNPs with desirable sizes and properties, which can be considered as an alternative for the more expensive femtoand picosecond lasers that are currently widely used. Additionally, it will be necessary to carry out experiments with real biological tissues by the photohyperthermia technique to confirm our calculations.

Data Availability Statement:
The data used in this research is available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.