format: text

text:
# SEMILEPTONIC 𝐵 DECAYS Recent Results from LEP and Comparison with Υ(4S) Data

## 1 Introduction

The subject of this review is the status of the studies on semileptonic (s.l.) $`B`$ decays. These decays represent a favorable laboratory to study the dynamics of heavy quark decays, to determine the size of the $`|V_{ub}|`$ and $`|V_{cb}|`$ CKM mixing matrix elements, crucial in the unitarity triangle test of the Standard Model, and to acquire informations on the structure of the $`B`$ meson itself.

The Aleph, Delphi, L3 and Opal experiments have analyzed together about 1.5 M s.l. $`B`$ decays recorded at Lep, at energies around the $`Z^0`$ pole from 1990 to 1995. This statistics is significantly lower than about 4 M s.l. decays recorded by Cleo and 7.5 M already logged by the BaBar and Belle experiments. However, the significant boost of the beauty hadrons, the confinement of their decay products in well separated jets and the production of almost all beauty hadron species make it possible to perform topological decay reconstruction with good efficiency on the Lep data sets, thus exploiting analysis techniques complementary to those used at lower energy $`e^+e^{}`$ colliders. Since operator product expansion (OPE) predictions apply to sufficiently inclusive observables, it is advantageous to reconstruct s.l. decays inclusively, with only mild cuts on the lepton and hadron energies. Further, the Lep kinematics allow the two extreme kinematical regions at small and large momentum transfer to be accessed, since the $`B`$ decay products gain enough energy from the $`B`$ boost.

## 2 Recent Results

### 2.1 Inclusive s.l. Branching Fraction

The determination of the inclusive s.l. branching fraction BR($`bX\mathrm{}\overline{\nu }`$) is important for the measurement of $`|V_{cb}|`$. Together with the measurement of the charm multiplicity and of exclusive decays discussed later in this Section, they provide with a test of the s.l. decay width. The main experimental issue here is the separation of the prompt $`b\mathrm{}`$ signal from the cascade $`bc(\overline{c})\overline{\mathrm{}}(\mathrm{})`$ and the $`c\overline{\mathrm{}}`$ backgrounds. The Lep analyses use the charge correlation between the $`b`$ and the lepton, the decay topology, and double tagged events where both beauty hadrons decay semileptonically. The dominant source of uncertainty is due to the modeling of the $`b\mathrm{}`$ and $`c\overline{\mathrm{}}`$ spectra. The Lep average, obtained from the direct determinations, gives BR($`bX\mathrm{}\overline{\nu }`$) = 0.1056 $`\pm `$ 0.0011 (stat) $`\pm `$ 0.0018 (syst). After having rescaled this value by $`\frac{1/2(\tau (B_d)+\tau (B_u))}{\tau (b)}=1.021\pm 0.013`$, to account for the different beauty hadron species produced, this result can be compared with the most recent determination at the $`\mathrm{{\rm Y}}(4S)`$, obtained by Cleo using a lepton tagged method to separate the prompt lepton yield in $`B`$ decays from the backgrounds giving BR($`bX\mathrm{}\overline{\nu }`$) = 0.1049 $`\pm `$ 0.0017 (stat) $`\pm `$ 0.0043 (syst).

It is useful to analyze these results in relation with the average number of charm particles in $`B`$ decays, $`n_c`$, since the semileptonic, double charmed and charmless yields are correlated once the total decay width is fixed. The average $`n_c`$ values, obtained with three independent methods, are sumarised in Table 1.

A recent by Sld based on the extraction of the double charm yield from the $`B`$ decay topology is dicussed in these proceedings.

The scaled Lep+Sld averages and the Cleo values can be independently compared with theoretical predictions, obtained using heavy quark expansion to order $`1/m_b^2`$. These are shown in Figure 1 for a range of values of the renormalization scale $`\mu `$ and of the ratio of the quark pole masses $`m_c/m_b`$. The experimental results agree with these predictions at $`\mu 0.25m_b`$ and $`m_c/m_b0.33`$.

### 2.2 Charmless s.l. Branching Fraction

An accurate determination of the charmless s.l. $`B`$ branching fraction is important for the measurement of the $`|V_{ub}|`$ element. This has represented a significant challenge to both theorists and experiments since Cleo first observed charmless s.l. decays as an excess of leptons at energies above the kinematical limit for decays with an accompanying charm hadron. However, the limited fraction ($``$ 10%) of charmless decays populating this end-point region results in a significant model uncertainty for the extraction of $`|V_{ub}|`$. An analysis technique based on the invariant mass $`M_X`$ of the hadronic system recoiling against the lepton pair, peaked for $`bX_u\mathrm{}\overline{\nu }`$ at a significantly lower value than for $`bX_c\mathrm{}\overline{\nu }`$, was proposed several years ago and it has been the subject of new theoretical calculations. If $`bu`$ transitions can be discriminated from the dominant $`bc`$ background up to $`M_XM(D)`$, this method is sensitive to $``$ 80% of the charmless s.l. $`B`$ decay rate. Further, if no preferential weight is given to low mass states in the event selection, the non-perturbative effects are expected to be small and the OPE description of the transition has been shown to be accurate away from the resonance region. The requirement to isolate the $`bu`$ contribution to the s.l. yield from the $``$ 60 times larger $`bc`$ one, while ensuring an uniform sampling of the decay phase space to avoid biases towards the few exclusive low-mass, low-multiplicity states, make this analysis a major experimental challenge.

Aleph, Delphi and L3 have developed new analysis techniques based on the observation that $`bX_u\mathrm{}\overline{\nu }`$ decays can be inclusively selected from $`bX_c\mathrm{}\overline{\nu }`$ by the difference in the invariant mass and kaon content of the secondary hadronic system and in the decay multiplicity and vertex topology. These features have been implemented differently in the analyses by the three experiments: Aleph used a neural net discriminant based on kinematical variables, Delphi preferred a classification of s.l. decays on the basis of their reconstructed hadronic mass $`M_X`$, decay topology and presence of secondary kaons and L3 adopted a sequential cut analysis based on the kinematics of the two leading hadrons in the same hemisphere as the tagged lepton. Starting from a natural signal-to-background ratio S/B of about 0.02, Aleph obtained S/B = 0.07 with an efficiency $`ϵ=`$ 11%, Delphi had S/B = 0.10 with $`ϵ=`$ 6.5% and L3 had S/B = 0.16 with $`ϵ=`$ 1.5%.

All three experiments observed a significant data excess over the estimated backgrounds corresponding to 303 $`\pm `$ 88 events for Aleph, 214 $`\pm `$ 56 for Delphi and 81 $`\pm `$ 25 for L3. The characteristics of these events correspond to those expected for $`bX_u\mathrm{}\overline{\nu }`$ decays (see Figure 2). The inclusive charmless s.l. branching ratios summarized in Table 2 were obtained. The Lep average value is BR($`bX_u\mathrm{}\overline{\nu }`$) = (1.74 $`\pm `$ 0.37 (stat.+exp.) $`\pm `$ 0.38 ($`bc`$) $`\pm `$ 0.21 ($`bu`$)) $`\times `$ 10-3 = (1.74 $`\pm `$ 0.57) $`\times `$ 10-3. The Lep experiments have shown the feasibility of these inclusive analyses, due to the favorable kinematics and the decay reconstruction capabilities of their detectors. While more precise data on $`B`$ and $`D`$ decays will decrease the dominant $`bc`$ systematics of this measurement, future perspectives for inclusive charmless s.l. rate determinations are not yet clearly outlined. The hadronic mass analysis puts even further problems to symmetric and asymmetric $`B`$ factories at the $`\mathrm{{\rm Y}}(4S)`$ peak, due to the confusion between the decay products of the $`B`$ and $`\overline{B}`$ and of the reduced $`B`$ decay length, compared with Lep. A feasibility study, performed by BaBar, requires full reconstruction of one $`B`$ meson through an exclusive decay mode to solve the first problem, and predicts a reconstructed signal sample with $`M_X<`$ 1.7 GeV/$`c^2`$ of about 100 to 300 events for the full data statistics of 100 fb-1. A new technique based on the reconstruction of the di-lepton $`\mathrm{}\overline{\nu }`$ invariant mass of the decay has been proposed recently. At low-energy $`B`$ factories, the di-lepton mass measurement will profit from the $`\nu `$ reconstruction techniques that rely on the $`B`$ production at threshold, already successfully exploited by Cleo.

### 2.3 Exclusive s.l. $`B`$ Meson Decays

Exclusive s.l. $`B`$ decays have been studied at the $`\mathrm{{\rm Y}}(4S)`$ and at Lep, to establish the relative contribution of the individual channels to the s.l. decay width and to extract the relevant CKM elements. Charmless $`B\pi \mathrm{}\nu `$ and $`\rho \mathrm{}\nu `$ decays have been observed by Cleo and their rates measured.

