Ablation study of laser micromachining process with molecular dynamics simulation

Abstract

A three-dimensional molecular dynamics model is presented for the simulation of the creation of a micro-hole on a thin film metal substrate via laser ablation. For the presented analysis, molybdenum and aluminium specimens are selected and short pulses are assumed. The laser fluence takes several values between 0.5 and 20 J/cm². The proposed models include significant laser ablation phenomena such as plasma shielding. However, they are not computationally intense. In this study, the Morse potential is used for the interactions of the atoms of the specimens. The analysis is carried out in order to investigate the ablation rate, the ablation depth and the mean temperature of molybdenum and aluminium targets under their heating by the laser beam, for several different values of fluence. Results for molybdenum indicate that as fluence increases, it takes less time for the atoms to be ablated. For low-fluence pulses, more than one pulse may be required for the ablation of all atoms. For high-fluence pulses, the ablation is not uniform across the entire duration of the pulse and the specimen is overheated. A fluence value around 2–3 J/cm² is suggested for uniform ablation. From the analysis, it is evident that the evolution of ablation and system temperature is different for molybdenum and aluminium, for the same laser fluence. This is attributed to different crystalline structures and absorptivity of each material. It may be said that molecular dynamics prove to be a powerful tool for the simulation of nanomanufacturing processes and useful conclusions are drawn from the analysis.

Keywords

Material removal processes laser ablation simulation molecular dynamics

Introduction

The demand for components that possess features in the micro- or nanometre regime has increased steadily over the past years. This is mainly due to the fact that such components find use in a wide range of applications. The IT-related components are the leading and perhaps most important sector that incorporates these parts. Other sectors such as health and biomedicine, automotive industry and telecommunications follow closely. The applications are microfluidics, pumps and valves, micronozzles, optical components, micromolds and microholes on various materials just to name some.¹ For the production of microcomponents, advanced manufacturing techniques are required. Silicon-based products are mainly produced with lithographic processes; most of these processes are planar or 2.5 D and have limited application to materials other than silicon. Technologies for processing other materials from a few micrometres to a few hundreds of micrometres, namely, metals, ceramics and polymers, are in use. These micromachining processes are abrasive, conventional or non-conventional machining processes. Unwanted material is carried away from the workpiece usually in the form of a chip, while evaporation or ablation may take place in some machining operations.

Laser technology is often cited in micromachining-related scientific works for the production of small-scale parts and tools.² Multiple operations have been successfully realized through laser micromachining and short-pulsed lasers in particular, even by using simple equipment and applying plain manufacturing methods. It is thus proved that the constantly evolving technology may lead to nanometre precision in micromachining processes. Laser micromachining is involved in a broad range of industrial and everyday applications. These include advanced materials processing, large- and small-scale cutting and drilling for biomedical devices and microfluidic applications, micro- and nano-surface configuration, nanotechnology applications, pulsed laser deposition (PLD) of thin films and coatings, micro-joining, surgical operations and artwork renovation.^3–6

In order to remove an atom from a solid through a laser pulse, energy exceeding the binding energy of that atom needs to be applied. Ablation process may depend on the type of material and on laser intensity, wavelength, duration and number of pulses. To explain the laser ablation process, several different mechanisms have been proposed. These include phase explosion,^7,8 critical point phase separation,⁹ spallation¹⁰ and fragmentation.^11,12 Recently, molecular simulation models have been significantly developed to cater for researchers’ needs in enriching their knowledge on laser ablation details and particularities. These include molecular mechanics (MM), Monte-Carlo simulation (MC) and molecular dynamics (MD). The latter of these three methods is characterized as deterministic, as opposed to the other two stochastic methods. The broad use of MD has been encouraged due to its ability to simulate the actual behaviour mechanisms of materials, as a result of interactions in the fundamental molecular level. The MD method exhibits the ability to predict the temporal evolution of a system comprising interacting atoms combined with Newton’s Second Law of Mechanics to reveal information about dynamic attribute changes in molecules. MD has been used for the simulation of various manufacturing processes at nano-scale,^13–15 including laser ablation; an account of the works pertaining to the latter case will be discussed in another section of this article.

In this article, a three-dimensional (3D) MD code is presented for the simulation of laser micromachining, which accounts for the subsidiary phenomena that occur directly from the moment when a single laser pulse contacts the surface of a target. The results pertain to nanosecond laser pulses on spots of thin films made of molybdenum (Mo) and aluminium (Al). Configurations of the code are based on the fact that the laser ablation process is undertaken using short pulses. Furthermore, the attenuation of the phenomenon over time because of the expanding plasma plume and the gradually attenuating radiation thereby caused is considered. The objective behind the realization of this analysis is to investigate, at regular time intervals, the evolution of such parameters as the number of absorbed photons by the material, the rate of atom removal from the main material, that is, the ablation rate, the maximum depth of the crater formed due to the application of the beam, that is, the ablation depth and the mean temperature of the target under its heating by the beam.

Laser ablation and plasma shielding

