Infrared backwards laser melting of a silicon wafer

We investigate a method for melting a silicon wafer’s rear side with a pulsed infrared laser (1064 nm) impinging onto the front side. The targeted application for this method is deep laser doping. Our numerical model simulates the evolution of the two-dimensional temperature distribution in the wafer caused by pulsed infrared laser irradiation. The model incorporates the temperature dependent material properties of silicon and the enthalpy-based phase change by means of finite volumes. The simulation yields spacial temperature distributions of the wafer’s cross section at defined time steps. We obtain the laser parameters for a continuous melt depth of 40 μm in a 200 μm thick wafer from the analysis of the


Introduction
The concept to introduce impurities into silicon (Si) by laser melting was investigated for the creation of diodes in the 1960s [1] and also for fabrication of first solar cells in the 1980s [2].This method, commonly known as laser doping, bases on the diffusion of dopant atoms from appropriate dopant precursor layers into the laser induced molten Si.Subsequent epitaxial recrystallization forms the doped areas.The Institute for Photovoltaics (ipv) has successfully applied its patented laser doping process [3] in the last years [4], and recently proved its quality with 22% efficient interdigitated back contacted solar cells [5].Laser doping takes advantage of the diffusion constants of impurities in liquid Si, which are several orders of magnitude higher than in solid silicon [6].Consequently, laserinduced diffusion occurs on a timescale of about 100 ns [7], keeping the thermal budget for the entire wafer minimal.The laser, additionally, offers a high spatial selectivity, because doping occurs only in laser-irradiated areas.Although various combinations of process conditions, dopant precursors and laser parameters have been applied for laser doping, the basic principle is always the same.The laser impinges on the front surface and subsequently melts it.Liquid silicon strongly absorbs the irradiation within a short distance from the surface.Thus, the melting depth is limited by heat diffusion.Consequently, this limits the achievable doping depth.
This contribution presents a fundamentally different approach for laser melting, applicable for deep laser doping.An infrared (IR) laser with wavelength λ = 1064 nm impinges onto the front side of a crystalline Si wafer, but a e-mail: patrick.lill@ipv.uni-stuttgart.demelts the rear side.Hence, we call this method backwards IR laser melting.It promises melt depths above 40 μm.The method was first proposed by Bruel [8], who theoretically showed that the formation of a steep thermal front inside silicon is possible with IR lasers.Although the possibility to melt silicon is briefly addressed in the application section of his contribution, no calculations for this mode were presented.
The present work investigates backwards IR laser melting in detail and gives laser parameters for a continuously molten domain.Our simulation yields the temperature distribution in silicon by solving the transient heat transport equation with the finite volume method.We find that the optimal laser parameter set among our results continuously melts the rear side of a 200 μm thick wafer for m d = 40 μm while the front side is unaltered.The simulation yields a pulse duration p d = 1 μs, beam diameter d 0 = 10 μm and pulse energy Q p = 31 μJ as optimal parameter combination to obtain the necessary irradiance E = 80 MW/cm 2 .

Backwards laser melting
Figure 1 shows a schematic view of the backwards IR laser melting concept, applied for deep laser doping.Silicon is almost transparent for IR-radiation at room temperature, thus IR laser photons impinging on the front surface of a Si wafer pass through it.
At the rear, a thin antimony layer absorbs the transmitted IR radiation energy at the Si/Sb interface.Concerning later experiments, we use antimony as absorption layer for this simulation because it is highly absorbing and also a dopant for silicon.Subsequently, the overlying adjacent silicon heats up through heat conduction.The absorption coefficient α Si of silicon exhibits a strong temperature dependence [9], which triggers significant absorption in silicon itself.For sufficient energy input, silicon melts.Consequently, an interface between the newly formed Si melt and the overlying crystalline silicon develops.With an appropriate process control, this liquid/solid interface propagates forward into the bulk of the wafer, opposite to the direction of the impinging laser beam.Due to this self-sustained moving melt interface, heat conduction no longer limits the melting depth as it does for conventional laser doping.Consequently, backwards IR laser melting enables deep laser doping of the wafer's rear side.