The decay $`\overline{B_d^0}D^+\mathrm{}^{}\overline{\nu }`$ has received specific attention for the extraction of $`|V_{cb}|`$ from a study of its rate $`\frac{d\mathrm{\Gamma }}{dw}=K(w)F^2(w)|V_{cb}|^2`$ as a function of the the $`D^{}`$ and $`B_d`$ four-velocity product $`w`$. In this expression $`K(w)`$ is a phase space factor and $`F(w)`$ the hadronic form factor. The reconstruction of the decay can be performed fully exclusively, through the decay $`D^+D^0\pi ^+`$ followed by $`D^0K^{}\pi ^+`$ and, at Lep, also by partial inclusive reconstruction of the $`D`$ meson resulting in a significant increase in the selection efficiency. At Lep, decays can be efficiently reconstructed closer to zero $`D^+`$ recoil energy than at the $`\mathrm{{\rm Y}}(4S)`$. But using exclusive reconstruction and at the $`\mathrm{{\rm Y}}(4S)`$, where the $`B`$ meson is almost at rest, a better resolution on $`w`$ is achieved: Cleo obtains $`\sigma (w)`$ = 0.03 compared with $`\sigma (w)`$ = 0.07 to 0.12 from the inclusive Lep analyses (see Figure 3).

The main background that these analyses have to reduce and understand is due to s.l. transitions into charmed excited states $`D^{}`$ producing a $`D^{}`$ in their decay as discussed below. The extrapolation of the form factor $`F(w)`$ is based on a dispersion-relation parameterization. By combining the measurements by Aleph, Delphi and Opal, the Lep averages $`F(1)|V_{cb}|=(34.5\pm 0.7((\mathrm{stat}))\pm 1.5((\mathrm{syst}))\times 10^3`$ and $`\rho ^2=1.13\pm 0.08((\mathrm{stat}))\pm 0.15((\mathrm{syst}))`$ were obtained with a fit confidence level of 12% (see Figure 4). Cleo recently reported $`F(1)|V_{cb}|=(42.4\pm 1.8((\mathrm{stat}))\pm 1.9((\mathrm{syst}))\times 10^3`$ and $`\rho ^2=1.67\pm 0.11((\mathrm{stat}))\pm 0.22((\mathrm{syst}))`$ that gives 2.4 $`\sigma `$ higher extrapolated rate at $`w=1`$ with a larger slope. Since the intercept at $`w=1`$ and the slope are highly correlated, this prompts further investigation of this decay through new analyses of the experimental data and of the underlying uncertainties.

The study of s.l. $`B`$ decays into orbitally and radially excited charmed mesons has established the production of the narrow orbital excitation $`D_1(2420)`$, while the $`D_2^{}(2460)`$ and a broad state decaying into $`D^+\pi ^{}`$ have been observed by Cleo in hadronic $`B`$ decays. A recent analysis by Delphi has separately measured the narrow $`D_1(2420)`$ and broad, or non-resonant, $`D^{()}\pi `$ production (see Figure 5). Results are summarized in Table 3.

The $`D_2^{}(2460)`$ production rate is less than that of $`D_1(2420)`$: $`R^{}=\frac{BR(\overline{B}D_2\mathrm{}^{}\overline{\nu })}{BR(\overline{B}D_1\mathrm{}^{}\overline{\nu })}<0.6`$ in disagreement with the HQET prediction in the infinitely heavy quark mass limit $`R^{}1.6`$, thus implying large $`1/m_Q`$ corrections. Including these corrections restores agreement with the data for $`R^{}`$, however the overall yield of the broad $`D_0^{}+D_1^{}`$ states is not expected to exceed that for $`D_1(2420)`$. This apparent discrepancy with the preliminary Delphi result can now be explained as an experimental effect, if the feed-through of non-resonant $`D^{()}\pi `$ production is non negligible.

A first estimate of the slope of the $`\mathrm{\Lambda }_b`$ form factor has also been obtained at Lep. Delphi performed a fit to the reconstructed $`w`$ distribution and event rate for a sample of candidate $`\mathrm{\Lambda }_b\mathrm{\Lambda }_c\mathrm{}\overline{\nu }`$ decays obtaining $`\rho ^2`$ = 1.55 $`\pm `$ 0.60 (stat) $`\pm `$ 0.55 (syst), which is within the range of the current theoretical predictions.

## 3 Extraction of $`|V_{ub}|`$ and $`|V_{cb}|`$

### 3.1 $`|V_{ub}|`$

The value of the $`|V_{ub}|`$ element can be extracted from the inclusive charmless s.l. branching fraction BR($`bX_u\mathrm{}\overline{\nu }`$) by using the following relationship derived in the context of Heavy Quark Expansion:

$$|V_{ub}|=.00445(\frac{\mathrm{BR}(bX_u\mathrm{}\overline{\nu })}{0.002}\frac{1.55\mathrm{ps}}{\tau (b)})^{\frac{1}{2}}\times (1\pm .020(\mathrm{QCD})\pm .035(\mathrm{m}_\mathrm{b}))$$
(1)
where the value $`m_b`$ = (4.58 $`\pm `$ 0.06) GeV/c$`^2`$ has been assumed. The theoretical uncertainties are small, due to the absence of $`1/m_b`$ term in the expansion and of $`1/(m_bm_c)`$ dependence, and it is dominated by that on the $`b`$ mass. Inserting the Lep average BR($`bX_u\mathrm{}\overline{\nu }`$), $`|V_{ub}|`$ was determined to be: $`|V_{ub}|=(4.13_{0.47}^{+0.42}(\mathrm{stat}.+\mathrm{det}.)_{0.48}^{+0.43}(\mathrm{b}\mathrm{c}\mathrm{syst}.)`$ $`{}_{0.25}{}^{+0.24}(\mathrm{b}\mathrm{u}\mathrm{sys}.)`$ $`\pm 0.02(\tau (b))\pm 0.20(\mathrm{HQE}))\times 10^3`$. The $`bu`$ and HQE model systematics is slightly below 10% and the large uncertainties from the modeling of $`bc`$ decays can be reduced in future by more precise data. This inclusive determination of $`|V_{ub}|`$ can be compared to that extracted from the determination of the exclusive rate for the decay $`B\rho \mathrm{}\overline{\nu }`$ by Cleo giving $`|V_{ub}|`$ = (3.25 $`\pm `$ 0.14 (stat.) $`{}_{0.29}{}^{}{}_{}{}^{+0.21}`$ (syst.) $`\pm `$ 0.55 (model) ) $`\times 10^3`$. The two measurements are consistent within their uncertainties, which are mostly uncorrelated. It is expected that the large model dependence of the exclusive method will be reduced by computing the hadronic form factor from lattice QCD.

### 3.2 $`|V_{cb}|`$

There are two methods, to extract the $`|V_{cb}|`$ element from s.l. decays: i) from the inclusive s.l. $`B`$ decay width, once the charmless contribution has been subtracted, and ii) by the extrapolation of the rate for exclusive decays, as $`BD^{}\mathrm{}\overline{\nu }`$, at zero recoil. The two results can be used as a check of the underlying theory and to improve the $`|V_{cb}|`$ accuracy by their averaging, the systematic uncertainties being partially uncorrelated.

For the inclusive method the relationship:

$$|V_{cb}|=.0411\times (\frac{\mathrm{BR}(bX_c\mathrm{}\nu )}{.105}\frac{1.55\mathrm{ps}}{\tau (b)})^{\frac{1}{2}}\times $$
(2)
$$(10.024\times \left(\frac{\mu _\pi ^20.5}{0.1}\right))\times (1\pm 0.030(\mathrm{pert}.)\pm 0.020(\mathrm{m}_\mathrm{b})\pm 0.024(1/\mathrm{m}_\mathrm{b}^3))$$

can be used to extract $`|V_{cb}|`$. By taking the inclusive and charmless s.l. branching fractions presented above and the average $`b`$ lifetime $`\tau (b)`$ = (1.564 $`\pm `$ 0.014) ps, the result is $`|V_{cb}|`$ = (40.7 $`\pm `$ 0.5 (exp) $`\pm `$ 2.0 (th)) $`\times 10^3`$. This inclusive method is limited by the lepton spectrum model systematics from the s.l. branching ratio and by the theory systematics due to the values of $`m_b`$ and $`\mu _\pi ^2`$ (or $`\lambda _1`$) and to the use of only the first terms of the operator product expansion as highlighted in Eq. (2).

In the exclusive method, the quantity $`|V_{cb}|\times F(1)`$ is determined from the differential decay rate extrapolated to zero recoil. The extraction of $`|V_{cb}|`$ relies on heavy quark symmetry for the form factor normalization $`F(1)=1`$ in the heavy quark limit and requires the computation of the finite $`b`$-quark mass effect and of non-perturbative QCD corrections. The value $`F(1)=0.8800.024\frac{\mu _\pi ^20.5}{0.1}\pm 0.035(\mathrm{excit})\pm 0.010(\mathrm{pert})\pm 0.025(1/m_b^3)`$ has been adopted by the Lep working group giving $`|V_{cb}|`$ = (39.8 $`\pm `$ 1.8 (exp) $`\pm `$ 2.2 (th)) $`\times 10^3`$. This result is in agreement with that obtained with the inclusive method. There has been a recent lattice determination of $`F(1)`$; unquenched computations will provide smaller and better understood systematics. However the larger value obtained by the recent Cleo analysis suggests a cautious attitude and the need for more experimental data.

