Experimentally constrained multidimensional simulation of laser-generated plasmas and its application to UV nanosecond ablation of Se and Te

The spatiotemporal evolution of laser-generated plasmas is a complex phenomenon that couples thermodynamic, electromagnetic, and quantum mechanical processes. Progress in tackling this problem has been made over the past two decades by the implementation of increasingly sophisticated models, enabled by advances in computational techniques and resources [25–33, 49, 50]. In broad strokes, simulation of the plasma formed upon laser irradiation of a solid requires the integration of processes of ejection of ablated species from the material, with processes of plasma plume ignition and expansion. These events occur with significant mutual interference. Processes of material removal involve light–matter interactions that depend on the laser wavelength, pulse duration, and target material properties. Thermal evaporation, subsurface boiling, supercritical boiling, as well as electronic and hydrodynamic sputtering vary in their relative importance in multiple regimes [51]. It is generally assumed that for sufficiently low laser fluence, thermal evaporation models are adequate for describing ablation by nanosecond pulses [52], whereas additional mechanisms need to be considered in the high fluence regime [53, 54]. If temperatures approach or exceed the critical temperature of the target, for example, supercritical processes take place. This requires models that include both, surface and volume removal mechanisms [28]. In all cases, the physical properties of the target, such as reflectivity, absorption coefficient, and density, develop temperature and time dependencies, which need to be taken into consideration [55, 56]. Once the material has been ejected, models of plasma formation and expansion are needed. Numerous collisional and radiative processes, including single- and multi-photon absorption, impact ionization, electron impact excitation/deexcitation, inverse bremsstrahlung, and plasma reflection may affect the electron and ion populations [57]. Laser plasmas may also contain molecular clusters and nanoparticles, calling for inclusion of molecular absorption and scattering processes [58]. The evolution of the plasma then needs to be calculated using kinetic or fluid models. For fluid models, local thermodynamic equilibrium (LTE) [52] or two-temperature assumptions are generally required for practical simulations [57, 59]. Beyond fluid descriptions, electromagnetic effects are also known to be significant. The emission of fast electrons prior to lattice melting, for example, leaves the target with a positive charge [60]. Space charge caused by this emission, as well as ion and prompt electron ejection from early plasma formation, has been implicated in the acceleration and broadening of the kinetic energy distribution of ions [61, 62]. Analysis of the plasma expansion may therefore need space charge corrections. Finally, geometrical effects can also modify the plasma evolution. Changes in laser spot size alter the interaction volume between the laser pulse and the initial vapor, and may also affect the efficiency of evaporation due to heat dissipation changes in the target [31, 63]. Competing effects lead to changes in the forward peaking of the expansion and can result in important variations of the plasma plume angular distribution [64, 65].

This broad variety of processes presents a challenge to plasma modellers. The inclusion of more intricate mechanisms in 3D simulations have high computational cost and generally involves the introduction of parameters whose values are experimentally unknown and whose ab initio computation is impractical or involve large uncertainties. While the gradual incorporation of more refined processes will continue to yield useful insights as computational power grows, in this paper we employ an effective model that allows fast predictive simulations to guide materials processing. This approach acknowledges the difficulty of accurately capturing all relevant processes involved in laser ablation of arbitrary materials. Instead, we adopt a model with basic well-known process constituents, and incorporate an additional effective mechanism in the laser-plasma interaction—which is a convenient entry point. The high computational efficiency permits a broad sweep on the effective parameter introduced in the model. Its appropriate value for specific irradiation conditions and target material is determined by constraining basic simulation outcomes to Langmuir probe experimental data. This approach is compatible with introduction of other effective parameters related to mass removal process from the target or further aspects of laser-plasma interactions. This method for studying laser plasmas, heretofore inhibited by the need for numerous recurrent and time-consuming simulations, is enabled by our efficient ACM implementation framework.

The effective model couples laser-induced surface evaporation with fluid plasma expansion based on previous formulations [23–25]. The pertinent equations are expressed in 2D-axisymmetric coordinates.

The target material is irradiated at normal incidence with a laser pulse of temporal profile

$$\begin{equation}{I}_{\text{laser}}(t)={I}_{0}\enspace \mathrm{exp}\left[-{\left(\frac{{(t-{t}_{\mathrm{c}})}^{2}}{2{\sigma }_{t}^{2}}\right)}^{P}\right],\end{equation} \tag{ 1 }$$

with t_c = 15.0 ns, σ_t = 9.05 ns, and P = 7, corresponding to a pulse width of ∼25 ns. The set value of the parameter P yields a 'top-hat' temporal profile that approximates typical experimental laser pulses. The value of I₀ is used to set the peak irradiance desired for each simulation. The laser intensity on the target surface $I(\mathfrak{r},\mathfrak{z},t)$, is centered at the origin of the radial coordinate $\mathfrak{r}$, with a radial fall-off for growing $\mathfrak{r}$ prescribed by the function

