Atomistic-Continuum Study of an Ultrafast Melting Process Controlled by a Femtosecond Laser-Pulse Train

In femtosecond laser fabrication, the laser-pulse train shows great promise in improving processing efficiency, quality, and precision. This research investigates the influence of pulse number, pulse interval, and pulse energy ratio on the lateral and longitudinal ultrafast melting process using an experiment and the molecular dynamics coupling two-temperature model (MD-TTM model), which incorporates temperature-dependent thermophysical parameters. The comparison of experimental and simulation results under single and double pulses proves the reliability of the MD-TTM model and indicates that as the pulse number increases, the melting threshold at the edge region of the laser spot decreases, resulting in a larger diameter of the melting region in the 2D lateral melting results. Using the same model, the lateral melting results of five pulses are simulated. Moreover, the longitudinal melting results are also predicted, and an increasing pulse number leads to a greater early-stage melting depth in the melting process. In the case of double femtosecond laser pulses, the pulse interval and pulse energy ratio also affect the early-stage melting depth, with the best enhancement observed with a 2 ps interval and a 3:7 energy ratio. However, pulse number, pulse energy ratio, and pulse interval do not affect the final melting depth with the same total energies. The findings mean that the phenomena of melting region can be flexibly manipulated through the laser-pulse train, which is expected to be applied to improve the structural precision and boundary quality.