## 4 Open Questions and Model Systematics

As more data on s.l. $`b`$ decays have become available and the analysis techniques improved, the determinations of basic parameters of $`b`$ decays, such as the $`|V_{ub}|`$ and $`|V_{cb}|`$ elements, are limited by theoretical uncertainties and by the modelling of the signal and backgrounds. While this is promoting new phenomenological approaches, it also requires new measurements to better determine the input parameters and further constrain the models. This interplay between progresses in theory and new experimental insights addresses three main classes of questions: i) the definition of a coherent method to include bound state effects in the description of inclusive decays, ii) the accurate determination of the fundamental parameters, $`m_b`$, $`m_bm_c`$ and $`p_b^2`$ (or $`\lambda _1`$), and iii) the estimate of the effects of the missing terms in the OPE and of violations of quark-hadron duality.

At present, inclusive spectra are predicted by a variety of specialized models ranging from fully inclusive models, such as the ACCMM model, to those saturating the inclusive decay width by the contribution of several exclusive final states, like the ISGW and ISGW-2 models. While these models can describe the data quite precisely, after tuning of their parameters, their application tout court to different processes is often not justifiable. In inclusive models, the spectra are obtained for a free quark decaying into partons and the non-perturbative corrections are included by convoluting the parton spectra with a function encoding the kinematics of the $`b`$ and spectator $`\overline{q}`$ quarks inside the heavy hadron. The model free parameters are derived from the shapes of inclusive observables.

Figure 6 summarizes the results from fits to the energy spectrum of the lepton in s.l. $`bX\mathrm{}\overline{\nu }`$ decays from Lep and Cleo data, showing a good consistency. The main issue here is to establish the validity of the ACCMM model, with these fitted values, to other decays, such as $`bX_u\mathrm{}\overline{\nu }`$ and $`bs\gamma `$.

This difficulty may be overcome by a description of the Fermi motion in the framework of QCD. This introduces an universal shape function, $`f(k_+)`$. At leading order and in the large $`m_b`$ limit, the light-cone residual momentum $`k_+`$ can be expressed as the difference between the $`b`$ quark mass and its effective mass $`m_b^{}`$ inside the hadron: $`m_b^{}=m_b+k_+`$. However, the functional form of the function $`f(k_+)`$ is not a priori known, except for its first three moments. The effect of different ansatz for the description of the Fermi motion of the $`b`$ quark inside the heavy hadron has been studied in the case of $`bX_u\mathrm{}\overline{\nu }`$ both at the parton level and for the physical observables after full detector simulation. It was found that the Fermi motion description contributes with an uncertainty on the branching fraction determination of $`\pm `$ 5 - 15%, increasing as the hadronic mass $`M_X`$ cut to distinguish $`bu`$ from $`bc`$ transitions is lowered and depending also on the other selection criteria.

Since $`m_b`$ enters at the fifth power in the expression of the decay width, $`\mathrm{\Gamma }m_b^5|V_{CKM}|^2`$, it is important to determine its value with an accuracy of better than $`\pm `$ 100 MeV/$`c^2`$ to guarantee a few percent error contribution in the extraction of $`|V_{ub}|`$ and $`|V_{cb}|`$ from inclusive decays. The dependence of the shape of inclusive spectra on the $`b`$ quark mass represents a further source of systematics from $`m_b`$. Recent estimates from $`\mathrm{{\rm Y}}`$ spectroscopy, QCD sum rules and NNLO unquenched lattice computations have shown a remarkable agreement. On the basis of these results, the Lep working groups have adopted $`m_b(1\mathrm{G}\mathrm{e}\mathrm{V})`$ = (4.58 $`\pm `$ 0.06) GeV/$`c^2`$, where $`m_b(\mu )`$ is defined such as $`\frac{dm_b(\mu )}{d\mu }=\frac{16}{9}\frac{\alpha _s(\mu )}{\pi }+\mathrm{}`$ and the quoted uncertainty defines a 68% confidence region. The first moments of the lepton energy and hadronic mass spectra can be used to determine the parameters $`\overline{\mathrm{\Lambda }}=m_Bm_b+\mathrm{}`$ and $`\lambda _1`$. An analysis of the Cleo lepton energy spectrum limited to the region $`E_{\mathrm{}}>`$ 1.5 GeV gave $`\overline{\mathrm{\Lambda }}`$ = (0.39 $`\pm `$ 0.11) GeV and $`\lambda _1`$ = (0.19 $`\pm `$ 0.10) GeV2. Cleo reported the preliminary results from the first combined study of the first two moments of the hadronic mass and the lepton energy spectra. The hadronic mass moments were found to be in good agreement with the previous result, while the lepton energy spectrum gave unlikely and incompatible values. However, the model dependence in the extrapolation to the full spectrum and the unknown contributions of higher order terms $`1/m^3`$ prevent the derivation of any conclusions from these disagreements and highlight the importance of further data.

The predictions for the inclusive analyses presented before and the extraction of the CKM matrix elements both rely on the basic assumption of duality, i.e. that the rates computed at the parton level correspond to those for the physical final states, after integrating enough hadronic channels. The validity of this assumption has been studied in a QCD model at the large $`N_c`$ and small velocity limit and in the (1+1) dimension ’t Hooft model. It was found that the sum over exclusive channels corresponds to the inclusive OPE prediction to a very good accuracy soon after having integrated the first resonant states. It is interesting to comment here on the consistency of the determinations of $`|V_{ub}|`$ and $`|V_{cb}|`$ with inclusive and exclusive methods discussed above. The agreement between the measured values is a test of the accuracy of the theoretical methods and also of the validity of the quark-hadron duality assumption, up to the combined measurement uncertainties, i.e. $``$ 8% for $`bc`$ and 25% for $`bu`$. With improved analyses and computational techniques and using the larger data sets, already becoming available at Cleo III, BaBar and Belle, these tests may reach a sensitivity of $`510\%`$ within a few years.

## 5 Acknowledgements

I wish to express my gratitude to U. Aglietti, I. Bigi, E. Barberio, G. Buchalla, G. Martinelli, M. Neubert, P. Roudeau, N. Uraltsev and W. Venus for their suggestions. The Lep results presented in this review have been provided by the LEP Working Groups on Electroweak Physics, $`B`$ lifetimes, $`V_{ub}`$ and $`V_{cb}`$.


text:
# Molecular-dynamics thermal annealing model of laser ablation of silicon

## I Introduction

The use of photons for the controlled removal of matter, i.e., laser ablation,chase94 ; haglund98  has given rise to a large number of applications for materials processing and growth; in microelectronics, examples include surface micromachining, surface cleaning, and the pulsed laser deposition (PLD)chrisey  of thin films, which allows a wide variety of materials to be grown.chrisey ; belouet96  In spite of this, a complete picture of the mechanisms underlying ablation for a broad range of pulse durations and wavelengths, and for such diverse materials as metals, polymers and semiconductors, is still lacking.

The great complexity of phenomena involved in laser ablation, determined by the coupling of the laser pulse with the optical, elastic and thermal properties of the material, has been an obstacle in understanding. On a mesoscopic scale, the processes brought about by ultrashort — femtosecond — pulses have seldom been addressed.linde98 ; linde00 ; schmidt00  On the other hand, evidence of a competition between photomechanicaljacques95 ; zhigilei00  and photothermalmiotello95 ; zhigilei00  mechanisms has begun to emerge in the light-induced removal of matter with picosecond and nanosecond pulses. In the stress confinement regime,jacques95 ; zhigilei00  ejection of matter proceeds from photomechanical effects: energy is deposited in the material in a time shorter than that needed for the generation of acoustic waves, and ablation is induced by the relaxation of the pressure gradient that builds up within the absorbing volume. In the thermal confinement regime,zhigilei00  energy is supplied to the irradiated volume at a rate faster than that at which it can be pumped out by thermal diffusion. Consequently, the material is rapidly brought near its critical point as matter removal is driven by photothermal effects:miotello95  the metastable superheated liquid undergoes a prompt transition to a gaseous phase consisting, to a large extent, of molten micron-sized droplets. This mechanism is known as explosive boiling or phase explosion;miotello95  other photothermal processes include vaporization and normal boiling (see Ref. kelly96,  and references therein). Finally, there have been repeated reportsdabby72 ; bhattacharya91 ; singh90 ; craciun98 ; craciun97  of a fourth process known as subsurface superheating (SSSH). In this scenario, heating of the material causes evaporation as a result of which a steep temperature gradient develops beneath the surface. Matter removal then takes place in a manner similar to a phase explosion following the relaxation of the pressure build-up in the bulk. It is not clear whether this mechanism is real as conflicting views can be found in the literature.miotello95 ; kelly96 ; craciun97 ; craciun98