$$\begin{equation}I(\mathfrak{r},\mathfrak{z},t)={I}_{\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{e}\mathrm{r}}(t)\mathrm{exp}\left[-{\left(\sqrt{\mathfrak{r}/{R}_{0}}\right)}^{A}\right](1-\mathcal{R}){\text{e}}^{-\alpha \mathfrak{z}}\end{equation} \tag{ 2 }$$

with R₀ representing the laser spot radius and a choice of A = 12, which approximately reproduces the radial fall-off of standard high-power lasers. The interaction of the laser pulse with the subsurface region of the target is modeled by the Beer–Lambert factor included in equation (2), where the distance from the target surface is given by $\mathfrak{z}$ and the material's absorption coefficient is denoted as α. A fixed target reflectivity $\mathcal{R}$ is assumed.

Heat transfer from the absorption volume into the bulk of the target is computed by solving the heat diffusion equation

$$\begin{equation}\frac{\partial T(\mathfrak{r},\mathfrak{z},t)}{\partial t}=\nabla \cdot \left[\frac{\kappa }{({C}_{\mathrm{p}}{\rho }_{m})}\nabla T(\mathfrak{r},\mathfrak{z},t)\right]+\frac{\alpha }{{C}_{\mathrm{p}}{\rho }_{m}}I(\mathfrak{r},\mathfrak{z},t),\end{equation} \tag{ 3 }$$

where κ is the thermal conductivity, C_p is the specific heat, and ρ_m is the mass density of the target material.

The pressure of the vapor that forms at the surface of the target p_vap is calculated from the Clausius–Clapeyron equation

$$\begin{equation}{p}_{\text{vap}}({T}_{\mathrm{s}})={p}_{0}\enspace \mathrm{exp}\left[\frac{{\Delta}{H}_{\mathrm{l}\mathrm{v}}({T}_{\mathrm{s}}-{T}_{\mathrm{b}})}{R{T}_{\mathrm{s}}{T}_{\mathrm{b}}}\right],\end{equation} \tag{ 4 }$$

where ΔH_lv, T_b, and T_s stand for the heat of vaporization, boiling point, and surface temperature, respectively, and R is the ideal gas constant. The vapor mass density ρ_vap,s at the surface is calculated using the ideal gas law

$$\begin{equation}{\rho }_{\mathrm{v}\mathrm{a}\mathrm{p},\mathrm{s}}=\frac{m{p}_{\text{vap}}}{k{T}_{\mathrm{s}}},\end{equation} \tag{ 5 }$$

where m is the atomic mass and k is the Boltzmann constant. Assuming a Maxwellian velocity distribution, the escaping atoms have normal velocity components whose mean is given by

$$\begin{equation}{\mathrm{v}}_{\mathrm{v}\mathrm{a}\mathrm{p},\mathrm{s}}=\sqrt{\frac{2k{T}_{\mathrm{s}}}{\pi m}}.\end{equation} \tag{ 6 }$$

Ionization of the vapor is governed by the set of Saha equations

$$\begin{equation}\frac{{X}_{\mathrm{e}}{X}^{z}}{{X}^{z-1}}=\frac{1}{{n}_{\text{vap}}}{\left(\frac{2\pi {m}_{\mathrm{e}}kT}{{h}^{2}}\right)}^{3/2}\mathrm{exp}\left(-\frac{{\mathrm{I}\mathrm{P}}^{z}}{kT}\right),\quad z=1,2\end{equation} \tag{ 7 }$$

where X_e and X^z are the fraction of electrons and ionic species with charge z, respectively. The number density of the vapor n_vap is calculated from the mass density of the vapor ρ_vap and the atomic mass. Ionization potentials for atoms with charge z are represented by IP^z and T is the plasma temperature. Matter and charge conservation are enforced by ${\sum }_{z=0}^{2}{X}^{z}=1$ and X¹ + 2X² = X_e.

The laser irradiance that reaches the target surface is attenuated via a Beer–Lambert relation by absorption in the plasma. Four terms are included as contributions to this absorption: (i) electron-neutral and (ii) electron–ion inverse bremsstrahlung, (iii) single-photon photoionization, and (iv) an effective absorption process proportional to the number density of neutral species in the plasma. The latter is meant as a stand-in term, accounting for more complex processes of light absorption in the plasma, such as multi-photon ionization, and secondary effects such as electron impact excitation, whose cross-sections are generally unknown or difficult to compute [57]. Because plasma absorption strongly impacts material removal from the target, this term may also provide an effective accounting of thermal processes beyond evaporation and non-thermal effects of laser-target interaction.

The total plasma absorption coefficient α_absorb may therefore be written as