The aim of every modelling procedure is to deliver a work that is as complete as possible. However, it is rather difficult to incorporate all phenomena taking place during laser ablation within a single computational model. Macroscopic observations give an initial and general view of the evolution of the ablation process. Initially, laser energy is absorbed by the target, with this energy subsequently transforming itself to thermal vibrations. In case the atoms of the target material acquire a considerable amount of energy, it is possible that various particles, for example, electrons, ions, neutral atoms or molecules, are emitted from its surface. Certain circumstances favour the formation of a plasma cloud in front of the target. Laser wavelength and pulse duration, along with the nature of the material, are critical in determining the succession and time scale of the events that occur during the ablation process.

A phenomenon that is usually neglected in the simulation of laser micromachining is plasma shielding, which plays a significant role in the analysis of laser ablation on metal targets. Plasma is formed because of the sublimation of the target material’s atoms, which are heated by the incident laser radiation. During the emission of the first and any subsequent pulses, the plasma is further heated and expanded. This leads to the creation of a protective shield for the material’s upper surface. Due to the plasma plume’s influence, the incident laser energy gradually decreases while the plume keeps expanding. Hence, the heating rate of the target decreases, having a similar impact on the ablation rate, the efficiency and the quality of the ablation process.

Bulgakov and Bulgakova¹⁶ and Liu and Zhang¹⁷ have thoroughly studied the plasma shielding phenomenon, emphasizing on the case of nanosecond laser pulses. Intermediate energy densities between 1 and 10 J/cm² not only cause the local vaporization of the material but also form a plasma plume whose temperature is in the 10⁴ K order of magnitude. It is worth noting that the optical depth of the plasma is governed by such parameters as the electron density, the temperature and the size of the plasma cloud. These parameters significantly vary over time, but possible attempts to model the evolution of the plasma’s optical depth demonstrate an extremely high level of complexity if they are to take time dependencies into consideration. The aforementioned model aims to estimate the energy absorption in the heated plasma. It can serve as a basis on which plasma shielding and its consequences can be quantified. According to this model, the density of the absorbed radiation, E_a, is a parameter that can exclusively characterize the absorption increase caused by the heating of the plasma. As the characteristic time of energy exchange between electrons and ions, for temperatures and ion densities of the respective 10⁴ K and 10¹⁹ cm⁻³ orders of magnitude, is much lower than the pulse durations of the nanosecond lasers, it is assumed that the E_a radiation density is absorbed instantaneously and uniformly among all evaporated particles. The absorption of the laser radiation by the target leads to a temperature increase, ΔT, described by equation¹⁶

Δ T = \frac{(γ - 1) E_{a}}{k N} = \frac{(γ - 1) m E_{a}}{k \cdot Δ z \cdot ρ}

(1)

where γ is the effective adiabatic coefficient, k represents the Boltzmann constant, N is the number of evaporated atoms from a uniform surface, m is the mass of a single atom and ρ is its mass density.

A formula used by Bulgakov and Bulgakova¹⁶ introduces the optical depth of the plasma plume, Λ(t), in order to determine the relation between the intensities of the incident laser beam, I, and the beam originally emitted by the laser device, I₀. According to this study, the numerical integration of the absorption, a, across the depth of the target, z, is used to calculate the plume’s optical depth

I (t) = I_{0} (t) \exp [- Λ (t)] = I_{0} (t) \exp [- \int_{0}^{\infty} a (n, T) dz]

(2)

It is possible to approximate the absorption with a linear product of the particle density, n, and a function f(T) whose value increases with temperature. This approximation is valid only when the temperature of the plasma is low and the plasma lies in a state of equilibrium

a = n \cdot f (T)

(3)

However, according to the analysis by Vasantgadkar et al.,¹⁸ for t_MLP being the time the maximum laser power is attained assumed to occur at the start of the pulse and t_FWHM being the full width at half maximum time, optical depth can be expressed as

Λ (t) = {(\frac{t - t_{MLP}}{t_{FWHM}})}^{2}

(4)

It is worth mentioning that Vasantgadkar et al.¹⁸ presented diagrams that show the evolution rate and maximum value of the target’s temperature, affected by the fluence of the laser beam, both by considering and ignoring the plasma shielding effect for various laser fluencies. If plasma shielding phenomenon is taken into account, then the peak absolute temperature of the target generally drops 4–6 times, compared to the case where this effect is not considered.

MD simulation in laser ablation

MD simulation has been used in laser applications, especially when target ablation by short-pulsed¹⁹ and ultrashort-pulsed lasers^20,21 is considered. The calculation of the interatomic forces during an MD simulation is made possible through an inter-particle interaction potential function, $U ({\vec{r}}_{1}, {\vec{r}}_{2}, \dots, {\vec{r}}_{N})$ , which depends on the displacement $\vec{r}$ of each of the system’s N particles. In the majority of situations, potential U must serve as the exclusive input data for the calculation of the desired quantities during the simulation. Especially as far as the forces are concerned, it is only necessary to compute the gradient of the potential.

Although the MD simulation procedure is advantageous in many aspects, it contains two very important limitations. The first one is the maximum number of atoms that a system can accommodate for the successful completion of the simulation process; the required simulation time is greatly affected by this number. The second one deals with the need to develop a mathematical formula that describes the interaction potential with maximum accuracy.