While it is generally assumed that the laser sputtering of metal surfaces involves phase explosion mechanisms (see, e.g., Ref. bennett95, ), both photomechanical and photothermal effects have been shown to be operative in organic solids.zhigilei00  In semiconductor materials such as silicon, however, understanding of matter removal processes induced by picosecond pulses has yet to come. In the present study, we give evidence of a SSSH effect in the laser ablation of silicon. This is achieved by carrying out detailed molecular-dynamics (MD) simulations of the interaction of a picosecond laser pulse with a Si(100) substrate, explicitly accounting for carrier-phonon scattering and carrier diffusion. It is assumed that the thermal annealing model (TAM)wood84  is valid so that thermal processes only are operative. In our computer code, ‘MADTAM’, we use the Stillinger-Weber (SW)stillinger85  potential to model the atomic interactions; preliminary results from an earlier version of our model can be found in Refs. lorazo00,  and lorazo00\_2, .

Anticipating our results, we find that, below the (predicted) threshold fluence for ablation, $`F_{th}`$, the solid restores internal equilibrium by the generation of pressure waves. Above $`F_{th}`$, our simulations reveal that matter removal is triggered by SSSH effects: intense heating of the material leads to the thermal confinement of the laser-deposited energy and a steep temperature gradient beneath the surface. Interestingly, the latter is not caused by the evaporation of atoms from the outer surface but, rather, originates from the process of carrier diffusion into the bulk. As a result, a strong pressure gradient builds up $`30\text{nm}`$ below the surface as the material is overheated up to a temperature corresponding to the critical point of silicon. Ablation is subsequently initiated by the relaxation of the pressure gradient: large pieces — several atomic layer thick — of molten material are expelled from the surface, a consequence of the nucleation of voids in the bulk. Before discussing these results, however, we give a detailed description of our model.

## II The model

The processing of silicon surfaces with femtosecond laser pulses involves ultrafast, non-thermal, phenomena which cannot be accounted for by the SW potential. For this reason, we remain within the picosecond regime where thermal mechanisms are involved; this is discussed in Sec. II.1. A description of the computational method, i.e., molecular dynamics, follows in Sec. II.2. Finally, the generation of electron-hole (e-h) pairs by photons, as well as their relaxation through carrier-phonon scattering and carrier diffusion, depicted in Fig. 1, are detailed in Sec. II.3 to II.5.

### II.1 Preliminary considerations

Various computational techniques allow the interaction of light with matter to be simulated.lietoila82 ; vandriel87 ; pronko95 ; seifert98 ; garrison85 ; vertes95 ; campbell98 ; zhigilei97 ; zhigilei98 ; zhigilei00 ; lorazo00  Among them, MD appears to be a most suitable approach for addressing the numerous, complex processes involved in the laser ablation and desorption of matter: no hypotheses are needed to account for the various phenomena and atomic motion can be followed, in real time, provided that the interactions between the atoms are correctly modeled.allen87 ; smith97  However, because laser ablation involves large length and time scales, and thus a great computational effort, attempts involving classical MD are scarce.garrison85 ; vertes95 ; campbell98 ; zhigilei97 ; zhigilei97\_2 ; zhigilei98 ; zhigilei00 ; lorazo00  In order to extend the range of MD simulations, Zhigilei et al. proposed a breathing-sphere model appropriate to the study of matter removal in organic solids.zhigilei97  In this work, the phase explosion due to overheating and the laser-induced pressure build-up were identified as the key processes responsible for ablation in the regimes of thermal and stress confinements, respectively.zhigilei97 ; zhigilei00

Another MD study of laser ablation was performed by Herrmann et al.,campbell98  who examined the interaction of femtosecond laser pulses with silicon; a truncated version of the SW potential was used to account for the rupture of interatomic bonds following the absorption of photons. The predicted values of the threshold fluence for ablation thus obtained were at least an order of magnitude higher than those observed experimentally.cavalleri99 ; cavalleri98  This discrepancy is likely due to the inability of the SW potential to describe the non-thermal mechanisms involved in the interaction of femtosecond pulses with silicon. Indeed, it is generally agreedlinde97 ; silvestrelli96 ; linde00\_2  that there exists, for silicon, a critical carrier density $`n_c10^{22}\text{cm}^3`$ separating two distinct regimes:

(i) With picosecond and nanosecond pulses, carrier densities remain below the critical value and the relaxation is thermal, i.e., the carriers relax by transferring their kinetic energy to the lattice by the spontaneous emission of optical phonons in a characteristic time of $`1\text{ps}`$. In this context, structural modifications such as melting and ablation are observed on a $`1001000\text{ps}`$ time scale. This model, known as TAM, is suitable for laser pulses down to approximately 10 ps.wood84

(ii) With femtosecond pulses, on the other hand, the number of free carriers may exceed $`n_c`$ and ultrafast melting proceeds within $`1\text{ps}`$ by the collapse of the lattice.stampfli94 ; silvestrelli97  Because melting occurs before the carriers have time to relax through the emission of optical phonons, the relaxation is non-thermal. This model, first proposed by Van Vechten et al.,vanvechten79  is known as the plasma annealing model (PAM).murakami84

Furthermore, if the light intensity at the surface exceeds $`10^4\text{GW/cm}^2`$, the generation of a high number of energetic carriers through a two-photon absorption mechanism leads to the complete ionization of the material by avalanche ionization and/or field ionization: a very dense ($`n10^{23}\text{cm}^3`$) and hot ($`10^6\text{K}`$) $`e`$-$`h`$ gas is generated near the surface.schuler96 ; private  One then speaks of a solid-to-plasma transition as the material is carried away by the expanding plasma.

Here, because we make use of the SW potential, it is assumed that thermal processes only are operative; this restricts us to 10 ps or longer pulses. Care is taken to avoid carrier densities in excess of $`n_c`$; in this context, the intensity remains well below $`10^4\text{GW/cm}^2`$.

### II.2 Molecular dynamics

The MD technique has been discussed extensively in the literature (see, for instance, Ref. smith97, ). Briefly, the set of $`3N`$ Newton equations of motion, where $`N`$ is the number of atoms in the system, are solved from the knowledge of the interatomic forces. If $`\text{r}_i`$, $`\text{v}_i`$, and $`\text{a}_i`$ are the position, velocity and acceleration of the $`i`$th atom, respectively, the time evolution of the system is obtained by integrating, at discrete time steps, the following differential equations:

$`{\displaystyle \frac{\text{r}_i}{t}}`$ $`=`$ $`\text{v}_i,`$
$`{\displaystyle \frac{\text{v}_i}{t}}`$ $`=`$ $`\text{a}_i={\displaystyle \frac{1}{m_i}}\left[{\displaystyle \frac{\text{U}(\text{r}_1,\mathrm{},\text{r}_N)}{\text{r}}}\right]_{\text{r}=\text{r}_i},`$ (1)
given a set of initial positions and velocities {$`\text{r}_i^\text{0}`$,$`\text{v}_i^\text{0}`$}. In the above equations, $`m_i`$ is the mass of the $`i`$th atom and $`U`$ the interatomic potential. Thermodynamic quantities such as temperature and pressure can then be computed in a straightforward manner.allen87

All calculations were performed using the program groF, a general-purpose MD code for bulk and surfaces developed by one of the authors (LJL). The original form of the SW potential was used to compute the forces between the silicon atoms.stillinger85  The equations of motion were solved using the velocity Verlet algorithm.allen87  The MD time step was set to $`\mathrm{\Delta }t=0.5\text{fs}`$ for all simulations.

### II.3 Target and laser pulse

Simulations are carried out for a supercell with approximate dimensions $`3\times 3\times 60\text{nm}^3`$ (containing a total of 31680 atoms) representing a small volume of a Si(100) target located at the center of the laser pulse.

In order to recreate the thermal and structural constraints of a macroscopic crystal, the supercell is first repeated in the lateral ($`x`$ and $`y`$) directions using periodic boundary conditions (PBC). The thermal constraint is imposed by coupling a few monolayers of atoms at the bottom of the supercell to a heat reservoir; this is done by renormalizing the velocities to an appropriate Maxwell-Boltzmann (MB) distribution. These atoms also serve to minimize the reflection of pressure waves at the bottom of the supercell. This latter phenomenon affects the results by, in some cases, adding an artificial contribution to ablation. Though other techniques were proposed to attenuate such pressure waves,zhigilei99  the method described above gives satisfactory results. Additionally, a few monolayers of atoms attached to their equilibrium positions are placed beneath the heat reservoir in order to mimic a semi-infinite crystal. The crystal is assumed to be initially perfect.