$$\begin{equation}{\alpha }_{\text{absorb}}={\alpha }_{\mathrm{I}\mathrm{B},\mathrm{e}-\mathrm{n}}+{\alpha }_{\mathrm{I}\mathrm{B},\mathrm{e}-\mathrm{i}}+{\alpha }_{\mathrm{P}\mathrm{I}}+{\alpha }_{\text{eff}},\end{equation} \tag{ 8 }$$

where α_IB,e-n and α_IB,e–i are the absorption coefficients for electron-neutral and electron–ion inverse bremsstrahlung, respectively, α_PI corresponds to absorption by photoionization, and α_eff is the effective absorption coefficient.

The inverse bremsstrahlung terms are evaluated by [25, 63]

$$\begin{equation}{\alpha }_{\mathrm{I}\mathrm{B},\mathrm{e}-\mathrm{n}}=\left[1-\mathrm{exp}\left(-\frac{hc}{\lambda kT}\right)\right]Q{n}_{\mathrm{e}}{n}^{0}\end{equation} \tag{ 9 }$$

$$\begin{align}{\alpha }_{\mathrm{I}\mathrm{B},\mathrm{e}--\mathrm{i}}& =\left[1-\mathrm{exp}\left(-\frac{hc}{\lambda kT}\right)\right]\frac{4{e}^{6}{\lambda }^{3}{n}_{\mathrm{e}}}{3h{c}^{4}{m}_{\mathrm{e}}}\\ & \quad \times {\left(\frac{2\pi }{3{m}_{\mathrm{e}}kT}\right)}^{1/2}\sum\limits _{z=1}^{2}{z}^{2}{n}^{z},\end{align} \tag{ 10 }$$

where h is Planck's constant, c is the speed of light in vacuum, and e is the elementary charge. The laser wavelength is represented by λ. The density of electrons and species with charge z are n_e and n^z , respectively. Ionization states up to z = 2 are accounted for, which is expected to be sufficient for relatively low laser irradiances. Q is the cross section for photon absorption by an electron during a collision with a neutral atom. It has dependencies on the electron and photon energies, as well as on the momentum transfer cross section of the atom [66]. Within the scope of the effective model, we adopt a fixed value of Q = 10⁻⁴⁰ cm⁵ for all simulations. This represents an estimate of the mean value of Q for the variety of temperatures and chemical species of interest in our study, at the laser wavelength used [66].

The absorption coefficient for photoionization α_PI is computed by [63]

$$\begin{equation}{\alpha }_{\mathrm{P}\mathrm{I}}=\sum\limits _{z=0}^{2}\sum\limits _{i={N}_{\mathrm{min}}^{z}}^{{N}_{\mathrm{max}}^{z}}{x}_{i}^{z}{n}^{z}{\sigma }_{i}^{z}(T),\end{equation} \tag{ 11 }$$

where ${N}_{\mathrm{m}\mathrm{i}\mathrm{n}}^{z}$ is the minimum energy level of an atom with charge z that can be ionized to charge z + 1 by absorbing a single photon of energy 5 eV (248 nm), ${N}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{z}$ is the maximum energy level considered. The fractional occupation of level i is calculated by ${x}_{i}^{z}=({g}_{i}^{z}/{U}^{z})\mathrm{exp}[-{E}_{i}^{z}/(kT)]$ where U^z is the partition function of charge state z and ${g}_{i}^{z}$ is the statistical weight of level i and charge z. The photoionization cross section of level i with charge z, ${\sigma }_{i}^{z}$, is calculated by [63]

$$\begin{equation}{\sigma }_{i}^{z}(T)=\frac{32{\pi }^{2}{(z+1)}^{2}{e}^{6}{k}_{\mathrm{e}}^{3}}{3\sqrt{3}{h}^{4}c{\nu }_{\text{laser}}^{3}{g}_{i}^{z}}{U}^{z+1}(T)\frac{\mathrm{d}{E}_{i}^{z}}{\mathrm{d}i},\end{equation} \tag{ 12 }$$

where k_e is the Coulomb constant, ν_laser is the laser frequency, and $\mathrm{d}{E}_{i}^{z}/\mathrm{d}i$ is the spacing between energy levels.

Finally, the effective absorption coefficient α_eff is assumed to be proportional to the number density of neutral species n⁰, and expressed as

$$\begin{equation}{\alpha }_{\text{eff}}=\begin{cases}0\quad & \quad \mathrm{i}\mathrm{f}\enspace T< {T}_{\mathrm{min}}\\ {C}_{\text{eff}}{n}^{0}\quad & \quad \mathrm{i}\mathrm{f}\enspace T > {T}_{\mathrm{min}}\end{cases},\end{equation} \tag{ 13 }$$

