Computational fluid dynamics simulations of silicon μ-structure formation during excimer laser modification of silicon-nanoparticle layers

The laser modification of silicon-nanoparticle layers with a nanosecond pulsed excimer laser leads to the self-organized formation of crystalline, μ-cone-shaped silicon structures with good electronic properties, which have allowed the demonstration of their potential for printed flexible electronics. With the current nanosecond laser process, silicon exhibits only short melting times, resulting in a method-defined substrate contact angle, instead of this parameter being defined by the substrate surface energy as expected for equilibrium conditions. This substrate material-independent non-equilibrium contact angle of the silicon melt was experimentally determined in this study to be Θ = 68 ± 10°. To gain deeper insight into the process of the sequential melting and the formation of the silicon μ-cone structures during laser modification, a two-dimensional computational fluid dynamics simulation was performed in COMSOL Multiphysics® solving the Navier–Stokes equation for incompressible fluids. The simulation uses an effective medium approach by applying the conservative level set method to describe the porous silicon-nanoparticle layer. Its sequential melting during the pulsed laser modification is modeled using a newly developed simulation methodology, which uses a time- and depth-dependent dynamic viscosity of the molten silicon. The simulation was carried out for different laser energy densities and verified using scanning electron microscopy images of corresponding laser-modified samples. The simulation results agree well with the experiment and were subsequently used to optimize the laser process.


Introduction
In general, the use of laser technology or other optical thermal post-processing of printed semiconductor thin films is very attractive for the rapidly growing area of printed, flexible electronics, as the required temperature profiles can be applied in an almost digital fashion, meaning local, fast heating and cool down, which results in a low substrate temperature impact. Therefore, here, the focus will be on laser tempering technology, where the careful parameter selection of such variables as wavelength, pulse duration, and intensity enables precise control of the semiconductor modification while minimizing the thermal stress on the underlying substrates [1]. However, due to the high number of variable laser and material parameters influencing the processing result, an empirical development of such a laser process may require a large number of experimental test runs and changes in equipment, motivating theoretical modeling capabilities to describe the laser-material interaction to help predefining a possible processing window for the desired result.
In consequence, this was the motivation for us to recently introduce a two-dimensional (2D) computational fluid dynamics (CFD) simulation [2], to describe the film formation during excimer laser modification of silicon-nanoparticle (Si-NP) thin films [3]. However, the original work by Caninenberg et al [2] assumes an instant melting of the entire Si-NP thin film at the beginning of the laser modification and it simulated the structuring of an initially porous-Si melt for a certain time period. To describe the experiment more accurately, the work at hand extends the CFD model to integrate the sequential melting of the Si-NP layer during the laser modification. For this purpose, the time-dependent melting depth information for each laser pulse is considered, as well as a newly developed simulation methodology for a time-and depth-dependent dynamic viscosity of the silicon melt. This approach treats the material as a liquid, even in its solid state, and sets the dynamic viscosity of Si as a function of a calculated temperature profile/melting depth, either to the actual value for molten Si or to a notional high value, so that the fluid behaves as a solid. The development of this model is described in the following, as well as its experimental validation using the laser modification of Si-NP thin films.