Most classical MD descriptions use either the embedded-atom method (EAM) or the pair potential approximation (PPA). The EAM approximation has been used by Nedialkov and Atanasov²² in investigating the ablation of Fe targets due to the incidence of femtosecond laser pulses, combined with the laser-induced deep drilling process. This semi-empirical expression helps to identify the factors that influence the final geometrical shape of the hole formed during laser ablation. The first one is the deposition of the ablated material at a specific height of the hole’s internal borders, where the hole becomes narrower. The second one refers to the further degeneration of particles that interact with the ones originally removed, referred as secondary ablation. EAM was used by Schäfer et al.²³ in an attempt to simulate the ablation of Cu targets using picosecond laser pulses. Yamashita et al.²⁴ investigated the heat transfer and shock wave formation processes during femtosecond laser ablation on a thick Al target, using 80,000 Al atoms for the MD simulation conducted. The 12-6 Lennard-Jones expression is a representative example of the PPA application. Potential is calculated by using a formula that depends on the material’s stress and strain parameters, that is, σ and ε, respectively, as well as on the distance r between two adjacent atoms. It is an adaptation of the Mie potential relation, for n = 12 and m = 6, expressed as¹⁵

u (r) = \frac{λ_{n}}{r^{n}} - {\frac{λ_{m}}{r^{m}} |}_{m = 6}^{n = 12} = 4 ε [{(\frac{σ}{r})}^{12} - {(\frac{σ}{r})}^{6}]

(5)

The Morse potential function (MPF) has been used in a wide range of applications, including the interaction analysis between laser beams and materials of metallic and organic nature. The Morse potential is expressed through the following formula²⁵

φ (r_{ij}) = D [e^{- 2 a (r_{ij} - r_{0})} - 2 e^{- a (r_{ij} - r_{0})}]

(6)

where r_ij is the distance between two adjacent atoms; r₀ is the equilibrium distance; D is the dissociation energy, which is the minimum energy value required for the breakup of a bond between two atoms; and a is the dissociation constant. The first derivative of equation (6) in terms of the distance r gives the interaction force between particles, F. The above formula has been constructed in such a way that the force is attractive (F < 0) for high distances, repulsive (F > 0) for low ones and nil (F = 0) for a given equilibrium distance r₀. For this distance, the potential force P(r) must acquire its minimum value.

Studies regarding the use of pair- or many-body interatomic potentials, considering micromachining in general and laser ablation in particular, can be found in the relevant literature.^26,27 It is concluded that many-body interatomic potentials give more realistic treatment of material properties. However, computational cost can increase by several times using a complicated potential function, compared to the application of pair-body potentials. A pair-body potential such as MPF is considered suitable for the description of cubic metals.²⁸ The computational simplicity of MPF is an advantage that has been exploited by researchers, for the simulation of laser ablation.^29–31

With concern to researchers’ studies on laser ablation aspects using the MPF, Cheng and Xu²⁹ investigated the disintegration mechanisms of Ni crystals during laser ablation, comparing and contrasting their differences according to the value of the laser energy. Nedialkov et al.³⁰ developed a theoretical MD model and applied it on three different materials, namely, Al, Ni and Fe, in an attempt to determine laser ablation mechanisms with pulse duration. Stavropoulos and Chryssolouris³¹ used Java platform in order to determine how notable laser ablation characteristics evolve over time, after the incidence of femtosecond beams on Fe surfaces. However, it is worth noting that the possible impact of plasma shielding in the intensity and final results is often neglected.

MD simulation methodology

The implementation of the MD method can take place in commercial- or open-source-specialized MD simulators^32,33 or in an in-house software.^27,29–31 In either case, code writing is required in order to incorporate all the features of the analysis. The models presented in this article were realized through MATLAB. It pertains to the simulation of a rather small but not too computationally intense system, which, however, considers the plasma shielding phenomenon. All parameters within the current analysis underwent a process of non-dimensionalization. The reason for this approach is the considerable difference in orders of magnitude between measurement units used both in the macroscopic and the atomic levels. Therefore, it is necessary to avoid possible confusion between measurement units belonging to different systems and to enclose all numerical values in such a range that does not risk leading to potential problems in precision of calculations.

A theoretical MD model capable of reliably simulating the laser ablation process of metals generally needs to fulfil two principal requirements. First of all, the initial position of each atom is necessary to be described; the cut-off radius that determines which atoms’ interactions will be included in the calculations is defined. Second, the energies of the atoms must be calculated and the boundary conditions that accompany the problem have to be determined. Next, all ablation stages up to the attainment of thermodynamic equilibrium have to be simulated. For the analysis presented in this article, two materials are considered, namely, Mo and Al. These are both cubic metals of body-centered cubic (BCC) and face-centered cubic (FCC) structures, respectively. These structures are respected during the creation of the initial configurations of Mo and Al atoms in the specimens. For the potential energy of a pair of particles, the MPF is given by formula (6), provided that their distance is less than a predetermined cut-off radius, r_c. For higher distances, the potential energy, P(r_ij), between the adjacent atoms i and j equals 0. The MPF parameters for many metals can be found in the relevant literature^25,28 and for Mo and Al are presented in Table 1.

Table 1.

The Morse potential parameters for Mo and Al.

Metal	a (Å⁻¹)	r₀ (Å)	D (eV)
Mo	1.5079	2.976	0.8032
Al	1.1646	3.253	0.2703

The kinetic energy of each atom is obviously calculated by multiplying its mass with the squared value of its vector speed and dividing the product by 2. Hence, the total energy is a result of the kinetic and potential energies’ addition. For each pair of i, j atoms, the following equation is in operation