Numerical model
Figure 2 shows the sample geometry of our simulation.A laser beam with spatial Gaussian intensity distribution irradiates the front surface of a crystalline Si sample with thickness L = 200 μm.Our computer program simulates the laser heating process by applying the thermal conduction theory.Due to the axisymmetric laser intensity distribution, it is feasible to calculate the transient temperature distribution in two dimensional cylindrical coordinates.This assumption yields the unsteady state heat equation (with for radius r, depth z and time t with temperature distribution T (r, z, t) and laser heat source q(T, r, z, t).The computational domain comprises the whole wafer thickness L = 200 μm and thin Sb layer (l = 200 nm) along the z-axis and runs from the laser spot center to width W = 80 μm.The temperature dependent material properties heat capacity c p , mass density ρ and thermal conductivity K are derived from literature as discussed below and listed in Tables 1 and 2.

Laser heat source
The laser beam has a spatial Gaussian (N r = 1) and temporal "Top-hat-like" or "Super-Gaussian" (N t = 50) intensity distribution I(r, t) according to: IR laser melting.The laser intensity I(r, z, t) inside the silicon consists of the part of the laser pulse transmitted by the air/Si interface with: and the part I re which is reflected at the silicon/antimony interface, Here, L is the Si wafer thickness, α Si is the temperature dependent absorption coefficient of Si and R the temperature dependent reflectivity at the respective interface.The temperature dependent complex index of refraction n − ik gives the reflectivity R a/b at the interfaces according to: where the extinction coefficient k is related to the absorption coefficient α and the laser wavelength λ.The indices a and b denote the materials (silicon, antimony and air) at the respective interface, and The simulation incorporates multiple reflections inside the silicon wafer as far as the reflected irradiance falls below of one percent of the incoming intensity at the front side.No reflection at the antimony/air interface at the rear side of the wafer is assumed, because no irradiation intensity reaches this interface.The laser heat source term

20104-p3
The European Physical Journal Applied Physics q(T, z, t) equals the absorbed laser power inside the respective material,

Phase change considerations
The enthalpy H exhibits a discontinuity at the melting temperature T m .This corresponds to a delta function of the specific heat c p = ∂H/∂T .The delta function is approximated with a Gaussian distribution to avoid numerical instabilities in our simulation.During the melting of the material, the latent heat H m intake increases the heat capacity within a temperature interval ΔT = 10 K, modeled by: Integration of this curve over the temperature interval yields the latent heat of fusion for the respective material.Due to this approach, the material in our simulation is fully molten only after the temperature exceeds T m + ΔT and fully solidifies only after T m − ΔT is passed.

Boundary conditions
The simulation implements the following boundary conditions at the borders of the computational domain.
3. The Dirichlet boundary condition at the outer border perpendicular to the radius with: is justified because the computational domain is sufficiently large.The surrounding material remains at room temperature T 0 = 300 K during the whole process.

Finite volume implementation
The mathematical model for the transient temperature distribution inside the sample is highly nonlinear.The reasons are the strong temperature dependence of the material properties and the phase changes.Hence, an analytic solution is not possible.We obtain an approximate solution by means of the finite volume method [10].Our model uses central differences and an constant grid discretization (Δr = const and Δz = const).The temperature T i,j is located in the center of each finite cell with area A cell = ΔrΔz.Due to Gauss's theorem, this approach always ensures the conservation of energy for each cell.Consequently, the differential equation ( 1) transforms into a system of equations: for r > 0, and for r = 0.The thermal conductivity K e lies at the Eastern, K w at the Western, K s at the Southern and K n at the Northern interface of the respective cell of Figure 2.However, the thermal conductivity K i,j (T ) of each cell is localized at the cell's center.Thus, the interface thermal conductivities are the harmonic means of the center conductivities.It holds for example that: The time integration is realized with the simple explicit Euler rule: which is only stable if the Courant-Friedrichs-Lewy condition [10]: is met.We developed a computer code which solves equation (12a) and (12b) for an infrared wavelength λ = 1064 nm and an antimony absorption layer.

Material properties
Tables 1 and 2 summarize the material properties of silicon and antimony, which are incorporated in the simulation program.
The expressions for c p , ρ and K of silicon are directly inherited from reference [11].Green [9] suggested a method to extrapolate α, k and n from the given values at 300 K with normalized temperature coefficients.
For antimony, approximations to literature data are used for c p [18], K [17] and ρ [15,16].Due to the lack of literature data for the temperature dependence of antimony's optical properties, we only distinguish between solid [19] and liquid [20] antimony.

Simulation results
The first two subsections of this section analyze the effects of different laser parameters on the temperature evolution on-axis (r-axis).The Section 4.1 discusses the simulated temperature evolutions of the Si/Sb interface (a) and the front surface element (b) of the sample (highlighted elements in Fig. 2).The Section 4.2 presents the temperature distributions throughout the sample along the z-axis at the laser spot center.The implications of the three-dimensional nature of the model (imposed by the Gaussian beam shape), especially with varying beam diameter are addressed in Section 4.3.An energy conservation evaluation in the Section 4.4 ensures confidence in the consistency of our numerical scheme.

Temperature evolution for constant pulse energy
Constant pulse energy pulse energy Q p = 31 μJ and beam diameter d 0 = 10 μm result in different irradiances E for varied pulse durations p d .Figures 3a and 3b present the simulation results of four pulses with varying pulse durations p d = 0.5 μs, 1 μs, 2 μs and 4 μs resulting in different irradiances E = 160 MW/cm 2 , 80 MW/cm 2 , 40 MW/cm 2 and 20 MW/cm 2 .The point when the laser is switched on defines the time axis' origin t = 0.
Figure 3a shows the temperature evolution of the laser spot center grid element at the Si/Sb interface (a).The temperature of the interface elements increases rapidly for all four pulse durations p d as the laser pulse is switched on.After briefly overshooting, the interface element temperature settles in the melting interval T m ± Δ T = 1687 ± 10 K for the three longer pulse durations p d = 1 μs, 2 μs and 4 μs. Figure 3b shows the temperature evolution of the laser spot center grid element at the front surface of the sample (b).Our numerical model ceases to be valid at the boiling temperature of silicon T b = 3538 K [14], thus the simulation for p d = 0.5 μs with corresponding E = 160 MW/cm 2 is stopped even before the laser pulse ends, because the temperature of the grid element (b) at the Si surface quickly exceeded T stop = 3500 K.It is evident, that the irradiance strongly affects the maximum temperatures of the grid element (b) at the Si surface.The temperature distribution along the depth axis, presented in Section 4.2, shows that no heat conduction from the rear side contributes to the surface temperature increase.Due to the low absorption coefficient of silicon α Si = 1.11 × 10 3 m −1 at room temperature T 0 = 300 K, only little absorption occurs inside the sample.Although the laser intensity is highest at the front side, the temperature increase is relatively slow and the maximum temperatures of the surface element (b) remain far below T m = 1687 K for the three longer pulse durations.However, it is also clear from Figure 3b, that the temperature at the surface (b) continuously increases during the laser pulse because of a thermal built-up.Thermal built-up means that more heat is supplied (pulse energy) than diffuses out into the neighboring elements.For the shortest pulse p d = 0.5 μs with highest irradiance E = 160 MW/cm 2 , the front surface heats up sufficiently that the strong temperature dependence of α Si results in an self-amplifying absorption process.As surface alteration by Si melting and evaporation is undesirable for our purposes, parameter sets leading thereto, are excluded from further discussions.
However, the self-amplifying absorption effect is even stronger directly at the Si/Sb interface.The silicon element at the interface (a) heats up (indirectly) via heat conduction, because antimony is highly absorptive for the incoming IR radiation (α Sb = 60.82 × 10 6 m −1 ), which in

20104-p5
The European Physical Journal Applied Physics turn leads to increasing absorption inside the Si element itself.Hence, the self-amplifying absorption commences, which explains the steep temporal temperature gradient at the Si/Sb interface element (a).

Temperature distribution along z-axis
Figures 4a-4c show the temperature distributions at the laser center along the z-axis at different points in time.The beam d 0 = 10 μm and pulse duration p d = 1 μs are kept constant while the irradiance E is varied.The onset of the laser pulse determines the start point for timing.To establish better clarity for the investigation of the melt formation mechanism and its dependence upon the irradiance, only distributions for the first few hundred nanoseconds are presented.
Section 4.1 reveals, that the irradiance has a strong influence on the temperature evolution of the front surface element.It becomes evident from the temperature distributions along z-axis in Figures 4a-4c, that the irradiance also determines the melt front speed at the Si/Sb interface.At t ≈ 330 ns, the melt depth d m is just a little larger than 5 μm (that is to say, the melt interface just crossed depth z = 195 μm) for the lowest irradiance E = 40 MW/cm 2 .Doubling the irradiance E sets a melt depth of d m ≈ 15 μm and a fourfold irradiance increase Additionally, Figures 4a-4c show that the front side temperature is not thermally coupled to the interface/rear side temperature.In all temperature distributions, the minimum temperature lies close to the melt front emerging from the rear side.This demonstrates that no outgoing heat flow from the heated elements at the rear side reaches the front side.On the other hand, the highest irradiance E = 160 MW/cm 2 result in the undesirable strong temperature increase at the front surface, which abortively stops the simulation because T stop = 3500 K is reached.
Figure 5 presents the mean melt front speed v m for different irradiances E (derived from melt depth vs. time curves, like Fig. 7).The mean melt front speed v m scales almost linearly with the irradiance E. Sections 4.1 and 4.2 explain why the melt front speed v m strongly depends on E. As shown in Figure 3a, the Si/Sb interface element temperature increases rapidly when the laser is switched on during the first 100 ns for all pulse durations with corresponding irradiances alike.However, the extent of the temperature overshoot is larger for higher irradiances, which is apparent from Figure 4a.Hence, the larger local temperature gradient results in a higher melt front speed.

Influence of the beam diameter
By focusing on the on-axis temperature profiles and spot center elements, the problem was kept one-dimensional in the previous subsections.In order to derive the threedimensional information, which is accessible with our numerical model, this section investigates the influence of the laser beam diameter d 0 on the melting mechanism and temperatures at Si/Sb interface (a) as well as Si front surface (b).
Due to the spatial Gaussian beam shape, the melt front inside the silicon is not parallel to the surface.This fact implies, that temperature profiles show in Figures 4a-4c as well as the mean melt front speeds v m from Figure 5 constitute upper boundaries for the off-axis (r > 0) conditions.In this section, the irradiance E = 40 MW/cm 2 is kept constant, to evaluate the implications of the Gaussian beam shape, especially the influence of the beam diameter (defined for the 1/e 2 peak intensity value).
The previous subsections show that the temperature of the front surface element (b) constitutes a limiting factor for the applicability of a specific laser parameter combination (Q p , p d and d 0 ) with corresponding irradiance E. Increasing the beam diameter d 0 , the pulse duration p d and irradiance E all enhance the thermal built-up at the front side.With regard to an application for deep laser doping, only parameter combinations which do not result in a significant temperature increase at the front surface are considered.
Figure 6 shows the maximum surface element (b) temperature at the laser spot center.Due to the constant irradiance E = 40 MW/cm 2 and constant pulse duration p d = 4 μs, the pulse energy Q p scales with increasing beam diameter d 0 (curve for d 0 = 2.5 μm is not shown because the surface element stays at room temperature).While the surface temperature remains below T = 400 K for d 0 ≤ 10 μm, the larger beam diameters result in a steep temperature increase in the surface element (b).The smaller the irradiated area, the higher is the temperature gradient to the surrounding elements.In turn, this leads to increased heat transport in r-direction and better cooling of the center element.Thus, for a given pulse duration and irradiance, the smallest beam diameter results in the lowest front surface element temperature.During the laser pulse (with a temporal "Super-Gaussian" shape), the thermal built-up increases continuously, thus applicable pulse durations decrease with increasing beam diameter.
To readily compare the effect of the different parameters (beam diameter combined with pulse duration) on the radial extension of the melt, we introduce a melt depth-towidth aspect-ratio R MD/MWHMD = m d /m HMD w , which compares the maximum melt depth d m (in z-direction) with the associated melt width at half maximum melt depth m HMD w (in r-direction).This aspect-ratio R MD/MWHMD constitutes a measure for the "shallowness" of the melt pool.Fundamentally, larger beam diameters cause a greater radial melt extent, which leads to a stronger thermal-built up.Smaller beam diameters exhibit a more effective cooling behavior of the on axis element due to larger heat fluxes in r-direction.
Figure 7 shows the melt depths m d in the laser spot center at the wafer's rear side as a function of time.Only parameter combinations for E = 40 MW/cm 2 which do

20104-p7
The European Physical Journal Applied Physics Table 3. Maximum melt depth and melt-depth-to-width aspect-ratio for 10 μm beam diameter.not reach the stopping temperature T stop = 3500 K (see Fig. 6) are plotted.Increasing beam diameters lead to a steeper slope (higher melt front speed) and larger maximum melt depth d m , because of diminished heat loss in r-direction due to enhanced thermal-built up at the rear side.The observation, that the aspect-ratio R MD/MWHMD enlarges with increasing beam diameter, reveals that the melt depth grows faster than the melt width.In fact, this implication affirms the self-amplifying absorption processes, because otherwise the melt width would grow at a similar rate than the melt depth, resulting in constant aspect-ratios R MD/MWHMD .Due to the strong temperature dependence of α Si , the irradiation transmitted to the rear side is absorbed at the melt front.Consequently, no "active heating" of the elements behind the first molten element occurs.Hence, the temperature gradient to the neighbor elements crucially influences the phase of the center element.
Table 3 summarizes maximum melt depths d m with associated aspect-ratios for R MD/MWHMD for constant beam diameter d 0 = 10 μs.The applicable irradiance E decreases with increasing pulse duration p d , due to the continuous thermal built-up at the front surface (see Fig. 6) during the laser pulse.In the upper part of the wafer, the self-amplifying absorption process becomes more and more effective, ultimately resulting in surface melting and evaporation.The high melt depth of d m ≈ 50 μm at E = 80 MW cm −2 for the pulses with p d = 2 μs and 4 μs are impeded by the thermal built-up at the front side, which limits the possible heating period to t ≈ 1.5 μs.Thus, these large melt depths could be accessible with carefully adjusting irradiance and pulse duration.

Energy conservation analysis
Confirmation of energy conservation for the simulation assures confidence in the numerical approach.To satisfy the energy conservation criterion, the incoming irradiation has to equal the sum of absorbed energy inside the  4 comprises the maximum D max and mean D mean deviation from energy conservation (expressed as percentage, relative to the total energy), evaluated over the time steps for each parameter combination.For all parameter combinations, the maximum deviations of D max ≈ 4 − 6.5% only occur at the beginning of the laser pulse.This behavior originates from the approximation to incorporate multiple reflections inside the solid silicon in our model.The mean deviation D mean (including the maximum values at the beginning of each pulse) lies below 1% for all parameter sets, except for small beam diameters d 0 = 2.5 μm with D mean < 1.5%.

Conclusions
The numerical results reveal that the irradiance E determines the melt front speed at the interface temperature.Higher melt front speed combined with longer moving time (longer pulses) enables larger melt depth at the rear side.On the other hand, the irradiance also strongly affects the front side temperature.Unfortunately, the thermal built-up at the front surface likewise increases for longer pulses and greater beam diameter.Ultimately, the front surface temperature constitutes a limiting factor for the application of the method for deep laser doping, because surface alteration by melting or evaporation of the upper part of the wafer is undesirable.
Additionally, this surface heating adversely affects the formation of a continuously molten domain at the rear side.Higher surface temperatures cause stronger absorption in the upper part of the wafer, which in turn leads to lower laser intensity at the melt front.The front side surface temperature is not coupled to the interface/rear side temperature via heat conduction.This means, careful adjustment of the laser parameters allow a large melting depth while avoiding a strong thermal built-up at the front surface.
Our simulation results promise a maximum continuous melt depth of m d = 40 μm in a 200 μm thick wafer for 20104-p8 P.C. Lill and J.R. Köhler: Infrared backwards laser melting of a silicon wafer the parameter combination p d = 1 μs, d 0 = 10 μm and Q p = 31 μJ for a spatial Gaussian and temporal "Super-Gaussian" laser beam with wavelength λ = 1064 nm .
The current mathematical model and computer code for simulation is easily adaptable for other laser wavelengths and absorbing dopant precursor layer materials.In order not to exceed the scope of this work, simulations with longer laser wavelengths, different laser intensity distributions (e.g., spatial "Super-Gaussian" or more sophisticated temporal shapes to access the larger melt depths d m ≈ 50 μm) and lower background temperatures will be included in a future publication.
Further investigations address the possible application of backwards IR laser melting method for a deep laser doping process with an appropriate experimental setup.
P.C. Lill thanks J.H. Werner for his generous support and T. Wurster, K. Ohmer as well as C. Sämann for helpful discussions.

Fig. 1 .
Fig. 1.Schematic view of the backwards IR laser doping process.(a) The IR laser beam passes through the silicon (Si) wafer and is absorbed by the antimony (Sb) layer on the rear side, (b) the Si melt front propagates towards the wafer front surface, (c) Sb atoms diffuse in the liquid Si and are incorporated into the Si lattice upon recrystallization.

Fig. 2 .
Fig. 2. Sample geometry used in the simulation.The radius r runs from the origin r = 0 at the laser spot center up to the right boundary r = W of the computational domain.The depth z runs from the front surface at z = 0 to the end of the thin absorber layer z = L + l at the rear side of the wafer.The interface grid element a and the front surface grid element b are highlighted because their simulated temperature evolution is addressed in section 4.

Fig. 3 .
Fig. 3. (a) Temperature evolution of the silicon/antimony interface element at the center of the laser spot with constant pulse energy Qp = 31 μJ and beam diameter d0 = 10 μm.The colors denote the different pulse durations p d = 0.5 μs, 1 μs, 2 μs and 4 μs.The dashed line indicates the melting temperature of silicon Tm = 1687 K, (b) temperature evolution of the sample front side element at the center of the laser spot.The simulation for p d = 0.5 μs is stopped shortly after the first half of the pulse because the temperature at the front surface element reached Tstop = 3500 K. Our numerical model ceases to be valid at the boiling temperature of silicon T b = 3538 K [14].Temperature distributions along the depth axis are presented in Figure 4.

Fig. 4 .
Fig.4.Temperature distribution at the laser spot center along the depth (z-axis).Beam diameter d0 = 10 μm and pulse duration p d = 1 μs are kept constant while irradiance E varies.The origin of the z-axis lies at the front surface of the sample (the laser impinges at z = 0 μm).The laser is switched on at t = 0 s.

Fig. 5 .
Fig. 5. Mean melt front speed vm at the wafer's rear side in the laser spot center along the depth (z-axis) for constant beam diameter d0 = 10 μm.The mean melt front speed vm scales almost linearly with irradiance E.

Fig. 6 .
Fig. 6.Maximum surface temperature of the front side surface element at the laser spot center with constant irradiance E 40 MW/cm 2 and constant pulse duration p d = 4 μs for varied beam diameter d0.The simulation stops at Tstop = 3500 K, because our numerical model ceases to be valid at the boiling temperature of silicon T b = 3538 K [14].

Fig. 7 .
Fig. 7. Melt depth m d of the wafer's rear side at the laser spot center along the depth (z-axis) for constant irradiance E = 40 MW/cm 2 and constant pulse duration p d = 4 μs.The maximum melt depth increases with increasing beam diameter because of the less effective radial cooling, while increasing R MD/MWHMD show that the melt depth grows faster than the melt width.

Table 2 .
Physical properties of antimony.
The front surface element heats up to Tstop = 3500 K around t = 1.5 μs, making larger melt depth unaccessible. *

Table 4 .
Deviation from energy conservation for constant irradiance E = 40 MW/cm 2 .