The simulation of a 10 ps laser pulse, Gaussian in time as well as in the $`x`$ and $`y`$ directions and of macroscopic width, is accomplished assuming that any spatial variation of the irradiance can be neglected. This is valid because the laser spot diameter is assumed to be much larger than the region under study. The irradiance is thus, in effect, spatially constant over the infinite $`x`$-$`y`$ plane by virtue of the PBC.

The relatively large number of photons contained in a pulse, typically a few tens of thousands, ensures a uniform spatial distribution of the energy at the surface. The temporal Gaussian distribution is simulated by the successive arrival of planes of photons spanning the entire surface of the supercell, the number of which being determined by the instantaneous irradiance; they are separated in time by intervals ranging from $`\mathrm{\Delta }t`$ to typically $`10\times \mathrm{\Delta }t`$.

### II.4 Absorption of light

The absorption of light in silicon proceeds by the excitation of valence electrons via intraband and interband transitions, the frequency of which depends on several parameters, such as the density of free carriers and the wavelength and intensity of the laser pulse. The absorption of photons of energy $`h\nu `$ larger than the bandgap $`E_G`$ can induce three types of optical transitions: (i) An interband transition following the absorption of a single photon by a valence-band electron; a bond is broken and the electron is promoted to the conduction band, leaving a hole behind. The newly created pair shares an energy $`(h\nu E_G)`$, each carrier receiving an initial kinetic energy determined by the set of selection rules. (ii) A nonlinear, two- (or more) photon, interband transition, occurring when $`n2`$ photons are simultaneously absorbed. If $`nh\nu \varphi `$, where $`\varphi =4.85\text{eV}`$ is the silicon work function,handbook  the electron can be photoemitted; the condition $`nh\nu <\varphi `$, i.e., $`\lambda 255\text{nm}`$ when $`n=1`$, thus ensures that no photoemission takes place. (iii) An intraband transition, or free-carrier absorption, according to which the photon is absorbed by an electron already in the conduction band.

Following their creation, the carriers, because of momentum and energy conservation rules, occupy thin energy shells and a small volume of k space. In a characteristic time $`\tau _{eh}`$ of typically a few tens of fs, the hole and electron subsystems achieve, through electron-hole scattering, a quasi-equilibrium state described by a Fermi-Dirac (FD) distribution at a common electronic temperature $`T_e`$ which is higher than the lattice temperature $`T`$. Thus, after $`100\text{fs}`$, the carriers and the lattice constitute two decoupled subsystems, most of the energy remaining stored within the former ($`T_e>T`$).

If $`n`$ represents the conduction band carrier density, the balance equation for $`e`$-$`h`$ pairs can be written asvandriel87

$$\frac{n}{t}+\text{J}=G+R,$$
(2)
where J is the current density, $`G`$ the carrier generation rate resulting from the absorption of photons, and $`R`$ the net recombination rate given by

$$R=\gamma n^3+\delta (T_e)n,$$
(3)
with $`\gamma `$ the Auger recombination coefficient and $`\delta (T_e)`$ the impact ionization coefficient. The carrier generation rate $`G`$ is the sum of one- and two-photon processes, $`G=G_1+G_2`$, and can be written as

$$G=\frac{(1\mathrm{\Gamma })\alpha I(z,t)}{h\nu }+\frac{(1\mathrm{\Gamma })^2\beta I^2(z,t)}{2h\nu },$$
(4)
$`(1\mathrm{\Gamma })I(z,t)`$ being the intensity of the laser at depth $`z`$ below the surface, $`\mathrm{\Gamma }`$ the reflection coefficient at the surface, and $`\alpha `$ and $`\beta `$ the one-photon and two-photon interband absorption coefficients, respectively. The absorption proceeds according to the Beer-Lambert law

$$\frac{I}{z}=(\alpha +\mathrm{\Theta }n)I\beta I^2,$$
(5)
where $`\mathrm{\Theta }`$ represents the free-carrier absorption cross-section.

It is instructive to compare the two-photon to one-photon carrier generation rates, that is the ratio

$$\frac{G_2}{G_1}=(1\mathrm{\Gamma })I(z,t)(\frac{\beta }{2\alpha }).$$
(6)
If we choose $`\lambda =308\text{nm}`$, each photon carries an energy $`h\nu =4.03\text{eV}`$, larger than the bandgap of silicon at 300 K, $`E_G1.12\text{eV}`$.sze85  For Si at 300 K, $`\mathrm{\Gamma }=\mathrm{\Gamma }_\text{0}0.59`$,wood84  $`\alpha =\alpha _\text{0}1.5\times 10^6\text{cm}^1`$,wood84  and $`\beta 40\text{cm/GW}`$.murayama94  For an irradiance as high as $`10^3\text{GW/cm}^2`$, a value at least an order of magnitude higher than the ones considered in this work, Eq. (6) gives $`G_2/G_10.01`$. Moreover, it has been shownlietoila82  that free-carrier absorption is negligible at short wavelengths, typically below 1000 nm. Therefore, the one-photon interband transition is the dominant mechanism for absorption in the present case and all other processes are ignored; as noted above, this allows the photoemission of electrons to be neglected.

In practice, upon arrival at the surface, the $`x`$-$`y`$ position of a photon is determined at random following a uniform probability distribution and reflected, again at random, according todeunamuno89

$$\mathrm{\Gamma }(T)=\mathrm{\Gamma }_\text{0}+\mathrm{\Gamma }_R\times T,$$
(7)
where $`\mathrm{\Gamma }_R=4\times 10^5`$. Upon melting of the surface (see Sec. III.1), the reflectivity is assumed to be that of l-Si, i.e., $`\mathrm{\Gamma }_l=0.73`$.wood84  If not reflected, the photon is absorbed in agreement with Eq. (5), neglecting two-photon and free-carrier absorption, at depthallen87

$$z=\alpha ^1\mathrm{ln}(\xi ),$$
(8)
where $`\xi `$ is a number generated from a uniform distribution; $`\alpha `$ is the absorption coefficient, given byjellison82

$$\alpha (T)=\alpha _\text{0}\mathrm{exp}\left(\frac{T\mathrm{\Theta }_\text{1}}{\mathrm{\Theta }_\text{2}}\right),$$
(9)
where $`\mathrm{\Theta }_\text{1}=300\text{K}`$ and $`\mathrm{\Theta }_\text{2}=4680\text{K}`$.

The silicon atom closest to the absorption site, and which possesses at least one valence electron, is then excited, i.e., its number of valence electrons is decremented by one. An $`e`$-$`h`$ pair is created with an initial position given by the $`(x,y,z)`$ coordinates of the excited atom. The relatively large absorption coefficient at $`\lambda =308\text{nm}`$, combined with a depth of 60 nm for the simulation box, ensures that more than 99% of the energy is absorbed within the supercell.

The determination of the exact initial kinetic energy of each carrier requires the knowledge of the allowed optical transitions. However, because very fast carrier-carrier scattering leads to a FD distribution at temperature $`T_e`$ within a few tens of fs, it is reasonable to assume the gas to be in a state of an instantaneous quasi-equilibrium. Further, because $`T_e`$ is relatively large, typically a few thousand degrees, the FD distribution is well approximated by a MB distribution. Thus, as an electron is promoted to the conduction band, its initial kinetic energy is calculated according to a MB distribution at the current temperature $`T_e`$. The problem of determining the set of optical transitions as well as the carrier density and temperature-dependent chemical potential, required to characterize a FD distribution, can thus be avoided. Using a Box-Muller transformation,numrecipes  a MB distribution at $`T_e`$ is obtained by calculating the electron initial kinetic energy, $`E_{k_i}`$, according to

$$E_{k_i}=k_BT_e\underset{i=1}{\overset{3}{}}\mathrm{ln}(\zeta _i)\mathrm{cos}^2(2\pi \xi _i),$$
(10)
where $`\zeta _i`$ and $`\xi _i`$ are random numbers taken from a uniform distribution. The hole thus receives an initial kinetic energy equal to ($`h\nu E_GE_{k_i}`$). The relationship used for $`E_G`$ (in eV) as a function of the lattice temperature $`T`$ and the carrier concentration $`n`$ isvandriel87 ; pierret87

$$E_G=E_G^\text{0}\frac{\kappa T^2}{(T+\mathrm{\Lambda })}\varsigma n^{\frac{1}{3}},$$
(11)
where $`E_G^\text{0}=1.17\text{eV}`$, $`\kappa =4.730\times 10^4\text{eV/K}`$, $`\mathrm{\Lambda }=636\text{K}`$ and $`\varsigma =1.5\times 10^8\text{eV cm}`$.

### II.5 Relaxation processes

A number of microscopic mechanisms allow the system to restore equilibrium following irradiation by a laser pulse; these are carrier-phonon scattering, carrier diffusion, impact ionization and Auger recombination.

Impact ionization is the process through which an $`e`$-$`h`$ pair is created when an energetic electron collides with a valence electron: the former transfers part of its kinetic energy to the latter which is promoted to the conduction band. Because impact ionization most likely originates from the generation of highly-energetic carriers through two-photon interband transitions (negligible here; see Sec. II.4),schuler96  it can be ignored.