Simulation description
The conducted computational fluid dynamics (CFD) simulations were performed with COMSOL Multiphysics ® (ver. 5.6) using the built-in CFD module, which simulates the fluid flow based on the finite element method (FEM). Among other computational approaches, FEM has proven beneficial to various applications in terms of simulation flexibility and stability [4][5][6][7]. COMSOL was chosen as it is a well-established simulation software for multiphysical problems, combining pre-processing, solving and post-processing in one software. Within the software environment, it is easy to transfer simulation results to other simulations. In our case, we used the simulation results of a 1D temperature simulation of the same laser process, as described in a previous publication [8], as input to the CFD simulation. Furthermore, the contribution at hand is an extension of a model recently developed by us in COMSOL [2], which solves the Navier-Stokes equation for incompressible fluids (equations (1) and (2)), in order to simulate the flow of molten silicon in a nitrogen gas (N 2 ) environment.
The assumption of incompressible fluids is valid for flow velocities below a threshold of 100 m s −1 [9,10], which is applicable for molten Si and N 2 in our case, as here the flow velocity is essentially induced by surface tension and there is no external pressure difference in the system.
Here, ρ represents the density,  u the velocity vector, t the time, p the system pressure, μ the dynamic viscosity and  F combines all volume forces acting on and in the fluids per unit volume, including contributions due to gravity and surface tension. A surface tension coefficient of σ Si melt = 0.72 N m −1 is assumed for molten Si [11].
The 2D model is assembled out of two domains, one which represents the porous Si-NP layer and another domain that consists of nitrogen on top of the Si layer (figure 1). The interface between solid or liquid Si and the substrate (boundary 1 in figure 1) is modeled with a wetted wall boundary condition [12], which requires information about the interaction between substrate and fluid, such as the surface energy of the corresponding materials or the resulting contact angle. A symmetric boundary condition is applied to the boundaries left and right of both domains (boundaries 2), which mirrors the boundary cells and prescribes no penetration and vanishing shear stresses [12]. The upper boundary (3) is set as an open boundary that allows fluids to enter and leave the domain [12]. One point (4) is chosen as a pressure point constraint, which sets the pressure level to atmospheric pressure.
The material properties of both layers are described using an effective medium approach by applying the conservative level set method. For this purpose, a modified depth (y) dependent level set function Φ(y) (equation (3)) is introduced, which allows describing both domains as a single phase with effective material properties X eff (equation (4)), such as an effective density r eff or a dynamic viscosity m eff [13]. The porosity of the Si-NP layer of approximately 65% [2,14] is expressed through the user-defined level set variable of Φ = 0.35, assuming a homogenous mixture of 65% N 2 and 35% Si as the initial value for the Si-NP domain. With this approach no heavily distorted mesh is needed to model the porous Si-NP layer, thus saving computational time.
Here, Φ(y) represents the depth-dependent level set variable, X eff the effective material property calculated from X N 2 and X Si , i.e., the density and dynamic viscosity for either silicon or nitrogen (cf. table 1). The CFD simulation of the laser modification of a 320 nm thick Si-NP layer in N 2 (simulated N 2 thickness 400 nm) is conducted for a section of 10 μm width to reduce the computation time to a reasonable length. This width, however, is large enough to simulate a section with multiple Si μ-structures (diameter < 700 nm), which is sufficient for further analysis. A mesh convergence study resulted in a sufficiently small mesh size of 10 nm for this CFD simulation, leading to a total amount of 72·10 3 elements. The default solver settings in COMSOL were chosen, which iterates to solve the non-linear Navier-Stokes equation system [12]. The direct linear solver PARDISO was used to solve a linearized version of the nonlinear system in each iteration, by applying the Newton method. The direct solver uses Gaussian elimination to solve the linearized matrix system and has been shown to be more robust than iterative solvers [15]. The chosen model settings resulted in a computation time of approx. 105 h. Our original CFD model 2 assumed an instant melting of the entire Si-NP layer at the beginning of the laser modification in a first-order approximation and it simulates the structuring of an initially porous Si-melt for a continuous time period, without considering solidification within the Si-NP layer. For the present CFD model, the sequential melting of the Si-NP layer during pulsed laser modification is considered by implementing the time-dependent melting depth for each laser pulse. This melting characteristic is extracted from a corresponding temperature simulation of the same pulsed laser modification, which solves the heat diffusion equation in a 1D model, as recently described [8]. The temperature simulation considers multiple consecutive laser pulses at a specific substrate location with different intensities, which will result from a scanned Gaussian laser pulse profile. Further, the CFD model in this work uses a time-dependent dynamic viscosity for Si, in order to simulate the depth-dependent melting and solidification of the Si-NP layer. The applied dynamic viscosity methodology has been previously described for the modeling of temperature-dependent melting and solidification of socalled phase change materials in thermal energy storage applications [16][17][18][19]. Such variable viscosity methods treat the material as a fluid, even in the solid state by modifying the viscosity to a notional high value, so that the fluid is forced to behave as a solid [20]. Here, the temperature dependency of the reported methodology is transferred to a time dependency, as the temperature and CFD simulation are decoupled.
For this purpose, a time-(t) and depth-(y) dependent smoothed Heaviside function B(t,y) (equation (5)) is defined for which a certain depth (y) in the simulation domain (i.e., in the Si-NP layer) is compared with the melting depth ML(t) at a certain point in time (t). This function B(t,y) is equal to 1 if the difference of a certain y and ML(t) is above 0, which means the considered y in the Si-NP layer is above the melting depth, i.e., Si is molten. For the case of B(t,y) = 0, the difference between the considered y and ML(t) is below 0, i.e., Si is solid. The function B(t,y) was implemented in COMSOL using the built-in function 'flc2hs' [21], as shown in Table 1. Summary of required material properties for the conducted 2D CFD simulation. Data collection for molten silicon (Si melt ) and nitrogen (N 2 ).
With n being the half-width of the transition zone from liquid to solid (here n = 10 ns). The function B(t,y) was then inserted in another time-and depth-dependent function A(t,y) (equation (7)), which equals zero when the material is liquid and equals 10 8 when the material is solid, and it was used to describe the time-and depth-dependent dynamic viscosity of silicon m ( ) t y , Si (equation (8)): Si Simelt with C and q being constants with assigned values of 10 5 and 10 −3 respectively [18,19], and m Si melt representing the dynamic viscosity of molten Si (cf table 1). Furthermore, an additional volume force  F a (equation (9)) was considered for the Navier-Stokes equation (equations (1) and (2)), which dominates all other forces in the Navier-Stokes equation for the case of A(t) = 10 8 , i.e., the fluid is in solid state, forcing the simulation to the trivial solution of =  u 0, which holds true for a solid medium [18].