where T_min is a cut-off temperature that excludes regions too cold for the vapor to be characterized as a plasma. In the simulations presented in this paper T_min was set to 1000 K.

The proportionality constant C_eff is used as an adjustable parameter to constrain the model by experiments. The value of C_eff for given irradiation conditions and target material is determined by matching the main features of simulated ion time-of-flight (TOF) kinetic energy distributions to experiments, as detailed in section 4.2.

Expansion of the ablated material is described by the equations of hydrodynamics, representing conservation of mass (equation (14)), momentum (equation (15)), and energy (equation (16)), where p is pressure, v is flow velocity, and ρ, ρ v, and ρE are mass, momentum, and internal energy density, respectively

$$\begin{equation}\frac{\partial \rho }{\partial t}=-\nabla \cdot \left(\rho \mathbf{v}\right)\end{equation} \tag{ 14 }$$

$$\begin{equation}\frac{\partial \rho \mathbf{v}}{\partial t}=-\nabla p-\nabla \cdot (\rho \mathbf{v}\otimes \mathbf{v})\end{equation} \tag{ 15 }$$

$$\begin{align}\frac{\partial }{\partial t}\left[\rho E+\frac{1}{2}\rho {\mathbf{v}}^{2}\right]& =-\nabla \cdot \left[p\mathbf{v}+\mathbf{v}\left(\rho E+\frac{1}{2}\rho {\mathbf{v}}^{2}\right)\right]\\ & \quad +{\alpha }_{\text{absorb}}{I}_{\text{laser}}(t)-{{\epsilon}}_{\text{rad}}.\end{align} \tag{ 16 }$$

In equation (16), energy is added to the plasma by α_absorb I_laser(t), with α_absorb given by equation (8), and is lost by bremsstrahlung radiation _rad. Assuming a Maxwellian velocity distribution for the electrons, _rad is evaluated using [67]

$$\begin{equation}{{\epsilon}}_{\text{rad}}={\left(\frac{2\pi kT}{3{m}_{\mathrm{e}}}\right)}^{1/2}\frac{32\pi {e}^{6}}{3h{m}_{\mathrm{e}}{c}^{3}}{n}_{\mathrm{e}}\sum\limits _{z=1}^{2}{z}^{2}{n}^{z}.\end{equation} \tag{ 17 }$$

The model does not invoke the Knudsen layer, where transition from kinetic (non-LTE) to fluid conditions (LTE) occurs. We assume LTE conditions and the validity of the fluid dynamics description throughout the entire simulation domain. In addition, we neglect the effects of finite viscosity (such as diffusion) in the boundary layer as well as in the flow bulk. Because such layers are extremely thin under our high-density laser plasma conditions, these assumptions are expected to be satisfactory when the dynamics of the system far away from the boundaries is of main interest. This approach has been used in hydrodynamic modeling of laser-induced plumes [50, 52, 68] and is common in computational fluid dynamics problems in general [69]. It can lead to an overestimation of the net evaporation flux by ∼15%–20% (due to neglect of material back diffusion) and absence of proper description of the nm-scale physics. However, these limitations are likely to be compensated by the effective model's ability to accommodate uncertainties in the initial and boundary conditions, and by the simulation's numerical efficiency, allowing useful modeling of large-scale gas dynamics properties of the studied system. This cm-to-m long-distance scale effective model, and efficient simulation, can be coupled to nm-scale physics models [28, 70] for evaluation of the impact of nanometer processes in the long-term plasma behavior.

The pressure p and temperature T of the expanding plasma are related to ρ and ρE by [71]

$$\begin{equation}p=(1+{X}_{\mathrm{e}})\frac{\rho kT}{m}\end{equation} \tag{ 18 }$$

$$\begin{equation}\rho E=\frac{\rho }{m}\left[\frac{3}{2}(1+{X}_{\mathrm{e}})kT+{\mathrm{I}\mathrm{P}}^{1}{X}^{1}+({\mathrm{I}\mathrm{P}}^{1}+{\mathrm{I}\mathrm{P}}^{2}){X}^{2}\right],\end{equation} \tag{ 19 }$$

where X¹ and X² are the fractions of singly and doubly ionized atoms. IP¹ and IP² denote the first and second ionization potentials of the atoms.