\begin{matrix} E_{tot, i} = D \\ \sum_{\begin{matrix} j = 1 \\ j \neq i \end{matrix}}^{N} [\exp [- 2 a (r_{ij} - r_{0})] - 2 \exp [- a (r_{ij} - r_{0})]] + \frac{1}{2} m_{i} υ_{i}^{2} \end{matrix}

(7)

However, the Morse coefficient D and the mass- and velocity-dependent terms may disappear during non-dimensionalization, thus facilitating the MD modelling process in general. Periodic boundary conditions are applied on the XY plane, across which the target material is supposed to extend. Free boundary conditions are applied on the upper surface of the material in order to allow the removal of individual particles or atom groups from the target. Finally, reflective boundary conditions on the lower surface ensure that the energy, volume and total number of particles within the atom system are preserved throughout the entire simulation.

The equilibration procedure is responsible for the thermodynamic equilibrium of the system; this is performed by defining the initial velocities of the atoms, solving Newton’s equations, checking the convergence of velocities, computing the lattice temperature and assigning new velocities to the atoms. During equilibration, the vector velocities of all atoms are chosen in such a way that the total momentum of the system is equal to 0. Furthermore, the velocity distribution in atomistic simulations is required to follow a certain statistic distribution in the case of thermal equilibrium, namely, Maxwell–Boltzmann.²⁹ At the initial stage, the velocities assigned to all the atoms are distributed according to Gauss distribution.³¹ Moreover, it was tested that the usage of the Gaussian distribution aids in attaining the desired system equilibrium under a standard environmental temperature value, such as 300 K. Testing was realized by introducing the velocity convergence function of Schommers,³⁴ for a system of N particles, with 3D vector velocities ${\vec{υ}}_{i}$ defined by their x, y and z constituents. The convergence function C(t) takes the following form

C (t) = \frac{\frac{1}{N} \sum_{i = 1}^{N} {[υ_{ix}^{2} + υ_{iy}^{2} + υ_{iz}^{2}]}^{2}}{{[\frac{1}{N} \sum_{i = 1}^{N} [υ_{ix}^{2} + υ_{iy}^{2} + υ_{iz}^{2}]]}^{2}}

(8)

Although a convergence value of 0.2 is recommended,³⁴ for the purpose of the current analysis, the criterion $| C (t) - 5 / 3 | ⩽ 0.25$ is satisfied for the desired thermodynamic equilibrium to be achieved. This leads to a slight precision degradation in the simulation, without significantly affecting its quality.

System’s temperature is expressed through the sum of the mean squares of the velocities of N atoms³¹

T = \frac{1}{N} \frac{m}{3 k_{B}} \sum_{i = 1}^{N} (υ_{ix}^{2} + υ_{iy}^{2} + υ_{iz}^{2})

(9)

As a criterion for continuing the analysis, a 3K difference between the system and the desired temperature needs to be satisfied; if the criterion is not met, then new velocities are applied. Actually, from the early steps of the presented simulations, both criteria described above were met, thus providing an appropriate velocity distribution at each time step. For example, the atomic velocity distribution for fluence 10 J/cm² in Mo, at two different time steps is plotted in Figure 1(a) and (b), along with the theoretic velocity distribution, namely, Maxwell–Boltzmann, for the corresponding conditions. It can be clearly seen that the aforementioned equilibration procedure can lead to an atomic velocity distribution that respects closely the Maxwell–Boltzmann distribution. Thus, it can be concluded that the thermodynamic equilibrium is achieved in every time step during the laser ablation simulation process.

Figure 1.

Atoms’ velocity distribution for laser fluence of 10 J/cm² at (a) 0.35 ns and (b) 0.75 ns simulation time.

During the course of the MD simulation, it is necessary to make the atom system always achieve equilibrium after the completion of each simulation time step. A classical Verlet algorithm or a predictor–corrector method would be relatively suitable to serve the desired purposes. However, the intended goal may be best fulfilled using a Leapfrog Verlet algorithm, thanks to its good numerical stability, its simplicity and convenience, as well as its satisfactory computational cost.

Before the application of the Leapfrog Verlet algorithm, it is vital to computationally model the laser radiation process on the target material, which is assumed to begin after the initial equilibrium condition has been fully accomplished. The procedure of approaching and simulating the radiation is mainly carried out in terms of energy. The beam is simulated both in the time and space domains, following different distributions for each. The time and space profile of the beam follows the Gaussian distribution. Beam’s modelling is complemented via the application of the Beer–Lambert law, which is favoured thanks to the fact that the beam modelling involves both its reflectance, due to the material’s optical properties, and laser beam waste index. This index is not going to be considered, because the beam’s radius does not change over time with regard to its initial value. Ignoring the index is also a result of using dimensions of a nanometre order of magnitude during this study.

For a $t_{p} = t_{2} - t_{1}$ laser pulse duration and a $t_{0} = (t_{1} + t_{2}) / 2$ time centre of the pulse, the term that expresses the Gaussian temporal distribution of the beam is expanded as follows²⁹

S_{p, Gauss} = \frac{1}{t_{p} \sqrt{π}} \exp {(- \frac{t - t_{0}}{t_{p}})}^{2}

(10)