Experimental framework for simulation validation
The synthesis of the used Si-NPs, their dispersion in organic solvents, and the related sample preparation are described elsewhere [8,[25][26][27]. A pulsed KrF excimer laser (ATLEX-300SI, ATL Lasertechnik GmbH) emitting at a wavelength of 248 nm was used to laser modify the Si-NP layers. The laser operated with a pulse duration of 6.2 ns at full width half maximum (FWHM), a repetition frequency of 50 Hz, and a scanning speed of 9.3 mm min −1 . The beam optics effectuate a top-hat-Gaussian line profile of 10 mm × 24 μm (FWHM). The effective laser energy density was changed between 0.83-2.08 J cm −2 . The laser modification was conducted within a nitrogen (N 2 ) flushed process chamber in order to avoid oxidation of the Si-melt. The resulting silicon microstructures were analyzed using a JSM 7500F SEM (JEOL Ltd) scanning electron microscopy (SEM). The evaluation of SEM images was done with the software ImageJ (ver. 1.52a).
For the simulation to run, material-dependent information about the liquid wetted wall boundary interaction needs to be available. In our case, this is the contact angle Θ for molten Si on the used substrates. The contact angle dictates which type of wetting occurs and typically can be calculated in thermodynamic equilibrium using Young's equation [28]. According to this equation, an increasing surface energy of the substrate leads to a decreasing equilibrium contact angle, i.e., substrate wetting. However, the equilibrium contact angle is typically reached after spreading times of ∼10 −2 s (non-reactive wetting) [29] up to 10 -104 s (reactive wetting) [30]. Time-dependent wetting of molten Si on sintered silicon nitride (SiN) was evaluated experimentally by Drevet et al [29]. They determined an initial contact angle of Θ 0 = 82 ± 3°, which decreases to its equilibrium value of Θ F = 49 ± 3°in about 150 s. For our case, dealing with ns-laser annealing of Si, the silicon is liquid for less than 130 ns after each laser pulse [8]. Hence, it can be assumed that the silicon solidifies with a non-equilibrium contact angle, independent of the substrate. In consequence, the contact angle is defined by the modification method, i.e., the ns-laser annealing, and not the substrate surface energy, as under equilibrium conditions. To measure this non-equilibrium contact angle, Si-NPs were laser modified on substrates with different surface energies, such as tungsten (W) on titanium (Ti), aluminum (Al), glass, and silicon nitride (SiN). W (thickness 120 nm) on Ti (1 μm) and Al (1 μm) were each sputter deposited on borosilicate glass substrates. A SiN wafer and a cleaned borosilicate glass substrate were chosen as non-metallic surfaces. All surfaces exhibit a comparable low surface roughness (R a < 5 nm), which was measured with a confocal microscope (μsurf custom, Nanofocus AG). The contact angle Θ of the resulting Si μ-cones was analyzed by cross-sectional SEM images, as exemplarily shown in figure 2(a) for a Ti/W coated glass substrate (cf figure S1 in the supporting information (SI) for Si μ-cones on other substrates). The experiment confirms our hypothesis, resulting in a contact angle that is almost independent of the substrate and thus also of the surface energy, averaging at Θ = 68 ± 10°( figure 2(b)). Only for Al, the contact angle is slightly lower. This is attributed to the low melting point of Al (933.4 K) [31], which is lower than the melting point of the Si-NPs (1683 K), and may result in partial Al melting during the laser processing [32,33].

CFD Simulation with melting depth profile
The discussed model is used to describe the film formation of a 320 nm thick Si-NP thin film on a Ti/Wmetallized glass substrate during pulsed excimer laser processing for different laser energy densities. The corresponding time-related melting depth profile for the considered 12 successive laser shots, as extracted from the discussed temperature simulation is illustrated in figure 3 for an energy density of 1.04 J cm −2 . This melting depth profile (ML(t)) is used to model a time-and depth-dependent dynamic viscosity of the Si in the CFD simulation (cf equations (5) -(9)). As the duty cycle between two laser pulses (20 ms) is several orders of magnitude larger than the Si melting time (< 130 ns), the sample returns to thermodynamic equilibrium between the pulses. In consequence, there is no continuous sample heating and resolidified Si can be assumed between the pulses [8,33]. Therefore, to reduce the computational time, the time between melting intervals was  For the selected Si-NP layer thickness of 320 nm and laser energy density of 1.04 J cm −2 , the Si-NPs melt sequentially, and five pulses are required to completely melt the thin film ( figure 3). As the temperature simulation assumes a poly-crystalline resolidification, the layer consists of only poly-crystalline Si (poly-Si) for the following pulses, which leads to a different melting behavior than before due to the significantly different material properties of the Si-NPs and poly-Si in terms of density and thermal conductivity [8]. Consequently, the next laser pulses cause either partial melting (pulses 6, 7, and 11 -12 in figure 3) and once more complete melting (pulses 8 -10). The maximum laser intensity is reached in the 9th pulse. In total, the Si is completely molten for 41.5 ns at the W-interface for this experimental framework.
In figure 4 the results of the related CFD simulation are shown as the volume fraction of N 2 , with dark red being 100% N 2 and dark blue 0% N 2 , i.e., 100% Si. The results are shown for characteristic times in the melting profile marked in figure 3 using the orange dots a-j. A video of the whole CFD simulation can be found in the SI. The porous Si-NP layer is initially considered as a homogenous mixture of 35% Si and 65% N 2 (yellow color in figure 4(a)). After the second pulse (t = 64 ns, figure 4(b)), some Si-agglomeration can be seen within the part of the Si-NP layer which has already been molten (melting depth ∼100 nm). This agglomeration phenomenon is known as the balling effect in material laser processing [34]. The un-molten Si-NP layer below shows a breakup of the initially homogeneously distributed Si/N 2 mixture, resulting in an agglomeration-like state of the Si. We assume that this break-up arises from numerical inaccuracies during the calculation of the material distribution using the level set variable (Φ), as the break-up does not appear if only one material is assigned to the lower domain (i.e., Φ = 1). Even though this phenomenon arises from numerical noise in the CFD simulation, it serves here to model the statistically random pore structure of the actual Si-NP layer, which acts as dewetting centers for liquified silicon and initiates the formation of the anticipated Si μ-structure. The next pulses show a sequential melting of the Si-NP layer, and the surface tension of molten Si leads to an ongoing and more pronounced formation of Si-agglomerates with increasing size (figures 4(c) -(f)). The subsequent pulses 8-10 (figures 4(g) -(i)) exhibit complete melting of the Si. In this process, the contact angle of the Si structures on the W interface converges sequentially towards the experimentally determined value of 68°, as set as the boundary condition for the W interface (cf figure 1). The resulting Si μ-cone size distribution depends on the melting characteristic and the total melting time of the Si, as discussed later. Any subsequent pulses after pulse 10 do not lead to further changes in the Si μ-cone structure ( figure 4(j)). Most formed Si μ-cones exhibit a proper interface to the W-substrate, whereas some show small inclusions of nitrogen, which might deteriorate the electric contact of these Si μ-cones. Further, some isolated, virtually floating circular Si structures (cf figure 4 at 0.15·10 4 nm) can be seen, which are related to the applied simulation methodology. The modeling of the Si solidification by changing the dynamic viscosity of the material in a certain depth to a notional high value (cf simulation description section), basically sets the velocity of the Si flow in this whole region to zero. This may lead to the formation of free-floating Si structures, which are physically not meaningful and requires further optimization in the simulation. However, these isolated events do not challenge the validity of the model, which is outlined in the following.

Experimental validation
To validate the accuracy of our modeling effort, Si μ-structures for different laser energy treatments were compared to our simulation results, as summarized in figure 5. At an energy density of 0.83 J cm −2 , the Si-NP layer melts completely only once and Si is only molten for a duration of 20.5 ns at the interface to W (cf figure S2a in SI for melting depth profile). This short melting time leads to an incomplete Si μ-cone formation and an impure Si/W-interface (figure 5(a).1, cf figure S2b in SI for intermediate steps of the structuring). The simulation result is in good agreement with the SEM cross-sectional image of a comparable sample ( figure 5(a).2). The non-ideal Si/W interface can be seen in form of smaller inhomogeneous Si structures at the bottom of the shown Si μ-cone. The CFD result for the already discussed laser energy density of 1.04 J cm −2 ( figure 5(b).1, as described in more detail above cf figure 3 and figure 4) shows mostly separated Si μ-cones and small Si residuals, both with a proper interface to the W-substrate and some small inclusions of nitrogen. The CFD result agrees well with the SEM image of a Si μ-cone in cross-sectional view ( figure 5(b).2). The highest simulated laser energy density of 2.08 J cm −2 leads to a complete melting of the thin film for multiple times and a total duration of 293.5 ns, for which the Si is molten at the interface to W (cf figure S3a in SI for melting depth profile). This melting characteristic leads to the formation of separated and slightly bigger Si μ-cones with proper interfaces to the W-substrate (figure 5(c).1, cf figure S3b in SI for intermediate steps of the structuring). Only one Si μ-cone shows a nitrogen inclusion and no other small Si residues are present between the Si μ-cones. We can conclude that the shape and size of the simulated (figure 5(c).1) and experimentally (figure 5(c).2) obtained Si μ-structures are in good agreement, which will be discussed further in the following.
In order to do so, SEM top-view images of generated Si μ-cones at laser energy densities of 0.83 J cm −2 , 1.04 J cm −2 , 1.25 J cm −2 , 1.66 J cm −2 , and 2.08 J cm −2 were analyzed with the software ImageJ regarding the diameter distribution and the number density of generated Si μ-structures. A detailed description of the analysis procedure is given in the SI (cf figure S4). An SEM top-view image of Si μ-cones is exemplarily shown in figure 6(a) for a laser energy density of 1.04 J cm −2 , other SEM micrographs can be found in the SI (cf figure S5). The Si μ-cones exhibit a log-normal diameter distribution ( figure 6(b), cf figure S6 for sizes distributions of other laser energy densities), as can be found for many random natural phenomena [35]. The diameter distribution can be fitted with a log-normal function (equation S1 in SI).
The results of the SEM-image analysis, in terms of geometric mean diameter D g and standard deviation STD (cf. equation S2 in SI) of Si μ-cones are compared to the results of the performed CFD simulations (figure 7). At first, the measured D g remains relatively constant at 500 ± 160 nm for a laser energy density range of 0.83 -1.25 J cm −2 and rises for increasing laser energy, up to 690 ± 152 nm at 2.08 J cm −2 . The results of the CFD simulations show the same trend and the determined mean diameters are in the same size range as the measured values. Especially the simulations for laser energy densities which lead to the formation of Si μ-cones with a proper interface to the W-substrate (1.04 J cm −2 and 2.08 J cm −2 , cf. figure 5(b).1 and c.1) are in very good agreement with the measurement.
Another parameter to describe the Si μ-cones, aside from contact angle and diameter, is their height. This parameter, however, can only be roughly determined from the cross-sectional SEM images, as the height of Si μcones not directly at the fracture edge may appear larger or smaller depending on the exact viewing angle and due to the distortion of the perspective. The measured height increases linearly from 314 ± 68 nm at 0.83 J cm −2 to 435 ± 111 nm at 2.08 J cm −2 (cf. figure S7 in SI), which is consistent with the observed increase in diameter. In comparison, the simulated Si μ-cones are clearly smaller in height (174 ± 60 nm at 1.04 J cm −2 and 242 ± 63 nm at 2.08 J cm −2 , cf figure S7 in SI). The highest Si structures (254 ± 54 nm) are simulated for 0.83 J cm −2 , which can be ascribed to the incomplete formation of Si μ-cones, as the simulated Si structures remain on top of only  partly molten Si, increasing their height. The clear difference in height between simulated experimentally obtained Si μ-cones cannot only be attributed to inaccuracies in the height measurement as indicated above. We suggest that the height variation is also due to the fact that the Si μ-cones in the experiment show some variation in the contact angle (Θ = 68 ± 10°, cf figure 2(b)), and some exhibit a slightly elongated shape, which can impact their size and cannot be described with the idealized COMSOL approach.
Finally, the ImageJ evaluation of the SEM images is used to correlate the Si μ-cone number density for the experiment and the simulation ( figure 8). Here, the number of Si μ-cones in the simulated section (10 μm) is projected onto an area of 1 mm 2 . In general, the growth of Si μ-cones (diameter and height) with increasing laser energy is accompanied by a linear decrease of the number density from 1.56·10 6 mm −2 at 0.83 J cm −2 to 0.52·10 6 mm −2 at 2.08 J cm −2 . This can be explained by the fact that the same initial mass of Si is either divided into many but smaller or larger and fewer Si μ-cones, in addition to potential evaporation phenomena. The simulation results show the same decreasing trend and comparable number densities as the experiment (figure 8), which confirms the good agreement between the simulation and experiment.

Conclusion
The self-organized structuring of molten silicon during pulsed excimer laser modification of siliconnanoparticle (Si-NP) layers was described using a computational fluid dynamics simulation and validated experimentally. The presented CFD simulation uses a newly developed methodology for a time-and depthdependent dynamic viscosity of molten silicon, to model the sequential melting of Si-NPs during the pulsed laser modification. This approach treats the material as a fluid, even in solid state, and sets the dynamic viscosity of Si as a function of a previously simulated melting depth, either to the actual value for molten Si or to a notional high value, so that the fluid behaves as a solid in the simulation. As the CFD simulation requires information about the interaction of fluid and substate, the contact angle between the Si melt and the substrate was experimentally determined. The contact angle of Si μ-cones was analyzed by means of cross-sectional SEM images, which revealed a substrate-independent, i.e., surface energy independent, contact angle of Θ = 68 ± 10°. This was interpreted as a non-equilibrium contact angle, as Si is only molten for short periods (< 130 ns), so it solidifies before reaching the equilibrium contact angle. In this case, the contact angle is determined by the modification method, i.e., ns-laser annealing, and not the substrate surface energy as would be expected for equilibrium condition.
The conducted CFD simulations show the formation of separated Si μ-cones with the experimentally determined contact angle of 68°, when Si is molten for a sufficient duration at the substrate-interface, i.e., at adequate laser energies. Lower energies and thus shorter melting times do not lead to the formation of separated Si μ-cones and lower contact angles. These results agree well with experimental observations. For higher laser energies, the Si μ-cones increase in size (diameter and height) and simultaneously separate further from each other, i.e., the number density decreases. This dependency of the Si μ-cone properties on the laser energy density could be well described by the CFD simulation, demonstrating a good agreement.