The simulation is carried out in a multidimensional setting using a state-of-the-art, open-source ACM framework [72]. The compressible solvers available in this framework (Riemann and all-Mach solvers [72]) were adapted to include the equations of state of the plasma (equations (7), (18) and (19)). The use of the dynamic ACM approach is essential since the initial (t < 30 ns) spatial resolution required near the target is tenths of microns, while the simulation domain is larger by multiple orders of magnitude. As the plasma plume expands (t > 30–60 ns), the fine near-target resolution is no longer required for most of the simulation domain which allows the grid to be coarsened. However, proper resolution of the moving plasma front is still necessary throughout the entire plasma dynamics until the front reaches the anticipated substrate location. Only the very narrow plasma front region needs to be resolved for this purpose, while the remainder of the computational mesh can stay coarse. The numerical efficiency of the ACM framework is further enhanced by an adaptive computation domain capability. This capability works in combination with ACM and allows for the size of the computational domain to automatically increase when the plasma front approaches its boundary. The technique is versatile (e.g. the domain size can be increased or decreased on demand by any predefined amount). In the present simulations, an expansion front nearing the computational domain boundary triggered the doubling of the domain size. Jointly with ACM, this method allows efficient and physically accurate simulations, starting from domain sizes of the order of millimeters and expanding to tens of centimeters, as needed. This combination of techniques thus makes the computation fast over the large distances from the target (>1–10 cm) and over long times (t > 1–10 μs). The vacuum background is characterized by a fixed number density of neutral species, set to match the experimental conditions described in the next section.

Measurements of laser plasmas analogous to the simulation were carried out inside a vacuum chamber with a background pressure of 3 × 10⁻⁶ torr. The focused beam of a KrF excimer laser (Lambda Physik LPX 305i), which has pulse duration of ∼25 ns and wavelength of 248 nm, was used to ablate solid metal targets at a 45° angle of incidence. The laser spot size on the target was adjustable using rectangular apertures that blocked the periphery of the excimer laser beam. Apertures were placed in the beam path between the laser output window and the focusing lens. For all apertures used, the spot area on the target was approximately rectangular with a 3:1 aspect ratio. A 10 × 10 mm square, planar ion probe consisting of a thin molybdenum foil backed by an alumina plank was placed at a distance d = 4.0 cm from the target surface, along the central axis of the plume. The probe was biased with −70 V, and the current resulting from the passage of the plasma pulse was determined from the voltage drop across a 10 Ω resistor, connected to ground.

Assuming quasi-neutrality conditions, the saturation current pulse I(t) measured by the probe under strong negative bias, can be used to infer the electron density n(t) at the probe location by [73]

$$\begin{equation}n(t)=\frac{I(t)}{vAe},\end{equation} \tag{ 20 }$$

where t is the time since the firing of the laser (i.e. the TOF of the fluid element being measured), v is the plasma flow velocity, A is the probe collection area, and e is the electron charge. A practical approximation for n(t) is obtained by using v = d/t in equation (20), where d is the distance between the ablation spot and the probe.

Se and Te ablation targets were synthesized by pressing elemental powders (Alfa Aesar, 99.999%) into 20 mm diameter discs. Each disc was sealed in an evacuated quartz ampoule (<10⁻³ torr) and sintered for 24 h at 80% of the melting temperature. The surface of each resulting high density, sintered target was polished before laser ablation. The Cu target was a nominally 99.999% pure one-inch diameter disc subjected to the same polishing step as the chalcogen targets.

The electron density at a fixed position in space provides general characteristics of the plasma expansion. Its time dependence, n(t), can be easily extracted from simulations and experiments. Figure 1(a) shows simulated n(t) for ablation of Te. The n(t) traces are for a distance d = 4.0 cm from the ablation spot, along the central axis of the plume. Simulations were run using a typical C_eff value. The curves show two peaks, with relative amplitudes that vary with the laser spot size. For the smallest spot size, a fast, high-density peak leads the expansion. A slower and broader peak of reduced density lags behind. As the spot size increases, the slow peak becomes dominant. As seen in figure 1(b), the same trend is detected by the Langmuir probe in the experimental plasma—albeit over a broader range of spot size variation. Traces with comparable characteristics are also seen for ablation of Se (not shown). The similarity in trends between simulation and experiment indicates that the effective model can reproduce essential features of observed laser plasmas. Quantitative correspondence between simulation and experiment is not expected here, because the value of C_eff has yet to be optimized as it will be explained in section 4.2—figure 1(a) represents simulations unconstrained by experiments. A sensitivity study of the impact of C_eff on the simulation (figure 1(c)) shows that the overall correspondence between simulation and experiment is preserved when C_eff is varied in the typical range expected to produce constrained simulations (0.4 × 10⁻²³ – 4.0 × 10⁻²³ m²). It is particularly meaningful that the 2D-axisymmetric simulation may allow evaluation of the effect of laser spot size, which has not been widely studied, either experimentally [65] or computationally [31, 63]. Spot size effects cannot be assessed in 1D simulations.

Zoom In Zoom Out Reset image size

A simulation can be constrained by determining the value of C_eff that causes the predicted n(t) to be consistent with its experimental counterpart. For constraining considerations it is instructive to convert the variable t into an energy axis defined by ɛ = m(d/t)²/2, where m is the mass of the corresponding ion. The quantity ɛ is often referred to as the 'ion TOF kinetic energy' [37].