In terms of the beam’s spatial distribution across the z axis, two important factors have to be taken into consideration: the plasma shielding effect, expressed by equations (2)–(4), and the Beer–Lambert law widely used to describe the attenuation of the beam’s intensity when penetrating into the target material. According to this law, the intensity of the radiation is dependent on the one that the beam has on the target’s surface (z = 0).

Assuming that the radiated surface has a circular cross section, $A = π {(r (z))}^{2}$ , and taking into account the fact that the energy of the radiation is a sum of all the energies that the beam’s photons entail, $E_{v} = N_{p} E_{p} = N_{p} hv = N_{p} hc / λ$ , it is possible to conclude that the number of photons changes with the depth inside the target material according to the formula below

N_{p} = \frac{P_{L} δ_{t} λ}{hc} \exp (- β z)

(11)

where P_L represents the laser power, δ_t is the time step, λ is the beam’s wavelength, h is the Planck coefficient and c is the speed of light in vacuum, given the fact that the above formula applies to the beam before its penetration inside the material, where the speed of light is lower.

The Gaussian spatial distribution of the beam in the XY plane is attained using a two-dimensional variable, which depends on the photons’ positions on the x and y axes, the mean values μ and standard deviations σ for each of the axes, and the correlation coefficient ρ. The general format of such an equation was presented by Stavropoulos and Chryssolouris,³¹ but may be further simplified according to various assumptions. First, if the geometry of the beam is circular, with radius r_b, it is obvious that the standard deviation is the same, regardless of the axis. Second, if the centre of the laser beam on the XY plane is considered as a reference point, the mean value μ for the Gaussian distribution in discussion is 0. A precision parameter a_c may be used in order to confront potential precision problems that arise from the use of the Gaussian distribution. When the above parameter is multiplied with the standard deviation value, it gives the diameter of the laser beam and leads to the simplification of the f(x, y) function

f (x, y) = \frac{a_{c}^{2}}{8 π r_{b}^{2}} \exp [- \frac{a_{c}^{2}}{8 r_{b}^{2}} (x^{2} + y^{2})]

(12)

The integration of equation (12) across the x and y axes is capable of providing the possibility for a photon to lie in a given (x, y) location. If R is the index of reflectivity, then absorptivity is evidently 1−R. The general equation that eventually provides the absorbed energy of the material is the following

\begin{matrix} S = (1 - R) [P_{L} δ texp (- β z)] \\ [\frac{a_{c}^{2}}{8 π r_{b}^{2}} \exp (- \frac{a_{c}^{2}}{8 r_{b}^{2}} (x^{2} + y^{2}))] \end{matrix}

(13)

Moreover, multiplying the number of absorbed photons N_AP with the energy of each separate photon, the total photon energy absorbed by the particles is provided. This energy is added to the one currently possessed by the particles before the radiation process starts.

In order to determine whether an atom will be removed from the simulation volume, the cohesion energy criterion is applied; if the sum of the atom’s kinetic energy and energy absorbed from the photons is higher than the atom’s cohesion energy, C, the atom will be removed. In the inverse situation, the atom will remain in the target material and a new velocity will be calculated for it. For Mo, cohesion energy equals −6.82 eV, while the value is considerably different in the case of Al, with a reliable estimate being −3.45 eV.

The new velocities of each atom are calculated by taking the new kinetic energy that it has acquired after laser radiation process and solving the standard kinetic energy equation in terms of the velocity, υ. The Leapfrog Verlet algorithm can subsequently be used in order to advance each atom’s characteristics by a single time step, which is chosen to be equal to a small segment of the entire duration of a nanosecond laser pulse. After the process of recalculating the new positions and velocities of the atoms, the following relation may be used to determine the inter-particle forces between two adjacent atoms in each of the three Cartesian axes

\begin{matrix} [\begin{matrix} f_{ix} \\ f_{iy} \\ f_{iz} \end{matrix}] = 2 aD \sum_{\begin{matrix} i = 1 \\ j \neq i \end{matrix}}^{N} [\begin{matrix} x_{j} - x_{i} \\ y_{j} - y_{i} \\ z_{j} - z_{i} \end{matrix}] \\ \frac{[- \exp (- 2 a (r_{ij} - r_{0})) + \exp (a (r_{ij} - r_{0}))]}{\sqrt{{(x_{j} - x_{i})}^{2} + {(y_{j} - y_{i})}^{2} + {(z_{j} - z_{i})}^{2}}} \end{matrix}

(14)

Using the above equation, it is now easy to calculate the accelerations by dividing the forces by the atomic mass. In Figure 2, a flowchart of the procedure is provided.

Figure 2.

Flowchart for the MD simulation of laser ablation.

Results and discussion

The model described in this article was built in order to account for the laser irradiation and ablation of Mo target, with some key parameters, for example, atomic radius, the Morse potential parameters and macroscopic cohesion energy, being amended for Al. The chosen materials are generally characterized by a stable crystalline structure, at least for ordinary temperatures up to 10⁴ K. The cross section of the Mo specimen is assumed fully rectangular, ignoring any possible geometrical imperfections that require complicated simulation calculations. The purpose of the micromachining process described is the creation of a hole in the specimen, caused by the removal of particles whose contained energy is sufficient to break the cohesion bonds that holds them in the crystalline structure. Atoms not removed after each step, move or vibrate at higher velocities without being detached from the simulation volume. Figure 3 shows the geometry of the initial simulation specimen.