In silicon, the dominant recombination mechanism is Augerlandsberg91  for which the characteristic recombination time, due to screening effects at high carrier densities, is $`6\text{ps}`$.yoffa80  Because the carrier density gradient, which develops at the surface following laser irradiation, is taken into account, in our model, by assuming a unidirectional flow of electrons and holes along the $`z`$ axis into the bulk (see Sec. II.5.2), the latter have usually diffused out of the simulation box in a time less than 6 ps. As a result, Auger recombination can also be neglected.

The important relaxation mechanisms, in the present context, are carrier-phonon scattering and carrier diffusion, which we explicitly take into account.

#### II.5.1 Carrier-phonon scattering

If the carrier density $`n`$ remains below the critical valuestampfli94 ; silvestrelli96  $`n_c10^{22}\text{cm}^3`$, relaxation is thermal,linde97  i.e., the carriers relax by transferring their kinetic energy to the lattice by the spontaneous emission of (mainly) optical phonons in a characteristic time $`\tau _{LO}\text{1 ps}`$ (also known as the carrier-phonon energy relaxation time $`\tau _E`$). The population of optical phonons thus increases for approximately 1 ps, after which the latter relax by emitting acoustic phonons in a characteristic time $`\tau _{TA}\text{10 ps}`$. Eventually, the energy is redistributed among all the vibrational modes of the lattice and a state of quasi-equilibrium is reached, characterized by a Bose-Einstein distribution at temperature $`T`$. The further evolution of the system can then be described through thermal processes only as $`TT_e`$, but $`T`$ is still higher than its initial value.

There are several types of carrier-phonon scattering processes. It is possible to identify the most important ones, as well as their associated mean scattering rates. The dominant mechanisms for energetic electrons are the $`f`$ and $`g`$ intervalley scattering transitions, which involve large changes in momentum and thus acoustic or optical phonons with wave vectors near the zone boundary.lundstrom90

The probability $`P_{if}`$ of observing a given intervalley transition, according to which a carrier in the valley $`i`$, with crystal momentum pi, scatters to a state with crystal momentum pf in the valley $`f`$, is proportional to the square of the intervalley deformation potential, $`D_{if}`$, which characterizes the strength of the scattering.lundstrom90  The probability can then be estimated as the weight associated with $`D_{if}^2`$. The probabilities thus obtained are listed in Table 1 and were calculated using the corresponding values for $`D_{if}`$.fischetti88  As one can see, nearly 70% of phonons that scatter with electrons are optical. For holes, however, the main mechanism is optical deformation potential scatteringlundstrom90  and we thus assume the holes to scatter only with optical phonons of 62.6 meV.lundstrom90  Moreover, the energy of the carrier determines whether scattering involves the absorption or the emission of a phonon. In Si, phonon emission is largely dominant for carrier energies exceeding $`E_s\text{50 meV}`$, while phonon absorption dominates below $`E_s`$.lundstrom90

To determine the rate at which the energy is transferred from the electronic to the ionic degrees of freedom, the mean carrier-phonon scattering rate as a function of carrier energy, $`1/\tau _{cp}^\text{0}`$, is also needed. This quantity was calculated in Si at room temperature using the full band-structure of the material by Fischetti et al.fischetti88  A fit to their data was embedded in our computer model; variations with temperature were ignored.

Because of screening effects, however, the observed carrier-phonon scattering rate decreases for carrier densities $`n10^{21}\text{cm}^3`$, as recently observed by Sjodin et al.;sjodin98  their results confirmed theoretical predictions by Yoffa.yoffa81  Screening effects were thus implemented in our model as suggested by Yoffa, with the new carrier-phonon scattering rate given by

$$\tau _{cp}^1(n)=\frac{1/\tau _{cp}^\text{0}}{(1+(\frac{n}{n_{sc}})^2)},$$
(12)
where $`n_{sc}1.2\times 10^{21}\text{cm}^3`$ is the critical carrier density to observe significant screening effects; $`\tau _{cp}^\text{0}`$ is the value of $`\tau _{cp}`$ without screening. Finally, the probability of an electron suffering a collision with a phonon was assumed to be Drude-like: the probability, $`P_s`$, that a carrier has scattered with a phonon after a time $`t`$, where $`t`$ is the elapsed time since the last scattering event, is thus

$$P_s(t)=1\mathrm{exp}(t/\tau _{cp}).$$
(13)
In practice, we proceed as follows: after being generated, a carrier is given an initial position and kinetic energy as described in Sec. II.4 and its motion proceeds until a scattering event takes place. At each time step, the time elapsed since the last scattering event is determined, $`\tau _{cp}`$ is computed and the probability $`P_s`$ is determined from Eq. (13). If a scattering event occurs, a phonon is either emitted or absorbed depending if the carrier kinetic energy is greater or lower than $`E_s`$, respectively.

For an electron scattering with a phonon, the intervalley process is determined at random according to the probabilities of Table 1. A quantum of energy equal to $`\mathrm{}\omega _\text{0}`$ is then given to (phonon emission) or removed from (absorption) the lattice according to a spatial Gaussian distribution in a radius of 5 Å from the carrier, the latter value being an estimate of the carrier wave packet width based on the uncertainty principle. Though a finite duration is associated with each scattering event,bordone96  the latter are assumed to be instantaneous for simplicity. Finally, if a phonon scattering event has occurred, the carrier energy is updated.

#### II.5.2 Carrier diffusion

Carriers diffuse into the bulk as a result of the density gradient, $`n`$, which is itself a consequence of the exponentially-decreasing absorption with depth. The variations with $`n`$ of the ambipolar diffusion coefficient $`D`$ and carrier mobilities, which characterize the dynamics of carriers in semiconductors, have been the subject of controversy; the available experimental data at densities exceeding $`10^{19}\text{cm}^3`$ are very scarce and concerns mostly germanium.auston74 ; jamison76 ; moss81  In fact, for $`n>10^{19}\text{cm}^3`$, it is not clear whether the ambipolar diffusion coefficient decreases due to additional carrier-carrier scattering, as proposed by Fletcher,fletcher57  or increases due to many-body quantum effects.vandriel82  At the same time, various models have been proposedvandriel82 ; kane92 ; velmre99  to explain the influence of carrier-carrier scattering, as well as quantum effects, on $`D`$.

Recent experimental datalinnros94 ; rosling94  for moderate ($`10^{15}10^{17}\text{cm}^3`$) carrier densities tend to support Fletcher’s model. In view of this, and the lack of experimental evidence on many-body quantum effects for $`n>10^{19}\text{cm}^3`$, we adopt Fletcher’s model for $`D`$ implemented in our code using the expression suggested by Berz et al.berz79

Thus, each $`e`$-$`h`$ pair is given an ambipolar diffusion coefficient which depends upon the local carrier density and local temperature, and which is updated at every time step, during which a carrier travels a distance $`\mathrm{\Delta }z=\sqrt{D\mathrm{\Delta }t}`$ away from the surface.

## III Results and discussion

All simulations were carried out assuming a laser pulse duration $`\tau _L=10\text{ps}`$ and photon energy $`h\nu =4.03\text{eV}`$ ($`\lambda =308\text{nm}`$). The pulse intensity has a temporal Gaussian profile, but is spatially uniform. The fluence, $`F`$, was varied between 0.01 and 0.75 J/cm2; corresponding intensities are in the range 1–75 GW/cm2.

In Sec. III.1, we present our results for fluences of 0.01 to 0.20 J/cm2. No matter removal is observed; a breathing mode of the surface reveals the propagation of pressure waves in the solid. At and above $`F_{th}=0.25\text{J/cm}^2`$, large pieces of molten material are expelled from the surface; the mechanisms by which they are ejected are examined in Sec. III.2.

### III.1 Below the threshold energy for ablation

Fig. 2 shows the carrier temperature $`T_e`$ during the first two picoseconds for a fluence $`F=0.05\text{J/cm}^2`$; $`T_e`$ is obtained by summing the kinetic energies of the electrons and holes in the conduction and valence bands, respectively. The time evolution of the carrier density $`n`$ and the surface, $`T_s`$, and bulk, $`T_B`$, values of the lattice temperature are shown in Fig. 3 for the same fluence; here, $`T_s`$ is obtained by averaging over the first $`\delta =1/\alpha 7\text{nm}`$ below the surface, where $`\delta `$ is the optical skin depth, while $`T_B`$ represents an average over the rest of the supercell. The pulse starts at $`t=0`$ and the initial temperature of the target is $`300\text{K}`$. Upon arrival of the first photons, $`T_e`$ rises rapidly to $`11000\text{K}`$ but $`n`$ remains low: the pulse provides sufficient energy to create highly-energetic e-h pairs while not significantly increasing the number of carriers. The newly-generated hot electrons and holes promptly begin to thermalize with the lattice: within $`500\text{fs}`$ to 1 ps, $`T_e`$ drops to $`1500\text{K}`$; this indicates a very fast initial cooling rate of the carriers through phonon scattering. For delays $`1\text{ps}`$, $`T_e`$ decreases at a slower rate, i.e., the rate at which energy is being transferred from the electronic to the ionic degrees of freedom diminishes significantly. The carrier-phonon energy relaxation time, $`\tau _E`$, is thus estimated to be a few hundred fs. These observations are consistent with recent probing of ultrafast carrier dynamics in silicon by Goldman et al.goldman94  who reported a value of $`\tau _E1\text{ps}`$. Equilibrium with the lattice is finally achieved after $`10\text{ps}`$ with $`T_e\text{300 K}`$.

