Innovatives Supercomputing in Deutschland
inSiDE • Vol. 10 No. 1 • Spring 2012
current edition
about inSiDE
index  index prev  prev next  next

Laser Ablation of Metals: Generation of Drops and Bubbles

With the invention of the CO2 gas laser in 1964 lasers became an important tool for industry. Now a beam with a wave length around 10 µm and several kW was available which could be focused down to a fraction of a millimeter. Especially for the processing of materials it became a valuable tool for example for cutting, drilling, welding, structuring and engraving.

With the development of ultra-short pulsed laser sources new possibilities arose: Short pulses are favorable since the interaction time of the laser with the material is reduced drastically, but the heat conductivity is unchanged. The result is a locally limited overheating with efficient materials removal but without warming of the whole workpiece. For example in eye-surgery the treatment with conventional lasers leads to huge pains, while the therapy with ultra-short laser pulses is completely pain free.

Figure 1: Snapshots of a 60 million atom simulation 4, 40, and 70 ps after a 100 fs laser pulse hits the surface. The colors indicate the kinetic energy of the atoms and are a measure for the temperature.

Each reduction of the the pulse duration increased the difficulties to predict the result of the materials processing theoretically. Pulses on the scale of nano- to micro-seconds lead to conventional melting and subsequent expelling of the melt. These classical phenomena can be modelled accurately by continuum mechanics methods, and solved with finite-element methods. They rely on empirical equations of states which are often not well known for extreme temperatures and pressures. A drawback is that continuum methods cannot give information on the atomistic scale, for example about the formation of gas or liquid drops. But such information is crucial since the drops may fall back into the hole and might clog it. At this place molecular dynamics simulations (MD) are the method of choice: By applying classical equations of motion and empirical potential interactions and integrating the equations at discrete time-steps samples up to billions of atoms can be handled. However, pulses on the scale of pico- down to femto-seconds lead to new challenges, for example for metals, because their properties can no longer be treated by classical physics. Thermal properties like heat conductivity and heat capacity are dominated by the free electron gas and cannot be modelled directly by classical molecular dynamics simulations. Furthermore, the laser beam is absorbed directly by the free electron gas. First the electrons are excited to no-equilibrium, second they equilibrate and then the energy is transmitted to the ion lattice which heats up and melts eventually.

A continuum model solution to this problem was presented by Anisimov et al. [1] in 1974: Take the standard heat conduction equation C dT/dt=grad (K grad T), (C: heat capacity, K: heat conductivity, T: temperature) split it up into separate equations for electron Te and lattice temperatures Tl and add a coupling term g (Te-Tl) to both equations. The method also solves the problem how to introduce the laser energy: Simply add it as a source term to the electron equation. The two-temperature model (TTM), as it is called, works if the electron relaxation is fast enough to reach an equilibrium which allows to define an electron temperature. The TTM model permits an adequate description of laser ablation with a very fast rise of the electron temperature and a subsequent transferral of the energy to the lattice. But it cannot tell us anything about the atomistic processes like gas and drop emission, or the generation of defects in the bulk around the hole.

The obvious solution is a combination of MD and TTM into a multi-scale method: Keep the equation for the electrons and solve this continuum problem with a finite difference method, and replace the equation for the lattice by the equations of motion of MD. The coupling term is introduced into MD as a friction or anti-friction term depending on the sign of the coupling. Only the coupling constant g has to be translated to atomistic observables.

Our simulation code IMD [2] has been upgraded by the TTM+MD model. To be able to define a temperature an average of the kinetic energy of about 100 atoms has to be taken. Stability criteria require that the electron equations have to be solved several hundred times during one MD time-step, but the time fraction required for the TTM part is still negligible.

Figure 2: Snapshot of an ablation simulation with inhomgeneous heating 60 ps after the laser pulse hit the sample. The size of the sample is 48 million atoms in a box of 80x80x180 nm³. The gas phase of the plume extends up to 1.5 µm above the surface. Colors indicate an electron density function.

The TTM+MD model has been used only to simulate samples with homogeneous illumination by the laser beam [3]. The reason is the following: the fast electron conductivity would remove a transverse temperature gradient immediately. Currently the only back-door to simulate the much more interesting inhomogeneous illuminated samples is to switch off the electron heat conduction which essentially means returning to classical MD without TTM and adding the energy supplied by the laser to the kinetic energy of the electrons. This is the so-called rescale model.

Figure 3: Snapshot of a simulation 80 ps after the laser beam hit the sample. The colors indicate the kinetic energy of the atoms and are a measure for the temperature.

Before we present results of simulations and quantitative analysis of the ablation plume and the formation of bubbles we shortly address the visualization and analysis of the results. Many programs exist to visualize macromolecules and atomistic solids, but none of them can handle multi-million atom samples on a standard workstation. This has been made possible by Megamol developed by the Visualization Research Center of the University Stuttgart [4]. A MD simulation produces only coordinates and velocities of atoms. To figure out which atoms belong to the same drop an algorithm has to be found which identifies the cluster. We successfully applied the DBSCAN (Density Based Spatial Clustering of Applications with Noise) algorithm. It depends on a minimal distance and a minimal number of points which are required to form a cluster. If a certain atom has enough neighbors then the points are density-connected. Points belong to the same cluster if they belong to chains of density-connected points. But how to identify empty bubbles? We superimposed a regular grid with mesh size of the order of the atom distance. If a grid point is close to an atom, it is removed. The remaining grid points represent a negative copy of the bubbles.