Figure 3.

3D view of the initial simulation specimen to be irradiated with laser.

In order to reduce the computational cost without negatively affecting the precision of calculations, a simple cubic simulation volume containing 15 atoms per dimension or 3375 in total was introduced. Taking into account the interatomic distance of 2.7252 Å in a Mo crystal, it can be concluded that the edge of the formed cube is equal to about 38 Å, assuming a totally homogeneous and isotropic specimen. The behaviour of the material during its irradiation by a single 10-ns laser pulse is simulated. The number of steps can be adjusted to be small in order to increase precision with a subsequent increase in computational time. However, the choice of 200 steps was deemed satisfactory after several trial procedures, as it allows for the adequate monitoring of the ablation procedure and attribute changes without unnecessarily increasing the calculation time and computational cost.

In Figure 4, the evolution of the ablation process is shown when 10 J/cm² laser fluence is applied. A total of 10 initial snapshots of the simulation volume’s format, during the execution of the code, are shown. In snapshots 2 and 4, the bold lines indicate the regions from where the removal of atoms is initiated. The time moments that correspond to each separate snapshot are shown in Table 2, which also gives a representative image of how the number of ablated particles and the temperature of the system evolve over time.

Figure 4.

Laser ablation evolution when irradiated by a single 10-ns laser pulse of 10 J/cm² fluence.

Table 2.

Number of ablated particles and system temperatures at different time steps of the laser ablation for 10 J/cm² fluence.

Snapshot no.	Simulation time step	Simulation time (ns)	Total number of ablated particles	System temperature (K)
0 (start)	0	0	0	300
1	7	0.35	24	3271.6
2	8	0.4	76	3676.38
3	9	0.45	145	4106.75
4	10	0.5	230	4551.1
5	12	0.6	517	5400.87
6	14	0.7	774	6305.48
7	16	0.8	1052	7131.34
8	18	0.9	1360	8170.38
9	20	1	1653	8980.71
10	22	1.1	1908	9923.43
11	25	1.25	2261	11475.87
12	28	1.4	2541	12532.31
13	32	1.6	2837	15387.19
14	36	1.8	3032	15112.57
15	40	2	3160	18819.54
16	45	2.25	3266	22846.09
17	50	2.5	3314	11832.62
18	60	3	3349	7373.97
19	70	3.5	3355	4704.24
20	100	5	3356	6340.41
21 (finish)	200	10	3369	5473.47

From both Table 2 and Figure 4, it can be seen that a rise in the system temperature above melting point of Mo coincides with the outset of particle removal from the simulation volume. A nearly linear rise in temperature and number of ablated particles occurs up to 0.9 ns from the simulation start. This moment is the turning point in the ablation rate, which starts to decrease. During the 2.25 ns moment, the maximum system temperature of about 23,000 K is attained; at later stages, the remaining atoms of the simulation volume start to cool. It is also obvious that over 90% of the volume’s particles have been ablated after 2 ns of laser irradiation, or a 20% of the total pulse duration. Thus, for the given fluence of 10 J/cm², the high absorptivity of the material, combined with the low thickness of the virtual target, contributes to the full and rapid ablation across the depth dimension, with the micromachining process being practically ineffective for the remaining 7–8 ns of the pulse duration.

The MD model for Mo was run several times for various fluence values. The values chosen for the different simulations range between 0.5 and 20 J/cm², in order to determine the different behaviours of the material when ablation takes place under contrasting circumstances. According to the simulation results, it was found that fluences lower than a threshold fluence of 0.4638 J/cm² are practically unable to induce the desired target ablation with a single laser pulse. The threshold fluence value was estimated according to the observation that the maximum temperature of the particle system, for a 0.5 J/cm² fluence, is equal to 3450 K, relatively higher than the melting point of Mo. It was also assumed that a mean target temperature in the region of 3200 K would be sufficient to provoke the fusion of atoms from anywhere on the target surface.

The number of ablated particles and the mean temperature of the specimen for different fluence values are depicted in Figures 5(a) and (b), respectively. It can be observed that both the number of ablated particles and the mean temperature of the specimen were subjected to changes. As the fluence increases, it takes less time for the majority of atoms to be ablated. Furthermore, the desired depth of material removal is reached sooner. A fluence of at least 5 J/cm² is needed for atoms to be certainly removed before the first nanosecond of the process elapses; the outset of this process is significantly delayed for lower fluences. Low values, in addition to this observed delay, do not favour the ablation of all atoms in the target, thus requiring for more than one pulse to act on the specimen in order to achieve their total ablation across the z axis. However, where the ablation rates are initially high for fluences over 5 J/cm², the drop in this rate quickly becomes radical, making the process inefficient after 1.5–2 ns. It is therefore suggested that a fluence value around 2–3 J/cm² best ensures a nearly uniform ablation rate combined with the removal of almost all atoms in the simulation volume, by the end of the single pulse.

Figure 5.

Variation of the (a) ablated particles and (b) mean temperature of the Mo simulation volume’s atoms over time, for different laser fluence values.