Figure 2(a) shows experimental traces of n(ɛ) for Se, Te, and Cu under the irradiation conditions illustrated in the inset. The Se plasma is dominated by ions with kinetic energy peaking at ∼7.6 eV, while a faster sub-population appears with energy in the 25–40 eV range—a tail of still higher energies extends up to ∼52 eV. The peaks for Te are broader and have greater peak energies compared to Se. The slow Te component peaks at ∼11 eV while the fast ions have maximum density at 79 eV. The peak ratios of the slow to the fast components are ∼1.8 for Se and 3.6 for Te. The n(ɛ) curve for Cu is centered at ∼19.1 eV.

Zoom In Zoom Out Reset image size

Simulated n(ɛ) curves for approximately equivalent conditions are shown in figure 2(b). Simulations using C_eff = 0 lead to overall densities that are much lower than experiments, with values for Se and Cu off by more than one order of magnitude. In addition, the two-peak structure is noted to be inverted for Te. Figure 2(b) also shows simulated n(ɛ) using values of C_eff that cause the maximum of each curve to match the maximum of the corresponding experimental trace. The best matches yield the values of C_eff shown in the figure and in table 2. While the width and exact peak locations retain substantial differences between simulation and experiment, we note the emergence of a two-peak structure that appears equivalent to the experimental chalcogen data. The peak ratios for both, Se and Te, evaluate to 1.6. This outcome matches the experimental peak ratio for Se, but is a factor of ∼2.2 lower than the Te measurement. The constrained Cu simulation produces a very narrow n(ɛ) curve, that is centered at 19.3 eV. This is <1% different from the measured one. In the next section, we describe the simulation predictions when the C_eff values determined by this constraining exercise are used throughout the calculations.

It is important to note that in applying the model to Se, Te, and Cu in this paper, we can only assert a partial equivalence between experimental and simulation environments. This is because of experimental limitations and computational symmetry requirements. The experimental laser spot is roughly rectangular due to the excimer laser electrode discharge geometry. The simulations, on the other hand, are axially symmetric, resulting in a circular laser spot on the target. The rectangular experimental spot implies a plume that deviates from axial symmetry. The amount of forward-peaking of the plume differs when viewed on the plane of the long side of the rectangle compared to the short side. By matching the simulation spot diameter to the short side of the experimental rectangle—as illustrated in the insets of figure 2—we are likely simulating a lower bound for the expansion velocities along the plume axis. Constraining the simulation by n(ɛ) from more symmetric experimental conditions, may allow updated C_eff values for improved predictions of overall plume behavior. It is also worth mentioning that the asymmetry of the experimental plume may be partially responsible for the broader experimental n(ɛ) curves compared to the simulation. Moreover, one can envision a variety of alternative methods to constrain the model and determine C_eff. Instead of the experimental n(ɛ), one could use, for example, the time-dependent angular distribution of emitted species, n(θ, t). Because n(θ, t) contains 3D information on the gas dynamics problem, it could allow simulations with significantly improved predictive capabilities. It could also enable testing of more sophisticated effective terms that include a time-dependent C_eff, additional nonlinear absorption effects, and mass removal processes that can better approximate the variety of phase transformations expected in the target material [28, 70]. These studies are underway in our group and will provide further opportunities to evaluate the range of applicability of the effective model approach.

Full scale simulations employing the C_eff values of table 2 can be conveniently described by dividing the temporal evolution of the plasma into three stages: laser-target interaction, early expansion and ionization, and long term expansion out to centimeter distances.

The laser-target interaction and early expansion of the plumes (expansion distances much less than the spot size) exhibit very low sensitivity to the dimensionality of the description. 2D-axisymmetric and 1D results are essentially the same along the expansion direction. This allows us to explore these early processes at the reduced dimensionality of the 1D description with the finest spatial resolution, without the extra computational cost of multiple dimensions. Figure 3(a) shows the laser irradiance predicted to arrive at the targets for the chosen laser output fluence of 4.0 J cm⁻². In the case of Cu, very little plasma shielding occurs. Essentially all of the laser energy reaches the target throughout the duration of the pulse.

Zoom In Zoom Out Reset image size

For Se and Te, on the other hand, strong plasma shielding is predicted. For Se, the onset of plasma shielding occurs just as the laser reaches peak power density, while for Te, shielding begins at ∼60% of the peak. In both cases plasma shielding is sustained for the remainder of the pulse, with near total shielding for Te. The normal velocity of the evaporated Te atoms (figure 3(b)) peaks quickly near 850 m s⁻¹ and then decays because the fully shielded target is cooling down for the remainder of the laser pulse. The velocity of the Se and Cu atoms increases more slowly with both reaching the maximum velocity just after the onset of plasma shielding for Se and at the end of the pulse for Cu.

