Unified Model for Laser Doping of Silicon from Precursors

Laser doping of silicon with the help of precursors is well established in photovoltaics. Upon illumination with the constant or pulsed laser beam, the silicon melts and doping atoms from the doping precursor diffuse into the melted silicon. With the proper laser parameters, after resolidification, the silicon is doped without any lattice defects. Depending on laser energy and on the kind of precursor, the precursor either melts or evaporates during the laser process. For high enough laser energies, even parts of the silicon’s surface evaporate. Here, we present a unified model and simulation program, which considers all these cases. We exemplify our model with experiments and simulations of laser doping from a boron oxide precursor layer. In contrast to previous models, we are able to predict not only the width and depth of the patterns on the deformed silicon surface but also the doping profiles over a wide range of laser energies. In addition, we also show that the diffusion of the boron atoms in the molten Si is boosted by a thermally induced convection in the silicon melt: the Gaussian intensity distribution of the laser beam increases the temperature-gradient-induced surface tension gradient, causing the molten Si to circulate by Marangoni convection. Laser pulse energy densities above H > 2.8 J/cm2 lead not only to evaporation of the precursor, but also to a partial evaporation of the molten silicon. Without considering the evaporation of Si, it is not possible to correctly predict the doping profiles for high laser energies. About 50% of the evaporated materials recondense and resolidify on the wafer surface. The recondensed material from each laser pulse forms a dopant source for the subsequent laser pulses.