The highest values of mean temperature are achieved when the fluence is maximized, as anticipated. In fact, a 30,000K peak temperature is encountered when the fluence is set to 20 J/cm². Nevertheless, this observation does not suggest that this fluence value produces optimum results. The objective is not to overheat the specimen, but to determine the intensity which allows uniform ablation of the target across the entire duration of the pulse. This can also be suggested by looking at Figure 5(a), where the 20 J/cm² fluence leads to the almost total ablation of the simulation volume before 25% of the pulse duration has elapsed. It is also observed that the rise in temperature is almost linear regardless of the fluence, before atoms start to emerge from the simulation volume due to ablation taking place. For example, a rising trend in temperature is apparent throughout nearly the entire simulation length, when the fluence equals 0.5, 0.8 or 1 J/cm²; this can be explained by the fact that it takes a long time before atoms start being ablated from the surface. In addition, the rise in temperature per time step is directly proportional to the fluence, prior to the start of particle removal. After its linear rise, the temperature momentarily drops and then rises again before a new set of atoms is removed. This is most probably the reason for the peak temperature value appearing earlier during the simulation as the fluence rises. Very high peak temperatures naturally lead to an increasing number of atoms being ejected from the surface due to their vaporization, resulting in a rapid drop in the mean temperature of the system after a considerably high peak value has been reached.

Figures 6(a) and (b) demonstrates the differences in the evolution of ablation and the mean system temperature, respectively, between Mo and Al, for the same laser fluence value of 10 J/cm². These differences are justified by the different crystalline structures and characteristics of the two materials, as well as the much higher absorptivity of Mo, namely, 57.14%, compared to Al, which is only 8%. As the absolute energy cohesion value and the melting and boiling points of Al are lower than these of Mo, it is expected that the removal of atoms from the volume will be initiated earlier. Moreover, the temperatures attained on the Al specimen are expectedly low compared to Mo, because the effective fluence of the laser beam, which is absorbed by the target, is approximately seven times lower, taking into account the aforementioned differences in absorptivity.

Figure 6.

MD simulation results on the evolution of (a) the number of ablated particles and (b) the mean specimen temperature, during the irradiation of Mo and Al thin films with 10 ns laser pulse of 10 J/cm² fluence.

Finally, it is worth noting that for a 10 J/cm² fluence, the simulation results for Al suggest that a single pulse is able to ablate all atoms across the target’s depth before 10 ns have elapsed, although the process is not as rapid as in Mo.

Conclusion

Laser micromachining was investigated through the construction of a 3D MD simulation. A model was realized in MATLAB for laser ablation of Mo and Al specimens, incorporating the plasma shielding phenomenon in the analysis. The specimens were modelled as perfect crystals, as homogeneous and isotropic materials of BCC and FCC structures for Mo and Al, respectively. The interatomic distance and potential are altered suitably for each material. Furthermore, several laser fluence values were tested in the models, namely, 0.5, 0.8, 1, 5, 8, 10 and 20 J/cm², for comparing the simulation results. From the simulation procedure, several conclusions are drawn:

In the simulations, the removal of an atom from the system coincides with the rise in the system temperature above the melting point of the material. There is a linear rise in temperature and the number of the ablated atoms, in the first steps of the analysis; then the ablation rate starts to decrease.

The ablation mechanism and depth, as well as the temperature evolution within the system, depend mainly on the laser fluence of the nanosecond pulse acting on the material.

For low values of the laser fluence, below 0.5 J/cm², no target ablation occurs for Mo. For fluence of 0.5–1 J/cm², the ablation process is rather slow. For fluence of 5–20 J/cm², almost the 90% of the atoms are ablated in the first 20% of the pulse duration.

Number of ablated particles and mean temperature of the simulation volume’s atoms increase with increasing fluence. However, large fluence value may lead to specimen overheating.

For fluences over 5 J/cm², ablation rates are initially high; in this case, the drop in ablation rate is fast after 1.5–2 ns, making the process inefficient.

Therefore, a fluence value around 2–3 J/cm² is suggested for ensuring a nearly uniform ablation.

The number of ablated particles and mean specimen temperature are higher for Mo than Al specimen, for the same laser fluence. The significant difference between the reflectivity coefficients of Mo and Al, crystal structure and properties lead to different results on the evolution and outcome of the ablation process for the two materials.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship and/or publication of this article.

References

Mamalis

Markopoulos

Manolakos

. Micro and nanoprocessing techniques and applications. Nanotechnol Percept 2005; 1(1): 31–52.

Jain

Sidpara

Balasubramaniam

. Micromanufacturing: a review – part I. Proc IMechE, Part B: J Engineering Manufacture 2014; 228(9): 973–994.

Balandin

Niedrig

Bostanjoglo

. Simulation of transformations of thin metal films heated by nanosecond laser pulses. J Appl Phys 1995; 77(1): 135–142.

Jain

Dixit

Paul

. Micromanufacturing: a review – part II. Proc IMechE, Part B: J Engineering Manufacture 2014; 228(9): 995–1014.

Zheng

Wang

Jiang

. Research status and application prospects of manufacturing technology for micro–nano surface structures with low reflectivity. Proc IMechE, Part B: J Engineering Manufacture. Epub ahead of print 25July2014. DOI: 10.1177/0954405414542113.

Prakash

Kumar

. Fabrication of microchannels: a review. Proc IMechE, Part B: J Engineering Manufacture. Epub ahead of print 17June2014. DOI: 10.1177/0954405414535581.