The distinct plasma shielding behavior of the three metal plumes stems from differences in their total plasma absorption. Figure 4 details the plume density, absorption coefficients, and plasma temperature as a function of distance from the target at 15 ns, the midpoint time of the laser pulse. The number density of neutral atoms n⁰ and electrons n_e are shown in figure 4(a).

Zoom In Zoom Out Reset image size

For the chalcogens, their lower boiling points and thermal conductivities cause the densities of evaporated neutral species at the surface of the targets to be three orders of magnitude higher than for Cu. This difference in initial density is certainly a determining factor in the subsequent plasma plume structure. Simulations using lower laser fluences of 1.2–1.5 J cm⁻² for the chalcogens result in initial densities of Se and Te neutrals in the 5 × 10²⁵–4 × 10²⁶ m⁻³ range, which are significantly closer to those predicted for Cu at the higher 4.0 J cm⁻² fluence.

It is important to note that in the context of the effective model, the high initial vapor densities near the target should be interpreted as effective densities, since they result from the admittedly simplistic description of equations (3)–(6). Significantly, the clearly overestimated initial vapor densities lead to greater laser absorption by the plasma, which quickly reduces target heating and brings the plasma densities into realistic ranges in the microsecond regime. This self-limiting coupling between material ejection from the target and plasma absorption lends itself well to an effective description, which should also consider the high surface temperatures expected as effective temperatures.

In this regard, an observation that demands careful scrutiny is that the solution of the heat transfer problem (equation (3)) leads to surface temperatures that exceed the thermodynamic critical temperature of the target materials. As a result, equation (4) can at best, only be considered an effective description of the problem. Near and above the critical temperature T_c, the liquid-to-vapor phase change described by the Clausius–Clapeyron equation does not represent the actual physical mechanism by which material from the target enters the gas phase. It is known from multiple experiments on laser ablation of metals that thermal-evaporation-only models, under conditions of moderate-to-high overheating, usually underestimate the amount of ablated material by an order of magnitude, with the typical thickness of the predicted ablated layer on the order of hundreds of nanometers, while the real ablation depth is on the order of microns. Hence, it is essential to evaluate the potential error or uncertainty introduced in our simulation due to supercritical processes not being explicitly taken into account. An estimate of this uncertainty can be obtained by comparing the amount of ablated material predicted by the effective model to the amount of ablated material that would be produced if all material overheated above 0.9T_c were removed from the target [53, 56, 78]. In the representative case of Te, the experimentally constrained effective model (C_eff = 2 × 10⁻²³ m², table 2) predicts an ablation depth of 10.7 μm (i.e. with p_vap evaluated by equation (4)). Estimates of complete removal of material above 0.9T_c indicate ablation depths between 0.48 μm and 3.3 μm. The lower bound supercritical estimate assumes all vapor forming with the critical density of Te (i.e. p_vap = ρ_c kT_c/m), whereas the upper bound would occur if the vapor formed with the solid density (i.e. p_vap = ρ_m kT_c/m). In these supercritical estimates, we used the effective model to evaluate the time-dependent recession velocity of the interface between material below and above 0.9T_c. All material above 0.9T_c was instantaneously shifted to the vapor phase. These results indicate that under the current constraining scheme, the effective model probably overestimates the amount of ablated material by factors ranging from 3 to 20. Caution should therefore be exercised when interpreting simulation outcomes if accuracy greater than what is warranted by the above margin is needed. Improvements are likely to be achieved by explicitly incorporating supercritical processes in the effective model and utilization of alternative methods to constrain the simulation as discussed in section 4.2.

The spatial dependence of the plasma absorption coefficients is shown in figure 4(b). The effective model shows early expansion dominated by α_PI and α_eff, with α_IB being of significantly less importance. Particularly high relative α_eff values emerge for the Cu case. The plasma temperatures are shown in figure 4(c). As expected, the temperature peaks in each case near the position where the plumes have the greatest α_PI and n_e. The Se plume has a peak temperature of 1.48 eV at 19 μm from the target, Te peaks at 1.40 eV at 32 μm, and Cu at 0.47 eV at 12 μm. The kinetic energy of singly-charged ions at the expansion front can be evaluated using the position of the sharp temperature increase at the plasma front and the elapsed time of 15 ns. Te has the greatest energy at 21.1 eV followed by Se and Cu at 9.68 eV and 3.02 eV, respectively.