Introduction
Femtosecond lasers play an important role in micro-and nano-fabrication because of the advantages of high energy intensity and small collateral damage [1][2][3][4][5].The interaction between femtosecond lasers and metallic materials can be divided into two stages [1,6,7].Firstly, upon femtosecond laser irradiation, the electron sub-system absorbs laser energy as a result of a sharp electron temperature increment.The lattice sub-system remains at a low temperature, resulting in thermal nonequilibrium between the electron sub-system and the lattice sub-system [8][9][10][11].Secondly, this nonequilibrium drives electron-phonon coupling energy transport.The whole relaxation process between the two sub-systems lasts several picoseconds and is accompanied by the combined effect of electron heat conduction and electron-phonon coupling [8,9].Eventually, the increment of lattice temperature triggers photo-mechanical and photo-thermal effects, leading to an abrupt phase change [10,[12][13][14].
In recent decades, a novel melting phenomenon known as homogeneous melting has been discovered.Distinguished from the conventional melting process (heterogeneous melting [23]), when subjected to femtosecond laser irradiation, the metal undergoes rapid heating and reaches a highly superheated state.Subsequently, the liquid nucleation process takes place homogeneously, leading to complete melting within a few picoseconds [23,24].
Its mechanism has prompted extensive theoretical, experimental, and simulation investigations [4,8,[24][25][26][27].The transition from the solid phase to the liquid phase of Al was observed via ultrafast electron diffraction (UED) by Siwick et al. [25].Rethfeld et al. proposed the theory of ultrafast thermal melting in laser-excited solids, specifically referring to the phenomenon of homogeneous nucleation [24].The formation mechanism of the two melting phenomena was further analyzed by Ivanov and Zhigilei using the molecular dynamics coupling two-temperature model (MD-TTM model) [8].The transition from heterogeneous to homogeneous melting at higher energy densities is visualized through UED by Mo et al. [4].The influence of electron heat conduction and electronphonon coupling on laser-metal interaction was elucidated by Ivanov et al., and the pulse duration and the laser fluence could affect the above processes [10].In general, under femtosecond laser irradiation, the metals are heated to a sufficiently high temperature where heterogeneous nucleation is suppressed and homogeneous nucleation predominates.Consequently, the system undergoes catastrophic homogeneous nucleation and completes the phase change within a few picoseconds.
The aforementioned research demonstrates that the control of ultrafast melting processes encounters significant challenges attributed to the imperative requirement for achieving electron dynamic regulation.Fortunately, the electron dynamics control (EDC), by the shaping femtosecond laser-pulse train, can realize the effective control of femtosecond laser fabrication at the electron level [1].Significant advancements have been made by precisely designing laser-pulse trains to modify thermophysical properties and subsequently alter the phase-change process [1,28].Li et al. discovered that the pulse train can alter energy transport and associated phase change, enabling control over the size distribution of metal nanoparticles [15].Povarnitsyn et al. unveiled two mechanisms for suppressing ablation irradiated by a double-pulses laser (DP) [29].Förster and Lewis conducted a comprehensive analysis of the mechanisms active during the different stages of the ablation process induced by a DP [13].Englert et al. controlled the ionization processes in high-band-gap materials via temporal asymmetric pulse shapes on intrinsic time and intensity scales [30].The incubation phenomenon during laser ablation with bursts of femtosecond pulses was reported by Gsudiuso et al. [31].Rong et al. found the advantages of femtosecond laserpulse train processing in reducing thermal-mechanical damage compared to single-pulse (SP) processing [17].
Despite the wide range of potential applications for laser-pulse trains, current research on the spatiotemporal shaping femtosecond pulses primarily focuses on the ablation phenomenon.There is limited direct research on controlling the melting process with spatiotemporal-shaping femtosecond pulses, particularly with pulse numbers exceeding two.However, the melting process directly determines the quality, efficiency, and precision of the processing.Roy et al. investigated the influence of Cu nanoparticle size and morphology on sintering temperature and sintering quality.It was found that the sintered product of Cu nanoparticles exhibited outstanding mechanical strength and high electrical conductivity [32].By employing ultrashort laser pulses in assisted powder-bed fusion, Ullsperger et al. discovered a more uniform melting pool shape, along with refined primary Si and eutectic structure during the powder-bed fusion process [33].Furthermore, advancements in first-principles research have shown that the material properties affected by electron temperature in the MD-TTM model are not suitable to be expressed by classical constant or linear equations [17,34,35].
In this work, firstly, the laser pulse is divided 1:1 by the combination of beam splitter and mirrors.The response of Al to the time interval and energy at energies slightly below the ablation threshold is investigated by adjusting the time interval and total energy of the double pulses.Secondly, the thermophysical properties of Al obtained by different methods are compared and applied to the MD-TTM model.After the model validation, the significance of temperature-dependent thermophysical parameters in laser-material interactions during the ultrafast melting process is elucidated through comparative analysis.Subsequently, a simplified method, similar to that used in [22,36], was employed to combine the 1D results to obtain 2D results for both the SP and DP with a 2 ps interval, thereby replicating the observed diameter increase in the melting region in the experimental results and analyzing relevant phenomena.Next, based on the verification of SP and DP simulation, the 2D melting region in Al was predicted under five pulses (FP) by employing the same mode.Furthermore, we also simulated the early-stage melting depth of Al when irradiated by a DP.The influence of different pulse numbers, pulse intervals, and pulse energy ratios on the longitudinal melting (propagation direction) phenomena of Al, including differences in melting depth and melting velocity, was investigated, along with their potential mechanisms.

Experimental Setup
Figure 1 illustrates the schematic diagram of the experimental setup.The femtosecond laser source is a Ti: sapphire chirped pulse-amplification (CPA) system (Spectra-Physics Spitfire Ace, Spitfire Pro-35F1KXP, California, CA, USA), operating at a repetition frequency of 1 kHz, with a full width at half maximum (FWHM) of 100 fs and a center wavelength of 800 nm.The laser pulse is split into two equal-energy pulses using a beam splitter (Thorlabs, Inc., Newton, NJ, USA) in a 1:1 split ratio.These two pulses are combined to form a DP train through the application of mirrors and a beam splitter.The time interval between the double pulses is precisely controlled by a one-dimensional translation stage (the minimum linear displacement is 1 µm, Zolix, Inc., Beijing, China).Subsequently, the pulses are focused onto the sample by a plano-convex lens ( f = 50 mm, Thorlabs, Inc., Newton, NJ, USA).The sample Al (10 mm × 10 mm × 1.0 mm, one side polished, orient <100>) is positioned on a hexapod six-axis positioning system (M-840.5DG,Physik Instrumente, Inc., Auburn, MA, USA).The microstructures are observed via scanning electron microscopy (SEM, Thermo Scientific Apreo 2C, Thermo Fisher Scientific Inc., Waltham, MA, USA).

The MD-TTM Model
The Al melting process under femtosecond laser-pulse train irradiation is simulated by the MD-TTM model [34,35].For the atomistic scale, the MD is used to provide a unique atomistic view of the laser-induced phase-change process.To ensure that the experimental results are reliable, we obtained multiple data by keeping the shutter open and moving the sample along a straight line.The average value of these data was then taken as the experimental result.In addition, we imported SEM images into a MATLAB program and uniformly adjusted their contrast ratio.Through this program, we can obtain the SEM images after processing, which have a clear light-dark dividing line.The dividing line is used as the boundary of the melting region.Finally, the length of the major axis and the minor axis is obtained by fitting the elliptic equation according to the dividing line.

The MD-TTM Model
The Al melting process under femtosecond laser-pulse train irradiation is simulated by the MD-TTM model [34,35].For the atomistic scale, the MD is used to provide a unique atomistic view of the laser-induced phase-change process.
where m i is the mass of atom i with velocity v i , F i is the resultant force acting on atom i, which is calculated from the EAM potential [37,38], and γ i is the friction term, which directly relates to the electron-phonon coupling factor G e−ph , where N is the atom number density, k B is the Boltzmann constant, and ∼ F ran is the random force from the Langevin dynamics.
where ∆t is the timestep and ∼ v i is a random vector.The sum of −γ i v i , and ∼ F ran is in direct proportion to the energy transfer of electron-phonon coupling.
For the continuum scale, the TTM model is used to describe laser energy input, thermal diffusion, and temperature evolution.To establish the coupling between the MD model and the TTM model, the lattice sub-system in the classical TTM model is substituted with the aforementioned MD model.

C e (T e )
∂T e ∂t = ∇[k e (T e )∇T e ] − (T e − T l )G e−ph (T e , T l ) + S(z, t), (4a) where C is the volumetric heat capacity, k is the thermal conductivity, and T is temperature.
The subscripts e and l, represent electrons and lattices, respectively.G e−ph is the electronphonon coupling factor.In this work, the thermophysical parameters (C e , k e , G e−ph ) of Al are determined by T e and T l , which are discussed in Section 3.2.The last term in the right hand of Equation (4a) is the laser source S(z, t), which is a Gaussian beam that decays along the laser irradiation direction.
where i is the pulse number of the current calculation, N laser is the total pulse number, R is the optical reflectivity (0.88 [37]), η i is the proportion of the ith pulse to the total energy, F ab is the total absorbed laser fluence, t is the time of the current calculation, t p is the pulse duration, t i is the pulse interval between the ith and (i − 1)th sub-pulses, L op is the optical penetration depth, which equals the reciprocal of the optical absorption coefficient, L ball is the ballistic transportation length (16 nm), and z is the depth.
It should be emphasized that the above MD-TTM model is only applicable to thermal equilibrium conditions, and the solution of non-thermal equilibrium processes needs to use the hyperbolic two-temperature model (HTTM) proposed by Sobolev [39].Nevertheless, the duration of the laser (100 fs) in this work is much longer than the relaxation time of aluminum [40], in this case, the HTTM degenerates into the TTM [39].Hence, it is justifiable to employ MD-TTM for analyzing the ultrafast melting phenomenon.In the case of ablation simulations, the relatively large [41] lateral size (laser spot radius direction, the top 150 nm surface region consists of 84.2 million atoms) can present more abundant phenomena, such as surface melting, rapid resolidification, and generation of voids and a nanocrystalline layer in the sub-surface region [41,42], while smaller sizes cannot.It is well known that atomistic dimensions are in the length scale of Angstrom (Å) [25].Therefore, simulations to match the actual manufactured size in micrometers would require tens to millions of atoms, which requires huge computational resources.Therefore, an appropriate simplification method is a premise used to achieve the research goal of this work.Fortunately, this problem can be solved by the simplified method proposed by Kumar et al. [36] and applied to the MD-TTM model by Zhou et al. [22].
Figure 2 reveals the whole process by using this method to simplify the 3D manufacturing problem into a 1D computational domain.First of all, the laser energy and melting region have a full-scale axisymmetric feature.A 2D computational domain is directly simplified from the 3D computational domain.Next, the 2D computational domain is divided into a series of 1D computational domains.The position of the 1D computational domain in the 2D computational domain determines the total energy it can obtain.It must be noted that this division will indeed affect the calculations of lateral and longitudinal heat transfer, but at the nanoscale subdivision setting, the temperature difference between the 1D computational domain and the 2D computational domain is less than 2% in Zhou et al.'s work [22].For the 1D computational domain, a similar method was used in our previous work [43].The laser irradiated on the Al surface along the negative -direction.Two vacuum regions [27,43] were added to the up surface and down surface of the material region to ensure that atoms did not escape from the computational domain due to material expan- For the 1D computational domain, a similar method was used in our previous work [43].The laser irradiated on the Al surface along the negative z-direction.Two vacuum regions [27,43] were added to the up surface and down surface of the material region to ensure that atoms did not escape from the computational domain due to material expansion.The whole system was divided into grids with independent T e and T l .The temperature of the material was initially set to room temperature: 298.15 K.The size of the differential grid was 4.05 nm × 4.05 nm × 1.0125 nm.Furthermore, according to the findings in [11,44,45], it has been observed that when the film thickness exceeds 100 nm, any further increase in thickness has a negligible impact on the surface laser-material interaction.This finding enabled us to accurately replicate the experimental results obtained from bulk samples (10 mm) through the simulation of 113.4 nm.The lattice sub-system was created using the classical MD method and equilibrated at 298.15 K for 50 ps [43] in terms of the canonical ensemble (NVT) to complete the initialization and relaxation, which are seen in Figure 2. The grids are colored in red and grey to represent the liquid region and solid region, respectively.The colors also agree with the colors in the atomistic snapshots in Section 3.For the convenience of readers, all the parameters in this work are listed in the Nomenclature.

Computational Flowchart
Benefiting from previous works [34,35,[46][47][48], the classical MD-TTM model has already been integrated into the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [49].On this basis, two new modules were developed and integrated into the LAMMPS program.The details are illustrated in Figure 3.At the beginning of each loop, the program first runs the developed Module 1 named Physical Parameters Calculation.The kinetic energy of atoms was collected to calculate the T l of each grid through Equation ( 6) [49] where i is the total kinetic energy of atoms in the grid.Secondly, the new thermophysical parameters of each grid are updated by solving Equations ( 8)- (10).Eventually, all the thermophysical and optical parameters of Al involved in the MD-TTM model were updated and input into Module 2 (the MD-TTM model) to complete the whole loop.It should be emphasized that the electron temperature and lattice temperature shown in this work are obtained by simulation.
Moreover, because the finite difference method (FDM) was utilized to solve the TTM model, it is necessary to meet the stability criterion.
where F o ∆ is the Fourier number, α = k e C e is the thermal diffusivity, and ∆t and ∆x, ∆y, ∆z represent timestep and spacesteps, respectively.If the conditions are not satisfied, the program will report an error.

Laser-Pulse Train
As shown in Figure 4, The spatiotemporal energy profiles throughout the overall femtosecond laser irradiation process are influenced by various factors, such as the pulse number, pulse interval, and energy ratio.When the simulation reaches the "Determine Pulse Mode and Parameters" process of "Module 2" in Figure 3, it will read the laser parameters and generate the desired spatiotemporal shaping laser-pulse train. =  +  +  0.5 (7) where  is the Fourier number,  = is the thermal diffusivity, and  and , ,  represent timestep and spacesteps, respectively.If the conditions are not satisfied, the program will report an error.

Judgement of Energy Conservation
In addition to the stability criterion, to meet the conservation of energy is essential.A series of additional 1D calculations with different aluminum film thicknesses and pulse settings are implemented.Table 1 shows the calculated results of T l .It should be noted that to facilitate the comparison, all temperatures are the overall average temperature at 50 ps.As seen from Table 1, it can be found that the temperature difference of multiple pulse modes is within 20 K for the cases with the same energy and thickness.Even for the 500 nm results with smaller T l , the error is within 2.5%.This is also shown in the standard deviation σ results.Therefore, the influence of the above errors on the simulation can be ignored.
As shown in Figure 4, The spatiotemporal energy profiles throughout the overall femtosecond laser irradiation process are influenced by various factors, such as the pulse number, pulse interval, and energy ratio.When the simulation reaches the "Determine Pulse Mode and Parameters" process of "Module 2" in Figure 3, it will read the laser parameters and generate the desired spatiotemporal shaping laser-pulse train.

Judgement of Energy Conservation
In addition to the stability criterion, to meet the conservation of energy is essential.A series of additional 1D calculations with different aluminum film thicknesses and pulse settings are implemented.Table 1 shows the calculated results of  .It should be noted that to facilitate the comparison, all temperatures are the overall average temperature at 50 ps.As seen from Table 1, it can be found that the temperature difference of multiple pulse modes is within 20 K for the cases with the same energy and thickness.Even for the 500 nm results with smaller  , the error is within 2.5%.This is also shown in the standard deviation  results.Therefore, the influence of the above errors on the simulation can be ignored.

Results and Discussion
Firstly, the response of Al to the time intervals and energies strictly below the ablation threshold is explored by adjusting the time interval and the total energy of the double pulses (Section 3.1).Secondly, multiple sets of thermophysical parameters obtained from different experiments and simulations are collected and compared.The impact of the laser-pulse number on the early stage of the melting process (first 50 ps) is conclusively demonstrated by the MD-TTM model incorporating temperature-dependent thermophysical parameters (Section 3.2).Thirdly, the 1D simulation results are assembled to obtain 2D results for an SP and a DP with a 2 ps interval by a simplified method proposed by Kumar et al., enabling an analysis of the lateral differences in the melting phenomenon [36].Additionally, the same method is employed to predict the results for FP (Section 3.3).Finally, the investigation explores the impact of pulse interval and pulse energy ratio on the longitudinal melting phenomenon in Al by the same model (Section 3.4).

Experimental Results of Al Response to Pulse Interval and Energy
The sample Al is irradiated by an SP and DPs with 1 ps, 2 ps, and 3 ps intervals, respectively, and varying energies.In the case of the DP, the initial laser-pulse energy is divided into two pulses of equal energy.
Figure 5 illustrates SEM images of Al irradiated by SPs with varying energies.The size of the melting region increases as the pulse energy increases.

Experimental Results of Al Response to Pulse Interval and Energy
The sample Al is irradiated by an SP and DPs with 1 ps, 2 ps, and 3 ps intervals, respectively, and varying energies.In the case of the DP, the initial laser-pulse energy is divided into two pulses of equal energy.
Figure 5 illustrates SEM images of Al irradiated by SPs with varying energies.The size of the melting region increases as the pulse energy increases.Figure 6 shows the SEM images of the melting region of four settings of an SP and DPs with 1 ps, 2 ps, and 3 ps intervals at a total fluence of 0.63 J/cm .In terms of the size of the melting region, both DPs with 2 ps and 3 ps intervals demonstrate larger melting compared to the SP (SP: 8.2 μm, the DP with a 2 ps interval: 8.7 μm, the DP with a 3 ps interval: 8.5 μm,).Interestingly, the melting region diameter is the smallest for the DP with a 1 ps interval (7.6 μm).This phenomenon may be attributed to electron excitation, as reported in [50][51][52].Although these factors may result in reduced energy absorption and subsequently a smaller melting diameter, SEM results for the SP and DPs with 2 ps and 3 ps intervals confirm that multiple pulses enhance melting diameter, even when employing the same or lower laser energy.In addition, when comparing Figure 6c,d, the melting region is smaller for a DP with a 3 ps interval than for a DP with a 2 ps interval.Figure 6 shows the SEM images of the melting region of four settings of an SP and DPs with 1 ps, 2 ps, and 3 ps intervals at a total fluence of 0.63 J/cm 2 .In terms of the size of the melting region, both DPs with 2 ps and 3 ps intervals demonstrate larger melting compared to the SP (SP: 8.2 µm, the DP with a 2 ps interval: 8.7 µm, the DP with a 3 ps interval: 8.5 µm,).Interestingly, the melting region diameter is the smallest for the DP with a 1 ps interval (7.6 µm).This phenomenon may be attributed to electron excitation, as reported in [50][51][52].Although these factors may result in reduced energy absorption and subsequently a smaller melting diameter, SEM results for the SP and DPs with 2 ps and 3 ps intervals confirm that multiple pulses enhance melting diameter, even when employing the same or lower laser energy.In addition, when comparing Figure 6c,d, the melting region is smaller for a DP with a 3 ps interval than for a DP with a 2 ps interval.
The sample Al is irradiated by an SP and DPs with 1 ps, 2 ps, and 3 ps intervals, respectively, and varying energies.In the case of the DP, the initial laser-pulse energy is divided into two pulses of equal energy.
Figure 5 illustrates SEM images of Al irradiated by SPs with varying energies.The size of the melting region increases as the pulse energy increases.Figure 6 shows the SEM images of the melting region of four settings of an SP and DPs with 1 ps, 2 ps, and 3 ps intervals at a total fluence of 0.63 J/cm .In terms of the size of the melting region, both DPs with 2 ps and 3 ps intervals demonstrate larger melting compared to the SP (SP: 8.2 μm, the DP with a 2 ps interval: 8.7 μm, the DP with a 3 ps interval: 8.5 μm,).Interestingly, the melting region diameter is the smallest for the DP with a 1 ps interval (7.6 μm).This phenomenon may be attributed to electron excitation, as reported in [50][51][52].Although these factors may result in reduced energy absorption and subsequently a smaller melting diameter, SEM results for the SP and DPs with 2 ps and 3 ps intervals confirm that multiple pulses enhance melting diameter, even when employing the same or lower laser energy.In addition, when comparing Figure 6c,d, the melting region is smaller for a DP with a 3 ps interval than for a DP with a 2 ps interval.In terms of the morphology of the melting region, the melting regions for an SP and a DP with a 1 ps interval exhibit high surface smoothness.Voids appear near the center of the melting region for DPs with 2 ps and 3 ps intervals.Specifically, the DP with a 3 ps interval shows more concentrated voids.In addition, in multiple pulse experiments, the areas of melting exhibit an elliptical shape, with the long axis corresponding to the polarization direction.Moreover, as the number of pulses increases, the long axis becomes relatively longer compared to the short axis.This phenomenon is widely seen in other multi-pulse studies [50,53], and to avoid the influence of polarization, we choose the short axis as the melting diameter.
In conclusion, the interval of the DPs affects the diameters and the surface morphology of the melting region.The diameter of the melting region generated by a DP with a 2 ps interval is the biggest with the largest void region.This will be discussed in detail in Sections 3.3 and 3.4.

Comparison of the Temperature-Dependent Thermophysical Parameters and the Effect of Pulse Number on Simulation Results
Figure 7 presents the thermophysical parameters of Al obtained through various methods.The solid, dashed, and dotted lines represent the results obtained from the full-run quantum treatment (Jiang et al. [54]), classical simplified expression (Ideal Electron Gas and Simplified Expression [55]), empirical formula (Chen et al. [46]), and ab initio QM (Zhigilei et al. [56], Rong et al. [17], and Bevillon et al. [57]), respectively.It is evident that C e is influenced by T e in Figure 7a.As T e increases, C e exhibits initial linear growth and then tends to approach a constant value.Figure 7b illustrates the different k e .In contrast to the previous parameter, k e is affected by the lattice temperature, which decreases as the lattice temperature increases.Regarding G e−ph , shown in Figure 7c, it is worth noting that the precise determination of G e−ph remains an unresolved issue.Based on the current findings, it can be confirmed that the G e−ph of Al exhibits an increasing trend with rising electron temperature.In summary, in this study, the temperature-dependent  is determined through the polynomial fitting.By taking  into account with the empirical formula, the temperature-dependent  and  are derived.The detailed expression is shown below: 1.For  In summary, in this study, the temperature-dependent C e is determined through the polynomial fitting.By taking T l into account with the empirical formula, the temperaturedependent k e and G e−ph are derived.The detailed expression is shown below: 1. For For 3.
Figure 8 shows the melting process of the first 90 nm from the surface of Al after irradiation with different laser-pulse settings, with a total fluence of 350 J/m 2 .Figure 8a-c represents the classical thermophysical parameters group, while Figure 8d-f corresponds to the temperature-dependent thermophysical parameters groups.During the 0 − 20 ps melting process in both groups, Al's surface undergoes a relatively rapid melting stage.Subsequently, from 20 ps to 50 ps, the melting region expands at a velocity lower than the speed of sound, while a clear melting region boundary appears.This two-stage process exhibits similarities to the phenomenon observed in our previous work [43].The appearance of obvious melting is delayed with the increase in the number of pulses.This delay is attributed to the longer thermal time of Al and the lower energy carried by individual laser pulses in the multi-pulse mode.Taking FP as an example, the total time required for the laser energy to be absorbed equals the sum of the total pulse intervals (4 × t interval ) and the time for a Gaussian beam (6 × t p ) to be fully absorbed, as shown in Figure 4-that is, 4.6 ps.Interestingly, the more obvious differences are evident in the visible liquid-phase region boundaries at 50 ps (black dot lines).In the case of the left group in Figure 8a-c, the impact of the change in pulse number is considered negligible.The variation in pulse numbers has minimal impact on the position of the melting region boundaries and the distribution of the liquidus point.However, for the temperature-dependent thermophysical parameter group in Figure 8d-f, more energy is transferred to a deeper position.Moreover, after the appearance of the boundary of the liquid phase region, there are still a large number of nucleation points below it, with liquid atoms comprising over 90% of the total number of atoms in some grids.
The spatial distribution of electron and lattice temperatures at different time points is illustrated in Figure 9.In each sub-figure, the solid lines in the temperature-dependent parameter group exhibit a smoother lattice temperature change with increasing depth compared to the dashed lines in the classical parameter group.This distinction arises from the increased electron thermal conductivity and electron-phonon coupling factor at higher electron temperatures, resulting in enhanced energy transfer towards the deeper region.For electron temperature in Figure 9a-c, as the pulse number increases, the peak time of electron temperature is delayed (SP¯0.5 ps, DP-0.8 ps, FP-4.6 ps), and the peak temperature of the surface is significantly reduced (classical parameters: SP-11, 538.10 K, DP¯8321.62 K, FP¯5052.08 K, temperature-dependent parameters: SP-9905.54K, DP¯7665.19 ps, FP¯5059.64 K).Regarding the lattice temperature in Figure 9d-f, an increment in pulse number leads to a higher lattice temperature.It is worth noting that the application of temperature-dependent parameters will increase the lattice temperature at a depth greater than 30 nm, which is consistent with previous results [17,43].This phenomenon is also the reason for the observed enlargement of the melting region shown in Figure 8.
required for the laser energy to be absorbed equals the sum of the total pulse intervals (4 × ) and the time for a Gaussian beam (6 ×  ) to be fully absorbed, as shown in Figure 4-that is, 4.6 ps.Interestingly, the more obvious differences are evident in the visible liquid-phase region boundaries at 50 ps (black dot lines).In the case of the left group in Figure 8a-c, the impact of the change in pulse number is considered negligible.The variation in pulse numbers has minimal impact on the position of the melting region boundaries and the distribution of the liquidus point.However, for the temperature-dependent thermophysical parameter group in Figure 8d-f, more energy is transferred to a deeper position.Moreover, after the appearance of the boundary of the liquid phase region, there are still a large number of nucleation points below it, with liquid atoms comprising over 90% of the total number of atoms in some grids.The spatial distribution of electron and lattice temperatures at different time points is illustrated in Figure 9.In each sub-figure, the solid lines in the temperature-dependent parameter group exhibit a smoother lattice temperature change with increasing depth compared to the dashed lines in the classical parameter group.This distinction arises from the increased electron thermal conductivity and electron-phonon coupling factor at higher electron temperatures, resulting in enhanced energy transfer towards the deeper region.For electron temperature in Figure 9a-c, as the pulse number increases, the peak time of electron temperature is delayed (SP-0.5 ps, DP-0.8 ps, FP-4.6 ps), and the peak temperature of the surface is significantly reduced (classical parameters: SP-11,538.10K, DP -8321.62K , FP -5052.08K , temperature-dependent parameters: SP-9905.54K , DP-7665.19 ps, FP-5059.64 K).Regarding the lattice temperature in Figure 9d-f, an increment in pulse number leads to a higher lattice temperature.It is worth noting that the application of temperature-dependent parameters will increase the lattice temperature at a depth greater than 30 nm, which is consistent with previous results [17,43].This phenomenon is also the reason for the observed enlargement of the melting region shown in Figure 8.In conclusion, the introduction of temperature-dependent thermophysical parameters can effectively reflect the phenomenon that multiple pulses increase the melting depth (first 50 ps).The increase in pulse number can significantly reduce the surface electron temperature.This causes the lattice sub-system to obtain more energy, thus increasing the lattice temperature.It is worth emphasizing that the enhancement effect of multiple pulses only occurs when the pulse energy is far below the ablation threshold (640 J/m ).When the pulse energy approaches the ablation threshold, the melting depth will be affected by the material size, surface expansion, and ablation phenomena and shows no increase or even a decrease [13,59].Furthermore, the mentioned enhancement effect only occurs in the early stage.Nevertheless, the current results are sufficient to demonstrate that the MD-TTM model incorporating temperature-dependent thermophysical parameters can effectively achieve the control mechanism of Al's thermophysical properties by pulse train.In conclusion, the introduction of temperature-dependent thermophysical parameters can effectively reflect the phenomenon that multiple pulses increase the melting depth (first 50 ps).The increase in pulse number can significantly reduce the surface electron temperature.This causes the lattice sub-system to obtain more energy, thus increasing the lattice temperature.It is worth emphasizing that the enhancement effect of multiple pulses only occurs when the pulse energy is far below the ablation threshold (640 J/m 2 ).When the pulse energy approaches the ablation threshold, the melting depth will be affected by the material size, surface expansion, and ablation phenomena and shows no increase or even a decrease [13,59].Furthermore, the mentioned enhancement effect only occurs in the early stage.Nevertheless, the current results are sufficient to demonstrate that the MD-TTM model incorporating temperature-dependent thermophysical parameters can effectively achieve the control mechanism of Al's thermophysical properties by pulse train.

The Influence of Multi-Pulse on the Phenomenon of Lateral Melting
Figure 10a,b depicts the morphology of the 2D melting region of the Al surface at 50 ps, under the same laser energy, for both the SP and DP with a 2 ps interval, respectively.Among them, the left is the simulation results, and the right is the experimental results.To facilitate a direct comparison between the simulated results and the experimental results presented in Section 3.1, the laser energy at the center of the simulation is set to 440 J/m 2 , which is below the ablation threshold (640 J/m 2 [37,38]).

The Influence of Multi-Pulse on the Phenomenon of Lateral Melting
Figure 10a,b depicts the morphology of the 2D melting region of the Al surface at 50 ps, under the same laser energy, for both the SP and DP with a 2 ps interval, respectively.Among them, the left is the simulation results, and the right is the experimental results.To facilitate a direct comparison between the simulated results and the experimental results presented in Section 3.1, the laser energy at the center of the simulation is set to 440 J/m , which is below the ablation threshold (640 J/m [37,38]).In the results of the SP, it is observed that melting occurs within ~5.1 μm of the farthest position from the center of the laser spot, and the fluence at the position is 160 J/m , representing the melting threshold for the SP.In the results of the DPs, this position expands to ~5.4 μm, with a fluence of 140 J/m .Combining the results from Section 3.2, it is found that DPs can reduce the melting threshold compared to the SP, thereby enlarging the diameter of the melting region.
According to the above results, we speculate the cause of the void.During the subsequent development of the melting process, the atoms with more kinetic energy in the central region flow toward the lower liquid level in the edge region.As the liquid melting region resolidifies into a solid state, the surface expansion decreases, resulting in the formation of irregular voids.The larger diameters of the melting region observed for the DP indicate a higher kinetic energy of Al at the surface during liquefaction.This phenomenon may be the cause of the generation of irregular voids near the center for a DP with a 2 ps interval, as discussed in Section 3.1.Combining this with the conclusion in Section 3.4 (that the early-stage melting depth is larger for a DP with a 2 ps interval compared to a 3 ps interval), the phenomenon of a more concentrated distribution of voids in the SEM image of the DP with a 3 ps interval in Figure 6d may be caused by the lower kinetic energy of liquid atoms in the melting region.
Figure 11 presents the predicted results for FP under the same conditions discussed in this section.The results of FP exhibit larger melting diameters and depths, consistent with the previous discussions (Section 3.2).In the results of the SP, it is observed that melting occurs within ∼ 5.1 µm of the farthest position from the center of the laser spot, and the fluence at the position is 160 J/m 2 , representing the melting threshold for the SP.In the results of the DPs, this position expands to ∼ 5.4 µm, with a fluence of 140 J/m 2 .Combining the results from Section 3.2, it is found that DPs can reduce the melting threshold compared to the SP, thereby enlarging the diameter of the melting region.
According to the above results, we speculate the cause of the void.During the subsequent development of the melting process, the atoms with more kinetic energy in the central region flow toward the lower liquid level in the edge region.As the liquid melting region resolidifies into a solid state, the surface expansion decreases, resulting in the formation of irregular voids.The larger diameters of the melting region observed for the DP indicate a higher kinetic energy of Al at the surface during liquefaction.This phenomenon may be the cause of the generation of irregular voids near the center for a DP with a 2 ps interval, as discussed in Section 3.1.Combining this with the conclusion in Section 3.4 (that the early-stage melting depth is larger for a DP with a 2 ps interval compared to a 3 ps interval), the phenomenon of a more concentrated distribution of voids in the SEM image of the DP with a 3 ps interval in Figure 6d may be caused by the lower kinetic energy of liquid atoms in the melting region.
Figure 11 presents the predicted results for FP under the same conditions discussed in this section.The results of FP exhibit larger melting diameters and depths, consistent with the previous discussions (Section 3.2).
Such results mean that it is possible to achieve an adjustment of the melting phenomenon through reasonable pulse design (pulse number, total pulse energy, pulse interval, and pulse energy ratio).Such results mean that it is possible to achieve an adjustment of the melting phenomenon through reasonable pulse design (pulse number, total pulse energy, pulse interval, and pulse energy ratio).

The Influence of Pulse Interval and Pulse Energy Ratio on the Longitudinal Melting Phenomenon
Due to the lack of experimental observation methods for longitudinal melting phenomena, this section exclusively relies on simulation to predict and analyze the results of longitudinal ultrafast melting.Figure 12 illustrates the longitudinal evolution of the melting region during the first 50 ps, irradiated by both the SP and DPs with varying pulse intervals.Consistent with the description in Section 3.2, Al mainly undergoes a rapid melting stage dominated by homogeneous melting and a slow melting stage dominated by heterogeneous melting [43].MD snapshots at 10 ps in Figure 12a and 50 ps in Figure 12b are presented to describe the evolution of melting depth at different melting stages.Furthermore, the average velocities for the five stages are plotted in Figure 12c.

The Influence of Pulse Interval and Pulse Energy Ratio on the Longitudinal Melting Phenomenon
Due to the lack of experimental observation methods for longitudinal melting phenomena, this section exclusively relies on simulation to predict and analyze the results of longitudinal ultrafast melting.Figure 12 illustrates the longitudinal evolution of the melting region during the first 50 ps, irradiated by both the SP and DPs with varying pulse intervals.Consistent with the description in Section 3.2, Al mainly undergoes a rapid melting stage dominated by homogeneous melting and a slow melting stage dominated by heterogeneous melting [43].MD snapshots at 10 ps in Figure 12a and 50 ps in Figure 12b are presented to describe the evolution of melting depth at different melting stages.Furthermore, the average velocities for the five stages are plotted in Figure 12c.
During the initial 0 − 10 ps, the melting velocities for all pulse intervals and pulse numbers are significantly higher compared to subsequent stages, necessitating the use of the right y-axis for distinction, as shown in Figure 12c.At 10 ps, there is little difference in melting depth between the SP and DPs, whereas the influence of pulse intervals becomes more pronounced.The maximum melting velocity is observed at the DP with a 0.5 ps interval (5379.6 m/s).As the pulse interval increases, there is an overall decreasing trend in the average velocity during the rapid melting stage.Surprisingly, at a 2 ps pulse interval, the melting velocity is higher than at a 1.5 ps pulse interval.As the melting process further evolves, between 10 ps and 40 ps, the melting velocities decrease to below the speed of sound.The depth differences resulting from pulse intervals gradually decrease.
To provide a more detailed analysis of the mechanisms underlying the variations in melting velocity, Figure 13 shows contour plots of the evolution of electron temperature and lattice temperature with different pulse settings: SP, DP with a 2 ps interval, and DP with a 3 ps interval.
For the SP, the electron temperature at the Al surface easily reaches peak value (13, 914.11K).However, due to the faster electron heat conduction compared to the lattice heat conduction, thermal rapidly diffuses from the surface "hot" electron region to the deep "cold" electron region, reducing the electron temperature, as seen in Figure 13a.Consequently, the thermalization time of the lattice is limited, leading to a smaller range of high lattice temperature regions, indicated by the red color in Figure 13d.During the initial 0 − 10 ps, the melting velocities for all pulse intervals and pulse numbers are significantly higher compared to subsequent stages, necessitating the use of the right y-axis for distinction, as shown in Figure 12c.At 10 ps, there is little difference in melting depth between the SP and DPs, whereas the influence of pulse intervals becomes more pronounced.The maximum melting velocity is observed at the DP with a 0.5 ps interval (5379.6 m/s).As the pulse interval increases, there is an overall decreasing trend in the average velocity during the rapid melting stage.Surprisingly, at a 2 ps pulse interval, the melting velocity is higher than at a 1.5 ps pulse interval.As the melting process further evolves, between 10 ps and 40 ps, the melting velocities decrease to below the speed of sound.The depth differences resulting from pulse intervals gradually decrease.
To provide a more detailed analysis of the mechanisms underlying the variations in melting velocity, Figure 13 shows contour plots of the evolution of electron temperature and lattice temperature with different pulse settings: SP, DP with a 2 ps interval, and DP with a 3 ps interval.During the initial 0 − 10 ps, the melting velocities for all pulse intervals and pulse numbers are significantly higher compared to subsequent stages, necessitating the use of the right y-axis for distinction, as shown in Figure 12c.At 10 ps, there is little difference in melting depth between the SP and DPs, whereas the influence of pulse intervals becomes more pronounced.The maximum melting velocity is observed at the DP with a 0.5 ps interval (5379.6 m/s).As the pulse interval increases, there is an overall decreasing trend in the average velocity during the rapid melting stage.Surprisingly, at a 2 ps pulse interval, the melting velocity is higher than at a 1.5 ps pulse interval.As the melting process further evolves, between 10 ps and 40 ps, the melting velocities decrease to below the speed of sound.The depth differences resulting from pulse intervals gradually decrease.
To provide a more detailed analysis of the mechanisms underlying the variations in melting velocity, Figure 13 shows contour plots of the evolution of electron temperature and lattice temperature with different pulse settings: SP, DP with a 2 ps interval, and DP with a 3 ps interval.For the DP with a 2 ps interval, although the peak electron temperature generated by the first laser pulse decreases compared to the SP (10, 811.35 K), due to the irradiation from the second laser pulse, the electron keeps a high temperature during the initial 4 ps, as shown in Figure 13b.This increases the lattice thermalization time.The lattice sub-system gains more energy, forming a larger range of high-temperature regions, indicated by the red and orange colors in Figure 13e.
For the DP with a 3 ps interval, the increasing pulse interval between the two subpulses allows for sufficient relaxation between the hot electron sub-system and the cold lattice sub-system.When the time the second sub-pulse arrives, the electron temperature reduces to 1212.68 K, and the thermophysical parameter changes caused by the first subpulse have less effect on the absorption of the second sub-pulse.In other words, as the pulse interval continues to increase, the melting phenomena of the two sub-pulses become more independent, leading to a decreasing trend in the surface melting phenomenon.
Furthermore, as indicated in Section 3.2, electron thermal conductivity and the electronphonon coupling factor are positively correlated with electron temperature.Combining the conclusions of [60], the effects of electron thermal conductivity and electron-phonon coupling factor are contradictory on lattice heating.The increase in electron thermal conductivity weakens the lattice heating, and a smaller melting depth is obtained.At the same time, the increase in the electron-phonon coupling factor enhances the lattice heating, and a larger melting depth is obtained.The peak velocity at 2 ps implies that the influence of changes in the electron-phonon coupling factor on lattice thermalization is greater than that of electron thermal conductivity.This facilitates the energy absorption of the lattice sub-system from the second sub-pulse, resulting in a greater melting velocity and melting depth.However, when the pulse interval is greater than 2 ps, during the irradiation of the second pulse, both the electron and lattice have sufficient time to relax and conduct heat deeper into the Al, leading to a "cooler" electron temperature at the surface.Consequently, the melting velocity decreases again, as shown in Figure 13c,f.Figure 14 illustrates the impact of the pulse energy ratio of sub-pulses for DPs on the melting process.During the initial 10 ps of the rapid melting stage, compared to the 1:1 and decreasing energy train, an increasing energy train results in a greater amplification of melting velocities, particularly evident at a 2 ps interval with a 3:7 energy ratio configuration.The simulation results in this section can be summarized as follows: (i) For the rapid melting stage (0 − 10 ps), the melting depth is influenced by the pulse interval and pulse energy ratio.A DP with a 2 ps interval and a pulse energy ratio of 3:7 exhibits the maximum melting depth.(ii) For the whole process of melting, it is observed that as the simulation time increases, the differences in melting depth gradually decrease.Predictably, due to the consistent total input energy, the final melting depth is independent of the pulse The simulation results in this section can be summarized as follows: (i) For the rapid melting stage (0 − 10 ps), the melting depth is influenced by the pulse interval and pulse energy ratio.A DP with a 2 ps interval and a pulse energy ratio of 3:7 exhibits the maximum melting depth.(ii) For the whole process of melting, it is observed that as the simulation time increases, the differences in melting depth gradually decrease.Predictably, due to the consistent total input energy, the final melting depth is independent of the pulse interval and energy ratio.This corresponds to the results presented in [13,59,61].

Conclusions
This study uses the femtosecond laser-pulse train melting experiment and the molecular dynamics coupling two-temperature model (MD-TTM model) incorporating temperaturedependent thermophysical parameters to investigate ultrafast melting processes of Al irradiated by different laser-pulse settings.The investigation encompasses variations in pulse number, pulse interval, and pulse energy ratio to explore both lateral and longitudinal melting processes.The research findings indicate that: (i) the MD-TTM model with temperature-dependent property parameters can effectively reflect the changes in Al's thermophysical properties irradiated by the femtosecond laser train.This model provides theoretical support for the results of femtosecond laser multi-pulse experiments.(ii) As the pulse number increases, the melting threshold at the edge region of the laser spot decreases, resulting in an increase in the melting region diameter.(iii) At the same laser energy (much less than the ablation threshold), the melting depth of Al in the early stage increases with the increase in pulse number.Specifically, melting depth and melting velocity in the early stage reach the maximum for a DP with an increasing energy ratio (3:7) and a 2 ps pulse interval.The pulse number, pulse interval, and energy ratio only influence melting depth and melting velocity in the early stage.This confirms that the independent control of melting phenomena can be achieved through pulse train.This will contribute to the improvement of processing quality in additive and subtractive manufacturing applications using femtosecond laser technology.

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

Materials 2024 , 22 Figure 1 .
Figure 1.Schematic diagram of the experimental setup of double-pulse laser processing.

Figure 1 .
Figure 1.Schematic diagram of the experimental setup of double-pulse laser processing.

Figure 2 .
Figure 2. Simplified dimensional schematic diagram of the Al melting process under laser irradiation and the design of 1D equivalent computational domain.

Figure 2 .
Figure 2. Simplified dimensional schematic diagram of the Al melting process under laser irradiation and the design of 1D equivalent computational domain.

Figure 3 .
Figure 3. Flowchart of the MD-TTM model incorporating dynamic thermophysical parameters.

Figure 3 .
Figure 3. Flowchart of the MD-TTM model incorporating dynamic thermophysical parameters.

Figure 4 .
Figure 4. Schematic diagram of the laser-pulse train.

Figure 4 .
Figure 4. Schematic diagram of the laser-pulse train.

Figure 6 .
Figure 6.SEM images of melting region for (a) an SP and DPs with pulse intervals of (b) 1 ps, (c) 2 ps, and (d) 3 ps, at a total fluence of 0.63 J/cm 2 .

Figure 7 .
Figure 7. Temperature-dependent thermophysical parameters of Al.Comparison of (a) electron heat capacity C e , (b) electron thermal conductivity k e , and (c) the electron-phonon coupling factor G e−ph obtained by different calculation methods [17,46,54-57].

Figure 8 .
Figure 8. Atomic snapshots of the melting process on the Al surface region of 90 nm, the left is the classical thermophysical parameters group for (a) SP, (b) DP, and (c) FP, while the right is the temperature-dependent thermophysical parameters group for (d) SP, (e) DP, and (f) FP.The black dotted line represents the boundaries between the solid and liquid phases.

Figure 8 .
Figure 8. Atomic snapshots of the melting process on the Al surface region of 90 nm, the left is the classical thermophysical parameters group for (a) SP, (b) DP, and (c) FP, while the right is the temperature-dependent thermophysical parameters group for (d) SP, (e) DP, and (f) FP.The black dotted line represents the boundaries between the solid and liquid phases.

Figure 9 .
Figure 9. Spatial distributions of (a-c) electron temperature  , and (d-f) lattice temperature  at different pulse numbers for (a,d) SP, (b,e) DP, and (c,f) FP.Dashed lines represent the classical thermophysical parameters group results and the solid lines represent the temperature-dependent thermophysical parameters group results.

Figure 9 .
Figure 9. Spatial distributions of (a-c) electron temperature T e , and (d-f) lattice temperature T l at different pulse numbers for (a,d) SP, (b,e) DP, and (c,f) FP.Dashed lines represent the classical thermophysical parameters group results and the solid lines represent the temperature-dependent thermophysical parameters group results.

Figure 10 .
Figure 10.The morphology of the 2D melting region (50 ps) of the Al surface at (a) SP and (b) DP with a 2 ps interval, the left is the simulation results, and the right is the experimental results.

Figure 10 .
Figure 10.The morphology of the 2D melting region (50 ps) of the Al surface at (a) SP and (b) DP with a 2 ps interval, the left is the simulation results, and the right is the experimental results.

Materials 2024 , 22 Figure 11 .
Figure 11.The morphology of the 2D melting region (50 ps) of the Al surface at FP with a 2 ps interval.

Figure 11 .
Figure 11.The morphology of the 2D melting region (50 ps) of the Al surface at FP with a 2 ps interval.

Materials 2024 , 22 Figure 12 .
Figure 12.The longitudinal evolution of the melting region during the first 50 ps under both SP and DPs with varying pulse intervals.(a) atomic snapshots at 10 ps, (b) atomic snapshots at 50 ps, and (c) melting speed in different stages.

Figure 12 .
Figure 12.The longitudinal evolution of the melting region during the first 50 ps under both SP and DPs with varying pulse intervals.(a) atomic snapshots at 10 ps, (b) atomic snapshots at 50 ps, and (c) melting speed in different stages.

Figure 12 .
Figure 12.The longitudinal evolution of the melting region during the first 50 ps under both SP and DPs with varying pulse intervals.(a) atomic snapshots at 10 ps, (b) atomic snapshots at 50 ps, and (c) melting speed in different stages.

Figure 13 .
Figure 13.Evolution of electron temperature T e and lattice temperature T l under different pulse settings (a,d) SP, (b,e) DP with a 2 ps interval, (c,f) DP with a 3 ps interval.

22 Figure 14 .
Figure 14.The impact of sub-pulses energy ratio of DP on the melting process.Sub-figures for DPs with (a-c) 1 ps, (d-f) 2 ps, and (g-i) 3 ps intervals, respectively.The first and second columns are atomic snapshots of 0 − 10 ps and 0 − 50 ps.The third column represents the melting speed vs. energy ratio at different stages.

Figure 14 .
Figure 14.The impact of sub-pulses energy ratio of DP on the melting process.Sub-figures for DPs with (a-c) 1 ps, (d-f) 2 ps, and (g-i) 3 ps intervals, respectively.The first and second columns are atomic snapshots of 0 − 10 ps and 0 − 50 ps.The third column represents the melting speed vs. energy ratio at different stages.

Table 1 .
The lattice temperature obtained by the simulation results with different thicknesses and pulse numbers.

Table 1 .
The lattice temperature obtained by the simulation results with different thicknesses and pulse numbers.
C e C e = 7.5209 × 10 2 + 8.8184 × 10 1 T e + 5.5628 × 10 −4 T 2 e − 5.7927 × 10 −8 T 3 e + 9.7671 × 10 −13 T 4 G e−ph Electron-phonon coupling factor, W/ m 3 K i Pulse number of the current calculation, none K e Total kinetic energy in the grid, J k B Pulse interval between the ith and (i − 1)th sub-pulses, s t p Laser-pulse duration, s v i Velocity of atom i, m/s