The carrier density $`n`$ in Fig. 3a reaches a peak at $`t7\text{ps}`$, thus well after $`T_e`$ has passed its maximum; this is because the e-h plasma loses energy to the lattice faster than it gains energy from the pulse. The decrease of $`n`$ for $`t>7\text{ps}`$ is due to diffusion of the e-h pairs away from the surface; moreover, the carrier density remains well below the critical value, $`n_c10^{22}\text{cm}^3`$, ensuring that the TAM is appropriate.

The pulse generates a high number of electrons and holes which, in turn, heat up the lattice; this is seen in Fig. 3b. For $`t>\tau _L`$, $`T_B`$ reaches a plateau at $`1250\text{K}`$; this value is below the melting temperature of c-Si, $`T_m=1685\text{K}`$.handbook  The surface temperature $`T_s`$, on the other hand, continues to increase beyond $`\tau _L`$, approaching $`T_B`$ asymptotically. Most importantly, $`T_s<T_B`$ for the whole duration of the simulation. This clearly indicates that the maximum lattice temperature is located beneath the surface, a direct consequence of carriers diffusing into the bulk as shown below.

The fact that the maximum temperature is found below the surface can be better appreciated from Fig. 4, where the lattice temperature $`T`$ is plotted as a function of depth $`z`$ at $`t=50\text{ps}`$ for a fluence $`F=0.20\text{J/cm}^2`$; $`z=0`$ corresponds to the surface. The temperature increases from a minimum of $`2000\text{K}`$ at the surface to a maximum of $`4000\text{K}`$ at $`z27\text{nm}`$; the latter is above the boiling point $`T_b=2753\text{K}`$, but below the critical temperature, $`T_c=5193\text{K}`$, of silicon.handbook  For $`z27\text{nm}`$, $`T`$ decreases. Thus, the maximum temperature is not occurring at the surface but, instead, at $`z30\text{nm}`$. These observations can be explained as follows: the pulse intensity decreases exponentially with depth and, consequently, more e-h pairs are created near the surface (see Sec. II.4). However, the time interval between two consecutive carrier-phonon scattering events being typically $`10\text{fs}`$ to a few ps, a carrier will diffuse, on average, a few nm away from the surface between two scattering events. Thus, a large proportion of the energy is not released at the surface but, rather, deeper in the bulk.

Clearly, the temperature in Fig. 4 does not decay exponentionally with depth, as one might have expected; rather, a steep temperature gradient has developed at $`z30\text{nm}`$. This is a clear indication of SSSH effects which are usually assumed to originate from a vaporization process at the surface.singh90  However, it is important to note that, in our simulations, no vaporization is observed; this contrasts, for example, with the results of Zhigilei et al. who identified evaporation as the mechanism of matter removal below the threshold energy for ablation.zhigilei97 ; zhigilei00  The absence of evaporation in the present simulations is not surprising if one computes the flux (in $`\text{atoms}\text{m}^2\text{s}^1`$) of atoms that desorb from the surface, $`\mathrm{\Phi }_e`$, given byohring

$$\mathrm{\Phi }_e=\frac{P_e(T_s^e)}{\sqrt{2\pi Mk_B}T},$$
(14)
where $`M`$ is the mass of a desorbed atom and $`P_e`$ is the pressure of the vapor in equilibrium with the solid at a temperature $`T_s^e`$. For $`F0.15\text{J/cm}^2`$, $`T_s2000\text{K}`$ (see Fig. 5), $`P_e2.11\text{Pa}`$,ohring  and we find $`\mathrm{\Phi }_e2.3\times 10^{22}\text{atoms}\text{m}^2\text{s}^1`$. The surface of the supercell has an area $`3\times 3\text{nm}^2=9\times 10^{18}\text{m}^2`$ and it is thus estimated that an atom desorbs every $`5\mu \text{s}`$, which is much longer than the duration of our simulations. Thus, the temperature profile in Fig. 4 cannot be explained by an evaporation process; rather, it is the result of carriers diffusing into the bulk.

Fig. 5 shows the surface temperature $`T_s`$ as a function of fluence. The threshold for melting, $`F_m`$, is thus estimated to be in the range 0.05–0.1 GW/cm2 since we must have $`T_sT_m`$ to observe melting. A more precise assessment of $`F_m`$ is obtained by computing the coordination number averaged over the first few atomic layers below the surface, $`\rho _\text{0}`$, as a function of time for various fluences. The results for $`F=0.09\text{J/cm}^2`$ are depicted in Fig. 6. A sudden change indicating the onset of melting is observed at $`t75\text{ps}`$: from $`3.88`$ (a value slightly below that of bulk c-Si, a consequence of including in the calculation undercoordinated surface atoms), $`\rho _\text{0}`$ approaches the coordination number of $`l`$-Si, i.e., $`6`$.wood84  Because no significant changes in the surface coordination number are observed at lower fluences, the predicted threshold is thus $`F_m0.09\text{J/cm}^2`$, melting taking place on a characteristic time of a hundred ps, in good agreement with experimental data.cavalleri99 ; bloembergen85 ; linde97

The mechanisms by which the system restores internal equilibrium below $`F_{th}`$ are revealed in Fig. 7 where the surface position, i.e. the $`z`$ coordinate of the topmost atomic layer, $`z_{top}`$, is plotted with respect to time for $`F=0.10\text{J/cm}^2`$; $`z_{top}=0`$ corresponds to the initial surface position. Also plotted is the bulk pressure $`P_B`$, averaged over the entire supercell. As one can see, oscillations with an amplitude of $`1\text{nm}`$ and a period $`30\text{ps}`$ are observed at the surface. This breathing mode is due to the propagation of pressure waves in the solid; the latter represent the elastic response of the material to the intense heating by the pulse. Upon inspection of Fig. 7, $`P_B`$ and $`z_{top}`$ are clearly out of phase: because silicon expands when heated, $`z_{top}`$ initially increases. At $`t20\text{ps}`$, $`z_{top}`$ is maximum and $`P_B`$ is minimum ($`P_B<0`$, tension); this is indicative of tensile stresses. In order to restore internal equilibrium, the material subsequently enters a compressive phase and contracts, i.e., $`z_{top}`$ decreases. At $`t35\text{ps}`$, $`z_{top}`$ is minimum and $`P_B`$ is maximum ($`P_B>0`$, compression); the solid is then under compressive stresses. It is also apparent from Fig. 7 that melting occurs at the surface: the average value of $`z_{top}`$ decreases with time, an indication that the material is more compact. Indeed, molten silicon is denser than crystalline silicon.handbook

The surface velocity, $`v_{top}`$, and acceleration, $`a_{top}`$, computed from the data in Fig. 7, are found to be $`100\text{m/s}`$ and $`10^{13}\text{m/s}^2`$, respectively. These values are confirmed by noting that, upon heating, a surface undergoes an expansion $`\mathrm{\Delta }z`$ given bytam91

$$\mathrm{\Delta }z=(1\mathrm{\Gamma })F\alpha _T/\rho C,$$
(15)
where $`\mathrm{\Gamma }`$ is the reflectivity, $`\alpha _T`$ is the linear thermal expansion coefficient, $`\rho `$ is the density and $`C`$ is the specific heat. Assuming $`\mathrm{\Gamma }0.59`$, $`\alpha _T2.6\times 10^6\text{K}^1`$,silicon  $`\rho 2.33\text{g/cm}^3`$, $`C0.713\text{J/gK}^1`$,silicon  and a fluence $`F0.10\text{J/cm}^2`$ at 300 K, Eq. (15) yields an oscillation amplitude $`\mathrm{\Delta }z0.6\text{nm}`$, a value consistent with the results of Fig. 7. With a pulse duration $`\tau _L=10\text{ps}`$, the velocity and acceleration of the surface become $`v_{top}\mathrm{\Delta }z/\tau _L60\text{m/s}`$ and $`a_{top}\mathrm{\Delta }z/\tau _L^210^{13}\text{m/s}^2`$, thus in good agreement with the above results. In addition, Fig. 8 gives the maximum value of $`v_{top}`$ and $`a_{top}`$ as a function of laser fluence. A maximum is observed at $`F=0.15\text{J/cm}^2`$. The subsequent decrease at $`F=0.20\text{J/cm}^2`$ might result from a stronger damping of the oscillations due to a thicker molten layer at the surface. Laser-induced oscillations at surfaces have found several applications. In dry laser cleaning, for instance, the removal of contaminants, i.e., adsorbed micron-sized particulates, is driven by momentum transfer: though of small amplitude ($`1\text{nm}`$), the surface oscillations take place on a very short time scale and thus exhibit very large accelerations.