A qualitative view of the long-term expansion of the plumes is obtained if simulations based on the 1D description are run for sufficiently long times as shown in figure 5. We note that each plume reaches 5 cm at different times, 7.68 μs, 5.97 μs, and 7.74 μs for Se, Te, and Cu, respectively. The number densities are overestimated, since the plasma cannot expand in the transverse dimensions. Overall, the chalcogen plumes appear denser than the Cu plume, consistent with the early expansion profiles of figures 3 and 4. The Se and Te expansions (figures 5(a) and (b)) show increasing trend for ionization of their interiors. This is in contrast to the Cu plume, which has meaningful concentrations of ions only in the frontal region. For all three cases, the leading edge of the plume is fully ionized, with the Te plasma being somewhat unique due to the presence of a high Te²⁺ number density at its leading edge for the used irradiation conditions.

Zoom In Zoom Out Reset image size

A richer and expressly quantitative perspective on the long-distance expansion is furnished by the 2D-axisymmetric results shown in figure 6. The distributions of neutral species for the chalcogens (figures 6(a) and (b)) are heavily weighted toward the near-target region and fall off rapidly in all directions. Cu⁰, on the other hand, is more homogeneous throughout the plume (figure 6(c)). The number density of singly charged ions is also quite distinct for Se and Te, with well defined peaks roughly located near a central bulge in the plume (figures 6(d) and (e)). This is in contrast to the Cu plume, which is dominated by Cu⁰, with Cu⁺ ions found exclusively at the leading edge of the expansion (figure 6(f)).

Zoom In Zoom Out Reset image size

Consistent with the greater degree of ionization of the interior of their plumes, the temperature profiles of figures 6(g)–(i) reveal that the chalcogen plasmas are hotter throughout. Their temperatures are approximately 0.5 eV on the front half of the plume and near 1 eV at the leading edge, where there is an appreciable concentration of doubly ionized species. The concentrations of Se²⁺ and Te²⁺ (not shown) have peak densities of 5 × 10¹⁷ m⁻³ and 1 × 10¹⁸ m⁻³, respectively. Doubly ionized species only exist in ∼1 eV temperature regions and are therefore absent in the Cu plume, which is cooler than the chalcogens, reaching its highest value of 0.38 eV only at the expansion front.

Figure 6 illustrates the key advantage of multidimensional simulations: plasma processes are allowed to couple to the translational degrees of freedom that are perpendicular to the direction normal to the target. The most obvious benefit is bringing the predicted number densities into much closer agreement with experiments. We note that although the 1D density vs distance profiles (figure 5) show good qualitative correspondence with the 2D-axisymmetric plumes along the axis of symmetry (figures 6(a)–(f)), the 1D number densities are not realistic. The maximum Te⁺ density, for example, is 3 × 10²² m⁻³ (1D) vs 2 × 10¹⁹ m⁻³ (2D-axi). The lower densities brought about by plume expansion in the directions transverse to the symmetry axis are in significantly closer agreement with experimental density values, as can be inferred from data in figure 2, for instance.

In order to evaluate the sensitivity of our results to the value of the cross section Q in equation (9), we ran Te simulations between Q/2 and 2Q, where Q = 10⁻⁴⁰ cm⁵ was the value used throughout this study, and represents an estimate of the mean value of Q for the variety of temperatures and chemical species of interest in our study [66]. For each value of Q, a new C_eff was determined. A comparison between results with optimized C_eff showed that the overall plume structure is unaltered in the Q/2–2Q range. The peak densities and spatial distribution of neutrals varied by no more than ∼11%. At expansion distances ∼6 cm, the variations in electron and ion densities were less than ∼8%. In expansions exceeding ∼10 cm, electron and ion densities were essentially insensitive to Q. In all cases, the velocity of the ionized leading edge was also unchanged in the Q/2–2Q range. Overall, these variations indicate that the effective model can accommodate to a good extent the uncertainties associated with the value of Q.

The gains of the effective 2D-axisymmetric description allow quantitative studies of changes in the geometry of the expansion, which may vary from the quasi-spherical character of Knudsen-layer-only processes (Mach number ∼1) to the highly forward-peaked profiles of hypersonic expansions [51]. The spatiotemporal dependence of the particle flux, the angular distribution of plasma species, and the effect of the laser spot size, can all be easily derived from these results to help guide experiments in materials growth and processing.

It is also of note that the effective model and our simulation framework are compatible with interaction of multiple plasma plumes. This is significant because there are many opportunities for exploring plasma phenomena triggered by the collision of two or more laser-produced plasmas. Changes in plasma shielding behavior and the ability to detect broadening as well as the formation of new shocks due to plume interaction are already built into our method. Furthermore, the repetitive nature of the pulsed process raises the question of interaction between subsequent pulses. In materials synthesis applications, where repetition rates are usually <50 Hz, ions and neutrals dissipate in a timescale shorter than the next pulse. It is known, however, that molecular, cluster, and nanoparticle constituents may remain in the gas phase for periods of time in excess of 1 s [20]. During actual deposition, subsequent pulses will interact with these remaining materials. Further developments of our model to include these larger constituents may suggest new directions for the use of laser-based plasma technologies in materials research.