Introduction
Edge isolation [1], texturing [2], ablation [3,4], cutting [5], and annealing of implanted silicon [6,7] are some of the common applications of laser processing for solar cell fabrication. Laser doping of silicon with a precursor layer that serves as a dopant source [4,[8][9][10] shows significant advantages compared with conventional furnace diffusion processes. Laser processing enables the selective doping of silicon with high spatial resolution.
Laser doping of silicon relies on melting the silicon's surface with a short laser pulse, which melts the silicon. Then, doping atoms (such as boron or phosphorus) contained in a thin precursor layer diffuse into the silicon melt. After the (defect-free) recrystallization of the Si melt, the doping atoms are contained in the solid Si. The detailed tailoring of the doping profiles requires not only a qualitative but also a quantitative understanding of the physics and chemistry of laser doping. Only under these conditions is it possible to tailor and improve the electrical and electronic properties of the doped surface layers.
In principle, there are three possibilities for the doping from the precursor layer into the silicon.
(i) Solid and Melt/Melt-Doping: Doping of the Si melt from the solid precursor; in this case, the doping precursor does not evaporate but is solved in the Si melt. The Si itself does not evaporate. For example, the doping of Si with the help of a pure, sputtered boron containing the precursor falls in this range. Boron melts at temperature T melt, B 2 O 3 = 2573 K [11] and evaporates at T evap, B 2 O 3 = 4203 K [12]; the evaporation of the precursor layer is not relevant for the doping process. Lill et al. [13] reported on this case. They used a boron precursor, a pulsed laser of 532 nm wavelength and a 42 ns pulse duration (at full width at half maximum (FWHM) and medium laser pulse energy densities. For these experimental conditions, their model [13] explained the measured doping profiles well. Unfortunately, the model is not valid if either the precursor or the Si evaporate.
(ii) Gas/Melt-Doping: Doping of the Si melt from the gas phase; in this case, the precursor evaporates and doping of the Si melt takes place from the gas phase. This process resembles classic doping using a doping furnace and a gaseous doping source. In this case, the Si itself also does not evaporate. This kind of doping takes place if one uses either a pure, sputtered phosphorus layer [14] or a spin-on phosphoric acid solution [15]. Phosphorus evaporates at 550 K [12] and crystalline silicon melts at temperature T melt, Si = 1687 K [16]. Köhler et al. [14] also used a pulsed green laser of 532 nm wavelength (7 ns FWHM pulse duration). In their publication [14], they used a model that resembled the description of classic furnace diffusion from a phosphorus silicate layer; their model predicted the experimental diffusion profiles well. However, it is only applicable to laser energy densities below the evaporation threshold of Si. In a similar way, the model of Blecher et al. [15] for experiments using a continuous-wave (CW) laser with 532 nm wavelength is also not applicable to pulsed lasers with a pulse duration in the nanosecond regime. In this case, however, the doping profiles were well predicted for different scanning speeds.
(iii) Gas/Gas-Melt-Doping: Doping of the Si melt from a mixture of evaporated doping source AND evaporated Si. Thus, in this case, parts of the silicon's surface evaporate too. This kind of doping takes place in the case of medium to high laser energy densities. In the following, we show that this range starts already at about 2.8 J/cm 2 .
The present work presents a unified model, which is able to cover all three possibilities (solid/melt, gas/melt, gas/gas-melt) of laser doping, including the high energy regime in which the doping precursor as well as parts of the silicon's surface evaporate. We developed a Matlab code for solving the coupled heat and diffusion equations. The present model enables us to predict the melt width W melt and depth d melt of the silicon's surface. Furthermore, it allow us to predict the diffusion profiles over a wide range of laser pulse energy densities H. Figure 1 shows three regimes for the measured sheet conductance G sh of a laser-doped silicon surface. In this case, we dope the silicon from a sputtered boron oxide precursor. Details of the laser system will be discussed below.

Pulse Energy Density Regimes
To measure the sheet conductance G sh of the doped silicon surface, an n-type 6inch Czochralski-grown Si Wafer having a base resistivity ρ base = 1.5 Ωcm and thickness d wafer ≈ 158 µm is provided with a sputtered boron oxide (B 2 O 3 ) layer of 1 nm thickness, covered with 12 nm sputtered amorphous silicon (a-Si) layer. The surface is then irradiated with a pulsed laser and the sheet resistance R sh is measured. The sheet conductance depends on the sheet resistance through For H < 1.5 J/cm 2 , no sheet conductance G sh is measurable. The laser energies are below the melting thresholds of Si; therefore, doping is not possible. In regime I (1.5 J/cm 2 ≤ H ≤ 2.8 J/cm 2 ), the sheet conductance G sh of the doped silicon increases linearly, which indicates that the amount of incorporated boron atoms in the doped layer also increases linearly. For 2.8 J/cm 2 < H ≤ 4.0 J/cm 2 , in regime II, the sheet conductance G sh saturates. This regime is the most important one for technological applications. We will show below that in this regime, parts of the doped silicon surface evaporate. In regime III (H > 4 J/cm 2 ), the sheet conductance G sh strongly decreases again. The reason for the behavior of G sh for H > 4 J/cm 2 lies in the enhanced evaporation of the doped silicon.

Laser Beam
Figure 2a describes the shape and intensity I distribution of the laser beam used in this work. A Q-switched, frequency-doubled Nd:YAG laser (532 nm wavelength) emits a pulsed laser beam with pulse durations between 41 ns ≤ τ p ≤ 110 ns at FWHM. The beam is shaped to a focused line. The width W p of the beam's short axis (measured in the x-direction at 1/e 2 = 0.135 of the maximum intensity I) is W p = 12 µm and the length l p of the long axis (in the z-direction) is l p = 280 µm. The intensity I, along the short (x-)axis, has a Gaussian distribution and a tophat distribution along the long (z-)axis.
We derive the precise value for the laser spot from the following procedure: A chemically and mechanically polished float zone Si-Wafer surface coated with an 80 nm-thick silicon nitride (SiN x ) is used for measuring the laser spot width. The pulsed green laser irradiates the surface with different pulse energy densities with no overlap between locally melted areas. Silicon melts and evaporates during irradiation. As a consequence, the evaporating silicon surface ablates the SiN x layer. Using a laser scanning microscope, we measure the ablated diameter D abt and plot D abt 2 against the logarithm of the pulse energy density. From the slope of the linear fit to these data points, we calculate the pulse width at 1/e 2 of the maximum intensity I. Figure 2b schematically shows the used specimen for investigating the surface topography of silicon wafer after the laser process. For the purpose of getting a flat surface, a chemically and mechanically polished n-type 4-inch float zone grown silicon wafer (FZ-wafer) (Si-Mat Silicon Materials, Kaufering, Germany) is used for these experiments. A 5% hydrofluoric acid dip removes the native SiO 2 layer on the silicon surface. The wafer is placed on an xz-translation stage. The laser beam irradiates the silicon surface with dif-ferent pulse energy densities H. The pulse energy density H is calculated from the average laser power P av , the laser focus area A spot , and the pulse repetition rate f according to The laser beam has a Gaussian intensity I distribution in x-direction and a tophat intensity I distribution in z-direction. (b) The wafer to be irradiated is placed on an xz-translation stage. The translation stage moves in x-direction with a scanning speed v scan = 300 mm/s; the laser irradiates the sample with a pulse repetition rate f = 12.5 kHz. Each laser pulse locally melts the silicon surface. Such a high scanning speed v scan ensures no overlap of the irradiated and melted areas; therefore, the individually treated spots are clearly separated. (c) The melted surface resolidifies after the pulse ends, leaving behind a deformed surface, as discussed later.
We use A spot = 12 × 280 µm 2 and f = 12.5 kHz. Each laser pulse melts the silicon surface locally. For the investigation of the surface patterns, we scan the laser with a scanning speed v scan = 300 mm/s by moving the xz-stage in the x-direction. This scanning speed is high enough to ensure no overlap between the locally melted areas on the silicon's surface. The overlap ratio O x between two subsequent irradiation pulses depends on the pulse repetition rate f , the scanning speed v scan , and the pulse width W p as From Equation (3), using a scanning speed v scan = 300 mm/s and a frequency f = 12.5 kHz, we obtain an overlap O x = −100%. Therefore, in x-direction, the distance between two neighboring pulses is just equal to the width of the pulses itself. Figure 2c shows the profile of the silicon surface after resolidification. The irradiated surface is investigated using a laser scanning microscope (LSM). Figure 3 illustrates the different phases of the laser doping process. Figure 3a shows the sample prior to irradiation. A Leybold Z550 sputtering machine (Leybold GBmbH, Cologne, Germany) deposited a two-layered precursor on a polished, chemically-cleaned 4-inch n-type Czochralski Si wafer (Cz-Si wafer) with base resistivity ρ base = 1.5 Ω.cm and thickness d wafer ≈ 160 µm. The machine uses a radio frequency sputtering technique (RF equals 13.56 MHz) and argon as sputtering gas. The sputtered precursor layer stack consists of a highly hygroscopic B 2 O 3 layer [17] and an a-Si layer. The a-Si layer protects the boron oxide layer from humidity during laser processing in ambient air. The thickness of the boron oxide layer measured by profilometer is 1 nm and the a-Si layer measured by ellipsometer is 12 nm thick. The boron oxide layer serves as the dopant source for boron atoms during laser processing.   Figure 3c shows the situation after one single laser shot: the doped molten silicon resolidifies and the evaporated materials are partially recondensed. In this case, they are recycled in the next shot. Figure 3d demonstrates the situation after the second laser pulse. The subsequent laser pulse irradiates and melts the silicon surface again: the not-yet-irradiated precursor layer and the recondensed material. The recondensed, recycled materials have to be considered in the model when one aims at predicting the incorporation and diffusion of boron atoms into the liquid silicon. Figure 3e finally delineates the doped silicon surface and the recondensed material covering it. In our laser-doping experiment, we irradiate four fields with four different pulse energy densities H 1 = 1.77 J/cm 2 , H 2 = 2.31 J/cm 2 , H 3 = 3.91 J/cm 2 , and H 4 = 4.25 J/cm 2 . After irradiation, using a diamond pen, the wafer with the four fields is cut into pieces of a size A sample = 8 × 8 mm 2 for the secondary ion mass spectroscopy (SIMS) measurements.

Experimental Results
This section validates the present model via a comparison with the experimental data. Figure 4 shows the deformed surface profiles captured by the laser scanning microscope of the FZ-wafer after single-pulse irradiation using laser pulse energy densities in the range 1.77 J/cm 2 ≤ H ≤ 4.97 J/cm 2 . For simplicity, not all of the captured surface profiles are included in Figure 4. However, the omitted data are finally included in the comparison of the measured and modeled data for the deformed width W def and the evaporation depth d evap .  i) Figure 4. Laser scanning microscope profiles of the float zone grown silicon wafer's surface after irradiation using single, well-separated laser shots with different laser pulse energy density H. The silicon surface deforms during laser irradiation, leaving behind a deformation profile of width W def . Figure 5 compares the measured deformation widths W def with the calculated melt widths W melt from our simulations. (a) For H = 1.77 J/cm 2 , the surface deformation is barely detectable. (b) At H = 2.38 J/cm 2 , the surface has two tiny satellite peaks with a relatively higher peak in the middle of the deformation surface profile. The tiny peaks stem from the capillary wave excited by the thermocapillary convection [19]. The peak in the middle is the result of the density anomaly of silicon during melting and resolidification [19]. (c,d) Increasing H boosts the thermocapillary convection within the silicon melt; as a result, the amplitude of the deformed surface peaks increases. For H > 2.8 J/cm 2 , as stated in Figure 1, parts of the silicon's surface evaporate. (e) For H > 4 J/cm 2 , the depth of the molten silicon in the evaporated zone significantly increases, causing a deep tubelike capillary. Through multiple reflections of the laser in the deep capillary, the absorption of laser radiation correspondingly increases [20]. (f-i) For even higher H, more and more silicon evaporates. As a consequence, the height of the two satellite peaks and the depth d evap of the missing, evaporated volume in the center increases. The deep evaporated surface hinders the creation of a middle peak via the solidification by the density anomaly like in (e). Figure 4a shows the the surface deformation for H = 1.77 J/cm 2 . The surface deformation is barely visible. Figure 4b represents the surface profile when using H = 2.38 J/cm 2 . The deformed surface profile has two tiny satellite peaks with a relatively higher peak in the middle of the deformation surface profile. The capillary wave excited by thermocapillary convection [19] leads to the two satellite peaks. The middle peak stems from the density anomaly of silicon during melting and resolidification [19]. Figure 4c,d demonstrate that the height of the three peaks of the surface profile increase with increasing the used laser pulse energy densities H = 3.46 J/cm 2 and H = 3.97 J/cm 2 . Here, it is worth mentioning that parts of the silicon surface evaporate when H exceeds 2.8 J/cm 2 , as already indicated in Figure 1. This evaporation is also confirmed by our simulations. Figure 4e represents the deformed surface profile when increasing the used laser pulse energy density H = 4.22 J/cm 2 . On the one hand, the height of the two satellite peaks significantly increases. On the other hand, the height of the middle peak is significantly decreased and sunken below the original surface level. Such a high laser pulse energy density evaporates the molten silicon from even deeper, causing a deep, tubelike capillary. During silicon resolidification, the density decrease is not sufficient to compensate the evaporated depth, leaving behind a small peak in the middle. Figure 4f shows that increasing the irradiating pulse energy density increases the height of the two satellite peaks even more. However, the middle peak disappears, as the evaporated tubelike capillary becomes deeper. The deeper-evaporated capillary in the middle of the melt hindered the density anomaly of silicon from creating a middle peak like in Figure 4e. Figure 4g-i illustrate that not only the amplitude of the two satellite peaks but also the depth d evap of the evaporated capillary in the middle further increases when increasing the used laser pulse energy density to be H = 4.22 J/cm 2 . The width of the deformation W def (the outer distance between the satellite peaks at the level of the wafer surface) increases with the pulse energy density H. For H > 4.22 J/cm 2 , the depth of the evaporated tubelike capillary d evap strongly increases with laser pulse energy density H. Figure 5a shows the measured melt depth d melt extracted from the measured concentration profiles in Figure 7 in comparison to the calculated d melt in the two cases, ignoring and taking the evaporation of silicon during irradiation into consideration. In regime I (H < 2.8 J/cm 2 ), silicon does not evaporate, and thus, both measured and calculated melt depth d melt agree well. However, for H > 2.8 J/cm 2 (regime II), which evaporates parts of silicon's surface, a significant overestimation occurs due to the fact that the excess absorbed energy deepens the melted volume without being consumed in the evaporation process. The estimated melt depth d melt seems to continue increasing linearly with pulse energy density H. A much closer aggrement with the measured melt depth d melt results from considering the effects of the evaporation of silicon during irradiation. A saturation-like relationship of the melted silicon depth d melt appears to happen for 2.8 J/cm 2 < H < 4 J/cm 2 . For H > 4 J/cm 2 (regime III), as represented in Figure 4e, the evaporation in the middle significantly increases and thus, the total melt depth d melt increases. Figure 5b also compares the measured surface deformation W def using LSM profiles in Figure 4 with the calculated melt width W melt in cases of considering and ignoring the evaporation of silicon during irradiation. The melt width W melt should be equal to the resolidified deformed width W def . In regime I (H > 2.8 J/cm 2 ), there is a reasonable agreement of the calculated width of the melted surface W melt with the measured deformed width W def . On the one hand, an overestimation of W melt results from ignoring that parts of the silicon surface evaporate during irradiation with H > 2.8 J/cm 2 (regime II). The overestimation is due to that the excess absorbed energy widens the melted volume without being consumed in the evaporation process. On the other hand, considering the evaporation effects in the model, yields a closer agreement of the melt width W melt with the measured deformed width W def . For H > 4.25 J/cm 2 (regime III), the calculated melt width W melt is not able to correctly predict the deformed width W def . A possible reason for that is that the real or effective thermal conductivity k Si of silicon is larger than the literature value. As a result of the enhanced absorption of radiation in the tubelike capillary due to multiple reflection, the melt mixes vigorously, which results also in an enhancement of the heat transfer process. Consequently, the melt width W melt increases as a result of the increased effective thermal conductivity k eff probably due to turbulent melt motion as a result of the absorbed energy. One more possible reason is the direction of laser radiation absorption in the capillary. The laser radiation is absorbed laterally on the walls of the vertical capillary. As a result, the expansion in the lateral direction increases more than that in the case of vertical absorption of radiation. Finding the effective thermal conductivity k eff for the mixing-enhanced heat transfer is beyond the scope of this work. Figure 5a,b validate the capability of the present model to predict the dimensions of the melt pool in the three regimes-regime I, II, and III (pulse energy density range 0 ≤ H ≤ 4.25 J/cm 2 ). Using H > 4.25 J/cm 2 in the laser doping induces defects in the doped silicon's surface, which in turn, impairs the electronic properties of the doped area. Therefore, for the laser doping process, the pulse energy density range H > 4.25 J/cm 2 is irrelevant.  (b) In regime I (H ≤ 2.8 J/cm 2 ), the calculated melt width W melt agrees well with the measured surface deformation width W def from the measured profiles in Figure 4. When parts of the silicon's surface evaporate for H > 2.8 J/cm 2 (in regime II), the simulation that considers the evaporation effects fits the measured values much better than the one that ignores evaporation. In regime III, for H > 4 J/cm 2 , the calculated melt widths deviate from the measured deformation width W def . A possible reason is that the model does not consider the enhancement of heat transfer through circulation as a result of the Marangoni effect. Another possible reason is that the absorption of laser radiation in the evaporated capillary is stronger in the x-direction than in the y-direction.

Melting and Resolidification
The correct prediction of the time-dependent temperature is crucial for the simulation of the diffusion of the doping atoms into the liquid silicon, during the heating and cooling periods of the melt. Not only is the overall duration of these two periods important, but also the prediction of the detailed time dependence, because the diffusion coefficient of the doping atoms depends sensitively (thermally activated) on temperature, which will be discussed later. Figure 6 compares the calculated temperature-time curves resulted from ignoring and considering the evaporation effects during laser melting of the silicon's surface. The simulation without considering evaporation of Si was done as follows: We assumed that the only existing latent heat H Si of silicon is that for fusion H melt, Si . Similarly, we assumed also the heat capacity C P Si of liquid silicon to be valid for temperatures even beyond the nominal melting (T melt, Si = 1687 K [16]) and evaporation (T evap, Si = 3538 K [21]) temperature of silicon.    Figure 6a represents the calculated temperature-time profile for H = 1.77 J/cm 2 when considering and ignoring evaporation of silicon. At H = 1.77 J/cm 2 (regime I), no evaporation of silicon takes place, and thus, the two curves agree well. The maximum reached temperature T max in both cases is the same. Not only the heating phase but also the cooling phase in both cases is the same. The figure shows the heating phase to be the duration within which the temperature of the center of the molten surface T surf increases from 2100 K to the maximum temperature T max . The maximum temperature T max is the highest temperature reached during the whole process. The maximum temperature T max increases with the used pulse energy density H to irradiate the surface. The cooling phase is the phase within which the temperature of the center of the melted surface T surf decreases from the maximum temperature T max down to 2100 K. Figure 6b compares the calculated temperature-time curve when considering the evaporation of silicon during irradiation of the silicon surface with H = 2.90 J/cm 2 (regime II). Ignoring the evaporation of silicon overestimates the calculated maximum temperature T max and the cool-down duration. The duration that the melt needs to completely resolidify is called the cool-down duration. Considering the evaporation of silicon keeps the temperature T surf less than or equal to the evaporation temperature T evap, Si . The excess energy absorbed evaporates the silicon, keeping the surface temperature T surf ≤ 3538 K [21]. The evaporation phase is the phase within which the temperature of the center of the molten silicon surface remains constant at the evaporation temperature of silicon T evap, Si . Figure 7a presents the measured and the simulated doping profiles for the specimen described in Figure 3. For doping at H 1 = 1.77 J/cm 2 and H 2 = 2.31 J/cm 2 , both in regime I, our model simulates the doping cases melt/melt and gas/melt (cases i and ii in Section 1) well. Amorphous silicon (a-Si) melts at 1420 K [22], while boron oxide (B 2 O 3 ) melts and evaporates at 723 K [18] and 1773 K [18], respectively. Subsequently, a-Si and B 2 O 3 melt before crystalline silicon c-Si melts (T melt, Si = 1687 K [16]). Thus, melt/melt doping occurs between the two melted species, a-Si and B 2 O 3 . Due to the evaporation of B 2 O 3 , the gas/melt diffusion takes place within less than 1 ns after the c-Si melts and continues until the melt resolidifies (more details are given in the Appendix C). The agreement between the measured and the simulated curves with boron oxide evaporation validates the model's ability to simulate the melt/melt and the gas/melt doping cases. A wrong prediction of the surface concentration as well as the doping depth occurs as a result of ignoring the effects of the evaporation. However, the silicon itself does not evaporate yet in this regime. Boron oxide evaporation and expansion (concentration loss) decrease the concentration gradient at the melted silicon surface. Consequently, ignoring the evaporation process of the precursor incorporates more boron atoms in the calculated doping profile. Figure 7b further proves the ability of our model to simulate the doping case iii (gas/gas-melt). In this regime II, for H > 2.8 J/cm 2 , not only boron, but also silicon evaporate, as stated by Figure 1. Doping with H 3 = 3.91 J/cm 2 (regime II) and H 4 = 4.25 J/cm 2 (regime III) includes all of the three abovementioned doping cases (case i, ii, and iii in Section 1). For such high pulse energy density, a part of the doped, melted silicon surface evaporates (T evap, Si = 3538 K [21]). The previously evaporated B 2 O 3 , followed by the recently evaporated doped Si, form a gaseous doping source for the melted nonevaporated silicon (corresponding to doping case gas/gas-melt). The agreement between the simulated profiles, which consider evaporation not only for boron but also for Si, and the measured doping profiles give direct proof for the validity of our model. Evaporation of the doping source and/or parts of the doped silicon surface decreases the amount of incorporated dopant inside the resolidified material. Consequently, the doping depth is decreased. Simulations that neglect the evaporation of materials overestimate the doping profile's height and depth. Figure 7a,b prove that without considering the evaporation effects of not only B 2 O 3 but also Si, it is not possible to correctly predict the doping profiles for the technologically interesting regime II. The next section describes the numerical simulation in detail.

Numerical Simulation
We exemplify our simulation by modeling the laser doping of silicon from a sputtered boron oxide precursor layer as the source for the boron doping. The boron oxide is covered by a thin a-Si layer. A purpose-developed Matlab code simulates the laser process through solving the 2-D heat equation (4) and the 2-D diffusion equation At any given time t, the model calculates the temperature T of the elements using the input heat (absorbed radiation)q, thermal conductivity k(T), heat capacity c p , and density ρ(T) of silicon. Simultaneously, the model calculates the time-dependent boron concentration C B of each element using the temperature-dependent diffusion coefficient D c (T). The diffusion coefficient D c (T) is described by the Arrhenius equation where D 0 is the pre-exponential factor, E a is the activation energy, and k is Boltzmann's constant.
The following subsections discuss the most important aspects of our model. Figure 8 schematically shows the spatial discretization of the considered cross-section of the silicon wafer and the precursor layer for the numerical simulation. The linefocused laser beam shown in Figure 2a irradiates the precursor layer and the silicon surface. The polished silicon surface reflects almost 37% (at wavelength λ = 532 nm and temperature T = 300 K) of the incident radiation [23]. Silicon has an absorption length l abs ≈ 1.27 µm [23] for wavelength λ = 532 nm at temperature T = 300 K. Using Lambert-Beer's law, for a wafer thickness d wafer = 160 µm, the laser energy transmitted through the wafer is negligibly small. Consequently, heating up and melting of the silicon surface are results of the absorption of the nonreflected part. The tophat distribution of the laser intensity I in the long axis minimizes the variation of the temperature of the heated silicon along the z-axis. Consequently, we model the laser doping process in two dimensions, x (scanning direction) and y (direction of wafer-surface normal). All the material properties used in this work are given in the Appendix A.

Numerical Discretization
An adaptive grid divides the considered volume into different sizes of elements. In the center of each element, a mesh point exists. The center of the observed cross section has finer-sized elements. The mesh size increases in the lateral direction and in the depth of the silicon wafer. The laser beam irradiates the surface in the center of the fine-size area for the numerical simulation. The advantage of adapting the mesh size is to achieve an accurate calculation of temperature and concentration in the region where the laser melts the silicon without significantly increasing the computational power. A detailed description of the implementation of the adaptive grid concept is presented in references [24][25][26].
The smallest element in the simulation area exists exactly underneath the center of the laser focus and has a size ∆x∆y = 25 × 125 µm 2 . In order to get a precise but fast simulation of the doping process, we optimized the mesh size. The biggest mesh size with which we got the same result as for the finer mesh sizes was 125 × 25 nm 2 . A smaller mesh size does not increase the precision, but only significantly increases the computing time. A bigger mesh size impairs the accuracy of the result accordingly. With the chosen mesh size, a simulation of a complete doping profile takes less than 7 h.  Figure 8. Scheme of the cross-sectional spatial discretization for numerical simulation of the silicon/precursor stack. A line-focused pulsed laser with a wavelength λ = 532 nm irradiates the precursor layer and scans the silicon surface in the x-direction. The discretization uses the adaptive grid concept to achieve increased calculation precision in the locations with high temperature gradients with a reduced calculation time and computational power [24][25][26]. The smallest element lies in the center of the laser focus and has a size of ∆x ∆y = 125 × 25 nm 2 .

Laser Absorption
The precursor layer system does not contribute to the absorption of the laser radiation. The absorption length l α of amorphous silicon for radiation with wavelength λ = 532 nm at room temperature (300 K) is around l α = 460 nm (absorption coefficient α a-Si = 21.88 × 10 3 cm −1 [27]). Using the Lambert-Beer law, the absorbed portion of the incident laser radiation in the 12 nm-thick a-Si layer is less than 2.5% of the incident radiation energy. Therefore, the model considers the precursor layer system as transparent for the incident laser beam.
Due to the laser pulse, the surface of the Si wafer heats up and also heats up the boron oxide. The boron oxide evaporates and boron atoms diffuse into the molten silicon surface.
For H > 4 J/cm 2 , the evaporated tubelike capillary at the center of the melted silicon enhances the absorption of the laser radiation. The absorption enhancement stems from multiple reflections of radiation within the evaporated tubelike capillary. In Figure 4e-i, during the irradiation with H > 4 J/cm 2 , the depth of the evaporated capillary in the center of the silicon melt is deep enough to cause multiple reflections of the incident laser radiation [20], and thus, enhance the absorption of laser in the evaporated volume. The simulation considers the enhanced absorption of laser inside the evaporated capillary by a locally reduced reflection R Si, l of the melted surface in the evaporated capillary. The depth of the evaporated capillary d evap in the deformed surface profile for H > 4.0 J/cm 2 , shown in Figure 4, gives a clue about the depth of the evaporated capillary. In the modeling, we use the measured d evap to calculate the equivalent surface reflectance R Si,l of the melted silicon surface in the evaporated capillary. The comparison between the calculated evaporation depth d evap for different values of R Si,l with the measured d evap enables us to adjust the R Si, l value for each pulse energy density H. For H = 4.22 J/cm 2 , no measurable d evap appears in Figure 4e. Therefore, we used the fit curve appearing in Figure A1b to find the corresponding R Si, l -value. The reflectance R Si, l for all H values is lower than the reported literature value (70% at T ≥ T melt Si [23]). Appendix B lists the calculated reflectance R Si,l versus pulse energy density H.
According to Figure 9, the best fit for the concentration profile at H = 4.25 J/cm 2 is obtained when a lower melted surface reflectance R Si,l is used in the evaporated capillary during the evaporation phase. For H = 4.25 J/cm 2 , using only the literature value of the reflectance of the molten silicon surface (70% [23]) underestimates the doping profile depth and overestimates the surface concentration. Using a reflectance R Si,l = 45% during the evaporation phase better predicts the correct concentration profile.

Heat Transfer
On the surface of the solid or liquid Si, convective heat transfer to the air is neglected for the following reasons: At the melting temperature T melt Si = 1687 K [16] of silicon, the thermal conductivity k Si, s = 22.8 W/(m K) of solid silicon and k Si,l = 62 W/(m K) [28] of liquid silicon are significantly higher than the thermal conductivity k air = 0.01 − 0.08 W/(m K) [29] of air.
To account for overheating (during melting with laser of pulse duration τ p = 42 ns) and undercooling (during resolidification), the model uses the melting and resolidification temperatures with 20 K differences from the literature values [30,31], which is a reasonable value for the solid/liquid interface velocity 8 m/s < v int < 10 m/s. Thus, the melting temperature in the model equals T melt Si + 20 K and the resolidification temperature equals T melt Si = T melt Si − 20 K.
On the time scale, upon laser irradiation, the first diffusion takes place between the melted a-Si layer (T melt a-Si = 1420 K [22]) and the melted boron oxide layer (T melt B 2 O 3 = 723 K [18]). Later, after the c-Si melts (T melt Si = 1687 K [16]), B-atoms diffuse only from evaporated B 2 O 3 into the molten c-Si surface because the energy needed to heat up and evaporate a 1 nm-thick B 2 O 3 element equals the energy needed to melt 1.5 nm of c-Si element. The detailed calculation is given in the Appendix. The finest silicon element considered in the numerical discretization in Figure 8 has a depth ∆y = 25 nm. Consequently, the B 2 O 3 element evaporates as soon as the underlying c-Si element melts. Figure 10 shows the different phases that occur during laser irradiation on the surface of the crystalline silicon (c-Si) wafer with the precursor layer (12 nm a-Si + 1 nm B 2 O 3 ). The irradiated system is represented in the form of the discretized elements shown in Figure 8. In Figure 10a, at the center of the laser beam, where the elements with the finest size in Figure 8 exist, the absorbed laser power increases the surface temperature T surf from room temperature until below the melting temperature T melt B 2 O 3 = 723 K [18] of boron oxide. In Figure 10b, the boron oxide melts at T surf ≥ T melt B 2 O 3 = 723 K, while the a-Si layer and c-Si surface are still solid. Figure 10c shows the start of the melt/melt diffusion between the melted boron oxide and the melted a-Si when T surf reaches the melting temperature T melt a-Si = 1420 K [22] of a-Si. More boron oxide melts as T surf increases. When T surf reaches the melting temperature T melt Si = 1687 K [16] of c-Si, the c-Si acquires the latent heat of fusion H fusion Si and melts. In the same time, B 2 O 3 takes up the latent heat of vaporization H vap B 2 O 3 and evaporates, as shown in Figure 10d. As soon as B 2 O 3 evaporates, in the simulation, a reordering of the matrix elements takes place: the evaporated elements go up and the melted elements move down. Figure 10e shows the case for H being high enough to evaporate the silicon (T evap Si = 3538 K [21]). In this case, the surface temperature reaches T surf = T evap Si = 3538 K [21]. At T surf = T evap Si = 3538 K [21], the a-Si layer and a portion of the superficial c-Si elements evaporate and gas/gas-melt diffusion takes place.
723 K < T surf < 1420 K Figure 10. The different phases of the crystalline silicon (c-Si) wafer cross-section and the precursor layers during laser irradiation. The spatial discretization represented in Figure 8 divides the irradiated system into elements. (a) The surface temperature T surf at the center of the beam increases from 300 K until below the melting temperature T melt B 2 O 3 = 723 K [18] of boron oxide. (b) When T surf ≥ 723 K, the B 2 O 3 melts, while the a-Si layer and c-Si surface are still solid. (c) Melt/melt diffusion between the melted B 2 O 3 and the melted a-Si starts when a-Si melts at T surf ≥ 1420 K [22]. (d) At T surf ≥ T melt Si = 1687 K [16], the c-Si melts and the B 2 O 3 evaporates. A matrix reordering of the melted and the evaporated elements takes place; evaporated elements go up, melted elements go down. (e) For pulse energy density H sufficient to evaporate silicon, the surface temperature T surf reaches T evap Si = 3538 K [21]. At T surf = T evap Si = 3538 K [21], a-Si and c-Si elements evaporate and gas/gas-melt diffusion of boron atoms occurs.

Diffusion
Depending on the maximum temperature in the melt, the pre-exponential factor D 0 (in Equation (6)) of the diffusion coefficient of boron in liquid silicon D c B has two values. The temperature T surf is the maximum temperature existing in the melt. At T surf < 2100 K, the pre-exponential factor D 0 equals 2.7 × 10 −4 cm 2 /s (the literature value in Ref. [32]), and at T surf > 2100 K, D 0 equals 8 × 10 −4 cm 2 /s (the reported value in the work of Lill et al. [13]). Appendix D presents a detailed calculation, showing that the diffusion coefficient D c B of boron in silicon must increase during the laser doping process. Figure 11 also confirms that the best fit of the measured doping profiles is obtained when two values for the pre-exponential factor D 0 are used: D 0, 1 = 8 × 10 −4 cm 2 /s (the reported value in the work of Lill et al. [13]) for temperatures T surf > 2100 K; D 0, 2 = 2.7 × 10 −4 (the literature value in reference [32]) for temperatures T melt, Si < T surf < 2100 K.
We choose 2100 K as the transition value for the following reason: this temperature is reached at the surface for H ≤ 1.2 J/cm 2 , which results in no visible surface deformation after irradiation. No surface deformation means that either the surface did not melt at all, or it melted without circulation. For this case, we take the literature value for D 0 . For higher H > 1.2 J/cm 2 , the silicon surface deforms as a result of the irradiation, indicating that Marangoni convection took place. This convection in turn increases the diffusion with a consequence of a higher D 0 . For that reason, the chosen transition temperature is the calculated maximum temperature T max in the center of the melt when a pulse energy density H = 1.2 J/cm 2 is used.  Figure 11. The best fit to the measured doping profiles occurs when using a higher pre-exponential factor D 0 value depending on the maximum temperature T surf . For T surf ≥ 2100 K, the higher pre-exponential diffusion coefficient D 0 = 8 × 10 −4 cm 2 /s is used. After the surface temperature cools down again to T surf ≤ 2100 K, the literature value D 0 = 2.7 × 10 −4 cm 2 /s is used again [32]. Exemplary, using only the literature value of D 0 underestimates the profile depth. A better fit is only achievable when using the combination of both values of D 0 with surface temperature T surf as the condition.

Recondensation
After each laser pulse, parts of the evaporated material (percursor material and for high laser energies also Si) recondense. Figure 12 demonstrates that the calculated doping profiles are fit best assuming a 50% recondensation of the evaporated materials after each pulse. For 40% recondensation, an underestimation of the calculated doping profile depth occurs. For 60%, an overestimation of the simulated doping profile occurs. For the same recondensation of 50%, Eisele and Köhler [30] also obtained the best fits for the measured doping profiles when using phosphosilicate glass as a precursor for laser doping.  Figure 12. A 40% and a 60% recondensation portion of the evaporated materials after each laser pulse yield an underestimation and overestimation of the concentration profiles, respectively. The best fit results from using a 50% recondensation portion.
For H > 4 J/cm 2 , where a deep tubelike capillary of silicon evaporates, the evaporated silicon and boron oxide recondense as a boron-rich a-Si layer. Considering that 50% of the evaporated materials recondense, the recondensed layer, under certain experimental conditions, reaches a thickness of more than 200 nm. A thickness of more than 25 nm (the thickness of the finest crystalline silicon element in Figure 8) contributes to the absorption of laser radiation and heat transfer calculations. For correct calculation of energy absorption and heat transfer, the relatively thick recondensed layer is divided into a corresponding number of sublayers, each having a thickness of ∆y = 25 nm.

Discussion and Conclusions
The present model numerically simulates the laser doping process of silicon with a low evaporation temperature precursor and also for high pulse energy densities H, where silicon evaporates. Previous models were either limited to precursor layers with a high melting and evaporation temperature [13], or they did not work when silicon evaporates due to high laser pulse energy densities [13,14,30]. Some models can only simulate laser-processing with a cw-laser [15]. A detailed comparison of simulated with experimental results of melt width and depth validates the capability of our model to predict the dimensions of the melt pool in the three regimes, regimes I, II, and III, until H = 4.25 J/cm 2 . The very good agreement between the calculated and the measured doping profiles for the investigated wide pulse energy density 1.77 J/cm 2 ≤ H ≤ 4.25 J/cm 2 range proves that the correct predictions of the doping profiles are only achievable when considering the evaporation effects of the precursor layer and Si. Finally, we showed that the best fits of the doping profiles occur when 50% recondensation of the evaporated materials, reduced reflectance in the evaporated capillary during the evaporation phase, and a higher diffusion coefficient are used.    [18] Appendix B. Aspect: Laser Absorption Figure A1a represents the measured depth d evap of the evaporated capillary tubes shown in Figure 4f-i. The higher the pulse energy density H is, the deeper the evaporated depth d evap becomes. Although silicon evaporates for H > 2.8 J/cm 2 (as stated in Figure 1), the evaporated capillaries firstly show up for H > 4.0 J/cm 2 (starting from the profile in Figure 4e). For 2.8 J/cm 2 < H < 4.39 J/cm 2 range, the middle peak in the center of the deformed profile resulting from the silicon density anomaly [19] does not allow an exact measurement of the evaporated depth d evap . From the fit curve in Figure A1b Figure A1b presents the calculated reduced reflectance R Si, l of the molten silicon surface in the evaporated capillary. Due to multiple reflections, the absorption in the capillary increases, therefore, the overall reflection R Si, l of the Si surface within the capillary decreases. In comparison with the measured depth d evap , simulating with the literature value of R Si, l (70% at T ≥ T melt, Si [23]) yields smaller values for the calculated depths d evap . The measured depth d evap serves as a target value in the determination process of the reduced reflectance R Si, l for each pulse energy density H. The model varies the calculated reduced reflectance R Si, l until the calculated depth d evap matches the directly measurable capillary depths d evap for H > 4.39 J/cm 2 . For H < 4.0 J/cm 2 , the reflectance R Si, l is R Si, l = 70%. From the fit curve plotted, we determine the calculated reduced reflectance R Si, l for H = 4.25 J/cm 2 , since we are unable to measure d evap from the profile in Figure 4e.

Appendix C. Aspect: Heat Transfer
The latent heat of fusion of Si is H fusion, Si = 1802 J/g [16] and the latent heat of evaporation B 2 O 3 is H evap, B 2 O 3 = 4308 J/g [18]. The density and heat capacity of liquid B 2 O 3 at evaporation temperature T evap, B 2 O 3 = 1773 K [18] are ρ B 2 O 3 , l = 1.498 g/cm 3 [36,37] and c p B 2 O 3 , l = 0.4472 J/(g K) [38]. The density of solid silicon at melting temperature T melt, c-Si = 1687 K [16] is ρ Si, s = 2.293 g/cm 3 [34]. The required energy per unit area E B 2 O 3 to heat up and evaporate 1 nm of B 2 O 3 per unit area is The required energy per unit area E Si to heat up and evaporate 25 nm of c-Si per unit area is Therefore, the energy needed to heat up and evaporate a 1 nm-thick B 2 O 3 equals the energy needed to melt 1.5 nm of c-Si. Consequently, the 1 nm-thick B 2 O 3 element evaporates as soon as a 25 nm-thick c-Si element melts. Thus, the diffusion process takes place between boron oxide vapor and the molten crystalline silicon surface.

Appendix D. Aspect: Diffusion
The best fits of the measured concentration profiles of boron after laser doping are obtained when two values for the pre-exponential factor D 0 are used: D 0, 1 = 8 × 10 −4 cm 2 /s (the reported value in the work of Lill et al. [13]) for temperatures T surf > 2100 K; D 0, 2 = 2.7 × 10 −4 (the literature value in reference [32]) for temperatures T melt, Si < T surf < 2100 K.
As a result of the Gaussian distribution of laser intensity I in the x-direction, a temperature gradient results between the hot center and the colder edges. The edges are just melted, or in other words, are at melting temperature. The temperature gradient induces a surface tension gradient, which pulls the hotter (molten) silicon in the center towards the colder edges. The surface tension gradient in the melt induces a circulatory motion in the melt. This phenomenon is called the thermocapillary effect or Marangoni convection [19]. Figure A2 schematically shows the circulatory motion induced in the melt as a result of the Marangoni convection. The element of the highest temperature in the melted silicon is the surface element at the center of the melt. The coldest elements (just melted) lie at the boundary between the liquid and solid silicon. The surface tension σ of the molten silicon depends on the temperature T. The temperature gradient across the molten silicon surface causes a surface tension gradient ∂σ/∂T (≈86 × 10 −6 N/(m K) [39]). The melt circulates from the hot center towards the colder edges as a result of the temperaturegradient-induced surface tension gradient. The circulatory motion enhances the mass transport inside the melted volume, and thus, leads to an apparently increased "effective" pre-exponential factor D 0, 1 . The literature values of D 0 are obtained from a silicon melt of a homogeneous temperature distribution. The reported values of D 0 are 2.4 ± 0.7 × 10 −4 cm 2 /s [40] and 2.7 × 10 −4 cm 2 /s [32].
The equation for calculating the circulation velocity v c [41] is where ∆r melt is the length, over which the temperature gradient extends, and η is the dynamic viscosity of silicon (η ≈ 0.75 mN s/m 2 [42]). Figure A2. Marangoni convection: The colors represent the temperature distribution inside the molten volume. Red color represents the maximum and yellow represents the melting temperature (the coldest liquid silicon). Due to the Gaussian distribution of laser intensity I in the x-direction, the hottest temperatures in the melt occur in the center of the molten surface. The coldest lie at the boundary between molten and solid silicon. The surface tension gradient induced by the temperature gradient mixes the melt. The temperature distribution in the melt is symmetric around the pool axis. As a result, the melt moves by circulating from the hotter center to the colder edges. This phenomenon is also called the thermocapillary effect. Figure A3 shows the approximated shape of half of the melt pool. The circumference is calculated according to the closest geometrical shape. The melt pools are simulated by ellipses. Figure A3a demonstrates that, at H = 1.77 J/cm 2 , the melt pool is a quarter of an ellipse. Its perimeter is L melt, total = (a melt + b melt )(1 + π/4), where a melt = ∆r melt = 3 µm and b melt = d melt = 0.75 µm. Consequently, the total circumference of the half melt pool at H = 1.77 J/cm 2 is L melt, total ≈ 6.7 µm. At H = 4.25 J/cm 2 , the half of the melt pool is approximated by two concentric ellipses. The inner ellipse represents the evaporated portion. The circumference of the pictured half of the melt pool is L melt, total = (a melt + b melt )(1 + π/4) + (a evap + b evap )(π/4 − 1), where a melt = ∆r melt = 5.25 µm, b melt = d melt = 1.73 µm, a evap = ∆r evap = 2.5 µm, and b evap = d evap = 0.48 µm. Therefore, the total circumference of the half melt pool at H = 4.25 J/cm 2 is L melt, total ≈ 11.8 µm. The simulation and validation of the time-dependent temperature during the heating and cooling in Section 4.2 allow us to calculate the maximum temperature reached during the irradiation. Furthermore, we calculate the heating, cooling, and evaporation durations t heating , t cooling , and t evap as well as the melt depths d melt to estimate the circulation velocity v c of the molten silicon. For a pulse energy density H = 1.77 J/cm 2 , the maximum temperature T max reached during the irradiation is T max = 2660 K and, for H = 4.25 J/cm 2 , the maximum temperature T max reaches the silicon evaporation temperature T max = T evap, Si = 3538 K. At H = 1.77 J/cm 2 , the heating phase duration (from 2100 K up to the maximum temperature T max ) is t heating = 20 ns and the cooling phase duration (from the maximum temperature T max down to 2100 K) is t cooling = 36 ns. The melt width 2∆r melt is about 6 µm (∆r melt = 3 µm). The calculated circulation velocity v c during the heating phase for an average melt depth d melt, heating = 340 nm is v c, heating = 12.6 m/s, and that during the cooling phase for an average melt depth d melt, cooling = 600 nm is v c, cooling = 22.3 m/s. The total circulated distance during both heating and cooling phases is L circ = 1.06 µm, which is about 15.8% of the total circumference L melt, total of half of the melted pool (L melt, total = 6.7 µm). For H = 4.25 J/cm 2 , the heating phase duration (from 2100 K up to 3538 K) is t heating = 12 ns, the evaporation phase duration (T max = T evap, Si = 3538 K) is t evap = 46 ns, and the cooling phase duration (from 3538 K down to 2100 K) is t cooling = 65 ns. The half melt widths during the three phases-heating, evaporation, and cooling-are ∆r melt, heating = 2.5 µm, ∆r melt, evap = 5.25 µm, and ∆r melt, cooling = 5.25 µm, respectively. The melt depths during the three phases are d melt, heating = 250 nm, d melt, evap = 663 nm, and d melt, cooling = 413 nm, respectively. Using average temperatures during the heating and cooling phases, the circulation velocity v c during the heating, evaporation, and cooling phases are v c, heating = 10.7 m/s, v c, evap = 26.5 m/s, and v c, cooling = 8.3 m/s, respectively. The result is a total melt circulation distance L circ ≈ 1.9 µm, which is about 16.1% of the total circumference L melt, total of half of the melt pool (L melt, total = 11.8 µm). The circulation action in both pulse energy densities, 1.77 J/cm 2 and 4.25 J/cm 2 , justifies the higher pre-exponential factor D 0 used in the model of Lill et al. [13] and in the present model. The reason is that the same maximum pre-exponential factor D 0 (D 0 = 8 × 10 −4 cm 2 /s) is used for the simulation of the whole H range, as the circulated portion L circ /L melt, total for H = 1.77 J/cm 2 and 4.25 J/cm 2 is almost the same-about 16%.