Kelly

Miotello

. Comments on explosive mechanisms of laser sputtering. Appl Surf Sci 1996; 96–98: 205–215.

Perez

Lewis

. Ablation of solids under femtosecond laser pulses. Phys Rev Lett 2002; 89: 255504.

Sokolowski-Tinten

Bialkowski

Cavalleri

. Transient states of matter during short pulse laser ablation. Phys Rev Lett 1998; 81: 224–227.

10.

Zhigilei

. Dynamics of the plume formation and parameters of the ejected clusters in short-pulse laser ablation. Appl Phys A 2003; 76: 339–350.

11.

Lorazo

Lewis

Meunier

. Short-pulse laser ablation of solids: from phase explosion to fragmentation. Phys Rev Lett 2003; 91: 225502.

12.

Perez

Lewis

. Molecular-dynamics study of ablation of solids under femtosecond laser pulses. Phys Rev B 2003; 67: 184102.

13.

Komanduri

Raff

. A review on the molecular dynamics simulation of machining at the atomic scale. Proc IMechE, Part B: J Engineering Manufacture 2001; 215: 1639–1672.

14.

Sun

Chen

Cheng

. Multiscale simulation on nanometric cutting of single crystal copper. Proc IMechE, Part B: J Engineering Manufacture 2006; 220(7): 1217–1222.

15.

Markopoulos

Kalteremidou

KAL

. Molecular dynamics modelling of nanometric cutting. Key Eng Mat 2014; 581: 298–303.

16.

Bulgakov

Bulgakova

. Thermal model of pulsed laser ablation under the conditions of formation and heating of a radiation-absorbing plasma. Quantum Electron 1999; 29(5): 433–437.

17.

Liu

Zhang

D-M

. Vaporization and plasma shielding during high power nanosecond laser ablation of silicon and nickel. Chinese Phys Lett 2008; 25(4): 1368–1371.

18.

Vasantgadkar

Bhandarkar

Joshi

. A finite element model to predict the ablation depth in pulsed laser ablation. Thin Solid Films 2010; 519(4): 1421–1430.

19.

Gamaly

Madsen

Duering

. Ablation of metals with picosecond laser pulses: evidence of long-lived non-equilibrium conditions at the surface. Phys Rev B 2005; 71: 174405.

20.

Ohmura

Fukumoto

Miyamoto

. Molecular dynamics simulation of ablation process with ultrashort-pulse laser. In: Proceedings of SPIE 4088: first international symposium on laser precision microfabrication, Saitama, Japan, 14 June 2000, p.84. Bellingham, WA: SPIE.

21.

Bulgakova

Burakov

Meshcheryakov

. Theoretical models and qualitative interpretations of Fs laser material processing. J Laser Micro/Nanoeng 2007; 2(1): 76–86.

22.

Nedialkov

Atanasov

. Molecular dynamics simulation study of deep hole drilling in iron by ultrashort laser pulses. Appl Surf Sci 2006; 252(13): 4411–4415.

23.

Schäfer

Urbassek

Zhigilei

. Metal ablation by picosecond laser pulses: a hybrid simulation. Phys Rev B 2002; 66: 115404.

24.

Yamashita

Yokomine

Ebara

. Heat transport analysis for femtosecond laser ablation with molecular dynamics-two temperature model method. Fusion Eng Des 2006; 81(8–14): 1695–1700.

25.

Chiang

K-N

Chou

C-Y

C-J

. Prediction of the bulk elastic constant of metals using atomic-level single-lattice analytical method. Appl Phys Lett 2006; 88: 171904.

26.

Oluwajobi

Chen

. The effect of interatomic potentials on the molecular dynamics simulation of nanometric machining. Int J Autom Comput 2011; 8(3): 326–332.

27.

Imamova

Atanasov

Nedialkov

. Molecular dynamics simulation using pair and many body interatomic potentials: ultrashort laser ablation of Fe. Nucl Instrum Meth Phys Res B 2005; 227: 490–498.

28.

Girifalco

Weizer

. Application of the Morse potential function to cubic metals. Phys Rev 1959; 114(3): 687–690.

29.

Cheng

. Mechanisms of decomposition of metal during femtosecond laser ablation. Phys Rev B 2005; 72: 165415.

30.

Nedialkov

Imamova

Atanasov

. Mechanism of ultrashort laser ablation of metals: molecular dynamics simulation. Appl Surf Sci 2005; 247(1–4): 243–248.

31.

Stavropoulos

Chryssolouris

. Molecular dynamics simulations of laser ablation: the Morse potential function approach. Int J Nanomanuf 2007; 1(6): 736–750.

32.

Eder

Bianchi

Cihak-Bayr

. An analysis method for atomistic abrasion simulations featuring rough surfaces and multiple abrasive particles. Comput Phys Commun 2014; 185(10): 2456–2466.

33.

Norman

Starikov

Stegailov

. Atomistic simulation of laser ablation of gold: effect of pressure relaxation. J Exp Theor Phys 2012; 114(5): 792–800.

34.

Schommers

. Molecular dynamics and the study of anharmonic surface effects. In: Schommers

von Blanckenhagen

(eds) Structures and dynamics of surfaces I: topics in current physics, vol. 41. Berlin, Heidelberg: Springer, 1986, pp.199–243.