Finally, Fig. 9 shows the pressure $`P`$ as a function of depth $`z`$ at $`t=0`$ and $`t=10\text{ps}`$; the fluence is $`F=0.10\text{J/cm}^2`$. As expected, the pressure is zero throughout the material before laser irradiation. By the end of the laser pulse, a pressure gradient has built up at $`35z40`$ nm; we show in Sec. III.2 that the latter is responsible for matter ejection at fluences $`F0.25\text{J/cm}^2`$. Below the threshold energy for ablation, however, the pressure gradient is not high enough to cause ejection of matter, and the system relaxes through the generation of pressure waves.

### III.2 Above the threshold energy for ablation

The threshold energy for ablation, $`F_{th}`$, is defined in this work as the minimum fluence for the onset of matter removal. No evaporation is observed (see Sec. III.1); rather, a collective ejection process is obtained at all fluences above $`F_{th}`$. Here, $`F_{th}=0.25\text{J/cm}^2`$, corresponding to a laser intensity $`I_{th}=25\text{GW/cm}^2`$.

Fig. 10 shows the surface position, $`z_{top}`$, as a function of time for a fluence $`F=0.35\text{J/cm}^2`$. As is the case for $`F<F_{th}`$, the material initially expands upon heating. However this does not result, this time, in a surface breathing mode. Instead, a dramatic change occurs at $`t\text{30 ps}`$: from a value of $`7\text{nm}`$, $`z_{top}`$ drops to $`35\text{nm}`$ within a few ps. This indicates that ablation has taken place over a depth $`h35\text{nm}`$, a value consistent with experiment.cavalleri99  No significant change of $`h`$ is observed at higher fluences.

A snapshot of the ablation process at $`t=33\text{ps}`$ for the same fluence is given in Fig. 11. The photons were incident from the top. The atoms that have absorbed at least one photon are colored in dark gray, others are in light gray. Evidently, most photons have been absorbed near the surface.

Two features should be noted: (i) Ablation is a bulk phenomenon, i.e., the ejection of matter originates in the bulk and not at the surface as one might have expected. More specifically, ablation results in the expulsion of a large piece — a few tens of nm thick — of molten material from the surface with an initial axial velocity of $`1000\text{m/s}`$. (ii) As opposed to a disordered network in the bulk, a few nm thick crystalline layer is visible at the surface. The thickness of the latter corresponds to the average distance, about 5 nm, carriers created at the surface travel before suffering a first collision. For $`z5\text{nm}`$, the temperature rises rapidly and, as shown below, the thickness of the layer removed during ablation depends upon the location of the pressure build-up in the bulk.

Insight into the mechanisms responsible for ablation is found in Fig. 12 showing the pressure, $`P`$, and the lattice temperature, $`T`$, as a function of depth for $`F=0.35\text{J/cm}^2`$. A steep temperature gradient has developed at the surface by $`t=10\text{ps}`$, Fig. 12a: from a value of $`2000\text{K}`$ at $`z=0`$, $`T`$ rises rapidly to reach a plateau slightly above the critical temperature, $`T_c`$, for $`z30\text{nm}`$. As shown in Sec. III.1, the occurrence of a maximum temperature below the surface is not a consequence of evaporation from the surface, nor does it result from a thermal diffusion process carried by phonons. Rather, the energy is transported into the bulk by the electronic degrees of freedom; this mechanism is, of course, absent in organic solids.zhigilei00

Intense heating by the pulse leads to the thermal confinement of the laser-deposited energy and, as shown in Fig. 12a, the material is heated far beyond its boiling point by the end of the pulse. At the same time, the lattice cannot undergo sufficient thermal expansion and a strong pressure ($`>0`$) in excess of 10 GPa develops at $`35z50\text{nm}`$, Fig. 12b. By $`t24\text{ps}`$ the latter has relaxed, leaving important tensile stresses in the material; these are responsible for the rapid expansion observed for $`t10\text{ps}`$ in Fig. 10. The pressure profile at $`t=24\text{ps}`$ in Fig. 12b also reveals a local minimum (in absolute value) at $`z35\text{nm}`$. As a consequence, strong tensile forces of opposite directions drive the material apart locally and nucleation of voids occur, as depicted in Fig. 13. The latter begins to form at $`t=27\text{ps}`$ and ablation follows a few ps later as the mechanical strength of the material is exceeded.

As seen above, the high compressive stress which develops within the absorbing volume is responsible for the subsequent tensile forces; thus, it is the driving force responsible for matter removal above $`F_{th}`$. Fig. 14 shows the variation of the maximum compressive stress, $`P_{max}`$, with laser fluence. Below the threshold for ablation, $`P_{max}`$ increases with $`F`$. Above $`F_{th}`$, however, $`P_{max}`$ saturates.

In a recent study of organic solids,zhigilei00  Zhigilei et al. reported the observation of void nucleation in the stress confinement regime and attributed it to the build-up of a strong pressure gradient within the absorbing volume. The condition of stress confinement is expressed by $`\tau _L<\tau _S\delta /\nu _S`$, where $`\nu _S`$ is the velocity of sound in the material. For c-Si, $`\nu _S9\times 10^3\text{m/s}`$ and,handbook  for a laser pulse at $`\lambda =308\text{nm}`$, $`\tau _S1\text{ps}`$. This value is $`<\tau _L`$ by a factor of 10; conditions of stress confinement are thus unlikely to develop in silicon following irradiation with 10 ps or longer pulses. The regime of thermal confinement, on the other hand, is expressed by $`\tau _L<\tau _{th}\delta ^2/D_T`$, where $`D_T`$ is the thermal diffusivity.zhigilei00  For silicon at 300 K, $`D_T0.86\text{cm}^2`$/s;sze85  this corresponds to $`\tau _{th}0.5\text{ps}`$, a value which is also about an order of magnitude lower than $`\tau _L`$. However, $`D_T`$ drops to $`0.10\text{cm}^2`$/s as $`T`$ reaches $`1500\text{K}`$,silicon  and it follows that $`\tau _L\tau _{th}`$: this suggests that conditions of thermal confinement are likely to arise when silicon is heated beyond $`T1500\text{K}`$ in a typical time of 10 ps or less. This phenomenon is responsible for the rapid heating of the material up to $`TT_c`$ in Fig. 12a. However, one should note that the presence of a steep temperature gradient at the surface rules out the possibility of phase explosion mechanisms which require that $`T/z0`$.miotello95 ; kelly96  Consequently, ablation is here the result of SSSH effects brought about by the thermal confinement of the energy.

Thus, the present simulations suggest that picosecond laser ablation of silicon is driven by SSSH effects. However, a more precise assessment of the matter removal mechanisms should involve the implementation of surface effects; these have been observed in a study of hot carrier dynamics at a Si(100) surface.goldman94  Accounting for surface effects is currently under development.

## IV Concluding remarks

A molecular-dynamics thermal annealing model (MADTAM) has been applied to the study of laser ablation of silicon with picosecond pulses on a 100 ps time scale. A detailed description of the thermal annealing model on a microscopic scale has been embedded into a molecular-dynamics scheme; this was accomplished by explicitly accounting for carrier-phonon scattering and carrier diffusion.

The model predicts a melting threshold $`F_m=0.09\text{J/cm}^2`$ and an ablation threshold $`F_{th}=0.25\text{J/cm}^2`$. Below $`F_{th}`$, oscillations of the surface reveal the propagation of pressure waves in the solid. Above $`F_{th}`$, important subsurface superheating (SSSH) effects are responsible for matter removal: as a result of the thermal confinement of the laser-deposited energy, the material is overheated to a temperature corresponding to the critical point of silicon and a pressure gradient builds up a few tens of nm beneath the surface. Ejection of matter is initiated by the relaxation of the strong compressive stresses after a few tens of ps, large pieces of molten material being expelled from the surface with high axial velocities.

We are currently investigating the role of surface effects. These are believed to be responsible for the formation of a surface space charge layer which is expected to influence the carrier dynamics and to compete with carrier diffusion.

###### Acknowledgements.

PL wishes to thank L. V. Zhigilei for the stimulating discussions. PL is also grateful to Dietrich von der Linde, Henry M. van Driel as well as to Massimo V. Fischetti for their patience in answering questions. This work was supported by the Natural Sciences and Engineering Research Council of Canada and the “Fonds pour la formation de chercheurs et l’aide à la recherche” of the Province of Québec. We are indebted to the “Réseau québécois de calcul de haute performance” (RQCHP) for generous allocations of computer resources.