Simulation of Gas and Drop Production

Samples with cross section of 101x101 nm² and 108 nm length, containing 60 million atoms have been simulated [5]. Five different laser intensities have been studied and five simulations have been carried out at each intensity to improve statistics (Figure 1-4). A single snapshot takes about six GB of disk space, thus about 6 TB were required to make a movie of one of the simulations. The ablation plume contain from 250,000 up to 850,000 atoms for the lowest to the highest laser intensity. We clearly find a bimodal distribution (Figure 5) with a gas fraction and a cluster fraction, both obeying a power law with exponents that agree with experimental results. Another important result is the angular distribution of the material: with increasing laser intensity the plume concentrates more and more to the direction of the incoming laser beam. Again the results can be fitted very well with the prediction of phenomenological laws derived from experiment.

Figure 4: Snapshot of an ablation simulation with inhomgeneous heating 60 ps after the laser pulse hit the sample. The size of the sample is 48 million atoms in a box of 80x80x180 nm³. The gas phase of the plume extends up to 1.5 µm above the surface. Colors indicate the spatial depth of the sample. Figure 5: Cluster size distribution for five different laser intensities. For the ablation yield Y(N) a power law can be found with different exponents for the low-mass and high- mass (N>10) ranges.

Simulation of Bubble Formation

Here the samples contained 7 million atoms in an area of 27x27 nm² and 217 nm length. If the laser intensity is below a certain ablation threshold then it only heats the material. This leads to the formation of fluid bulges, a process which is also used to pattern sur- faces. The challenge is to form stable bulges by recrystallization, a result which has not yet been achieved in simulations. The results of our simulations resemble the classical Ostwald ripening process (Figure 6, 7): First many small bubbles form everywhere. Then the size of the bubbles grows by coalescence and their number shrinks. Finally the bubbles disappear since the ablated layer bounces back to the bulk. The interesting observation is that there are two populations of bubbles: the smaller sized grows and vanishes after about 25 ps for all intensities while the larger sized has a maximum after about 30 ps and stays longer for higher laser intensities (Figure 8).

Figure 6: Evolution of the voids at 11, 13, 25 and 44 ps (from left to right).


A simulation of a 60 million atom sample for 200 ps (equal to 200,000 iteration steps) took more than 60 h on 512 nehalem cores. Thus running 25 simulations consumed of the order of one million CPU hours, including special runs and extra runs for movies. Our goal is to bridge the gap between simulations and experiment together with the Institut für Strahlwerkzeuge (IFSW). While the experimentalists have to try to get the laser spots smaller, we want to increase the sample size to a few µm. It is not necessary to increase the depth of the sample since the length is determined by the heat conductivity Thus only the cross section has to be grown by a factor of a hundred at least, leading to sample sizes of 10 billion particles with about a TB of storage required for a single snapshot. Such a simulation is well within the domain of the new Cray XE6 of the HLRS, however, it causes new challenges to visualization, animation and analysis of the results even if only the interesting part of the sample is extracted.

Figure 7: Series of bubble nucleation, coalescence and disappearance for 400 fs laser pulses. Colors indicate the temperature of the sample. The greenish part at the left is molten.

Figure 8: Evolution of the total bubble volume in time. For six pulses of the same intensity but different pulse duration.


We acknowledge collaboration with the S. Grottel (Megamol), and his colleagues from VISUS. The project has been supported by the DFG through project B5 of SFB 716. We are also thankful to the HLRS for providing us with computation time.


[1] Anisimov, S. I., Kapeliovich, B. L., Perel'man, T. L. Electron emission from metal surfaces exposed to ultra short laser pulses, Soviet Physics - JETP 39 (1974) 375-7

[2] Roth, J., Gähler, F., Trebin, H.-R., Trebin, A. Molecular dynamics run with 5,180,116,000 particles, Int. J. Mod. Phys., C 11 (2000), 317-22,

[3] Sonntag, S., Roth., J., Gähler, F., Trebin, H.-R. Femtosecond laser ablaation of aluminum, Applied Surface Science 255 (2009), 9742-4

[4] Grottel, S., Reina, R., Vrabec, J., Ertl, T. Visual verification and analysis of cluster detection for molecular dynamics, IEEE Trans. On Visual. and Comp. Graph, 13 (2000), 1624-31

[5] Sonntag, S., Trichet Paredes, C., Roth, J., Trebin, H.-R. Molecular dynamics simulations of cluster distribution from femtosecond laser ablation in aluminum, Applied Physics A 104 (2011), 559-65

[6] Karlin, J., Roth, J., Trebin, H.-R. Simulation of laser ablation in aluminum: Void production below threshold. Submitted to Proc. of the 11th Conf. on Laser Ablation (Cola11) in Playa del Maya, Mexico

• Johannes Roth
Institut für Theoretische und Angewandte Physik, Universität Stuttgart

• Steffen Sonntag
Robert Bosch GmbH

• Johannes Karlin
Institut für Theoretische und Angewandte Physik, Universität Stuttgart

• Hans-Rainer Trebin
Institut für Theoretische und Angewandte Physik, Universität Stuttgart

top  top