Squeezed Thermal Phonons Precurse Nonthermal Melting of Silicon as a Function of Fluence

Eeuwe S. Zijlstra,* Alan Kalitsov, Tobias Zier, and Martin E. Garcia Theoretical Physics, University of Kassel, Heinrich-Plett-Strasse 40, 34132 Kassel, Germany Center for Interdisciplinary Nanostructure Science and Technology (CINSaT), Heinrich-Plett-Strasse 40, 34132 Kassel, Germany Department of Physics, University of Puerto Rico, San Juan, Puerto Rico 00931, USA (Received 14 June 2012; revised manuscript received 18 October 2012; published 29 January 2013)


I. INTRODUCTION
The interaction of an intense femtosecond-laser pulse with a solid typically leads to a nonthermal state that usually exists for approximately 1-10 ps, where the electrons acquire a temperature of several 1000-10 000 K and the atoms move on the potential energy surface created by the hot electrons [1,2].Ultrashort laser pulses thus provide a direct means to transiently manipulate the electrons that are responsible for many material properties, including structural stability.As a result, depending on the laser fluence, the atoms may start to move along unconventional trajectories within ultrashort times.One of the possible ensuing effects is the melting of a crystal within a few hundred femtoseconds, which was discovered three decades ago and has since been observed in many semiconductors and semimetals, both experimentally [3][4][5][6][7][8] and in simulations [9].In particular, in silicon, gallium arsenide, and indium antimonide, theoretical studies have demonstrated that nonthermal melting, i.e., melting out of thermodynamic equilibrium due to laser-induced changes in the potential energy surface, is initiated through a laserinduced lattice instability of transverse acoustic phonons at the Brillouin-zone boundary [10,11].
In this article, we direct our attention to excitations that are not sufficiently intense to melt a material and we investigate which laser-induced structural modifications can be used to predict the melting transition by looking for collective atomic motions below the melting threshold.Besides being of fundamental interest, we expect that knowledge of this process should also have implications in the field of femtosecond-laser materials processing.In previous attempts to find a precursor to ultrafast melting, it has been proposed that, for fluences below the melting threshold, the softening of coherent optical À-point phonons provides a relevant monitoring parameter [12].However, their different locations in the Brillouin zone indicate that the atoms follow different pathways during ultrafast melting and À-point phonon softening, which are therefore not directly related.So, roughly 30 years after its discovery, it remains unclear which physical process precurses nonthermal melting.

II. ATOMIC PATHWAYS
In order to address this open question, we have performed ab initio molecular-dynamics simulations of silicon, in which we explicitly take into account the hot-electron plasma created by femtosecond-laser excitation for different fluences (see the Appendix for details about the method).We study supercells containing ðn þ 1Þ Â n Â n conventional unit cells with N ¼ 96, 288, 640, and 1200 atoms in total.Atomic velocities and displacements are initialized so as to reproduce a Maxwell-Boltzmann distribution with a temperature of 1 mhartree (316 K).In Fig. 1(a), we show the time dependence of the root-mean-square atomic displacements from their equilibrium positions during the first 210 fs after femtosecond-laser excitation for absorbed fluences of 7, 12, 17, and 22 mJ=cm 2 , where we assume that nearultraviolet light with a penetration depth of 10 nm heats the electrons to 40, 50, 60, and 70 mhartree, respectively.Each curve shows the average over 60 runs for a 96-atom supercell.We see that, for the two highest electronic temperatures, the atoms move away rapidly from their lattice sites without bound, which is an indication that the silicon melts (see Video 1).Noticeably, both curves already exceed the Lindemann stability limit (equaling approximately 15% of the Si-Si nearest-neighbor distance, i.e., 0.35 A ˚) within 150 fs.In contrast, the atomic motions for the two lowest curves stay bound and we observe an oscillatory behavior that has not been noticed before in silicon: At the electronic temperature of 50 (40) mhartree, the root-mean-square atomic displacement reaches a maximum after 146 (98) fs and decreases afterwards.As we will elaborate below, this oscillatory behavior is a direct manifestation of the thermal squeezing of the classical ensemble of acoustic phonons in silicon, where the mean-square atomic displacements and momenta rotatively dip below their thermal averages in the laser-excited potential.In Fig. 1(b), snapshots after 146 fs of a particular run for each of the four selected electronic temperatures are superimposed.Strikingly, the atoms move in essentially the same directions for all electronic excitation densities, demonstrating that the same atomic pathways are followed during the first stages of nonthermal melting at high fluences and of phonon squeezing at lesser fluences.

III. MICROSCOPIC MECHANISM
A squeezed state of a classical oscillator is characterized by a variable that deviates less from its average value than at thermal equilibrium, while its conjugated variable has a larger variance.This so-called thermal squeezing has first been experimentally observed in a pumped microcantilever, where the cosine component of the time-averaged shot noise could be deamplified while the sine component was simultaneously amplified [13], and has been exploited for subthermal noise measurements [14,15].Whereas in Refs.[13][14][15] the squeezing of an oscillator has been achieved by means of a modulation of the macroscopic spring constant, in solids, it is possible to manipulate microscopic bonds by a femtosecond-laser pulse [16,17] that leads to squeezing [18].The microscopic mechanism of femtosecond-laser-induced thermal squeezing is schematically illustrated in Figs.2(a)-2(e) for an ensemble of nearly degenerate phonon modes.Before laser excitation, the thermal distribution of the atomic displacements is in equilibrium with the harmonic potential [Fig.2(a)].Through a femtosecond-laser pulse, the potential softens almost instantaneously, rendering the initial distribution function too narrow (squeezed) for the actual potential, so that, on average, the atoms start to move outward [Fig.2(b)].After a quarter of a phonon period, the distribution reaches its maximum width, which is wider than the equilibrium distribution of the laser-excited potential [Fig.2(c)].Thereupon, the distribution narrows again and, after half a period, when all oscillators have approximately finished a 180 phase shift, the displacements are back near their initial absolute values, albeit with opposite signs.The narrowness of the initial distribution is, however, not fully regained due to anharmonicities, phononphonon interactions, and a temporally diminishing constructive interference between phonon modes with different frequencies [Fig.2(d)].Depending on the strengths of these effects, further oscillations may be observed [Fig.2(e)].From the general nature of this mechanism, it follows that phonon squeezing is a common phenomenon that must emerge after intense femtosecond-laser excitation below the melting fluence in materials that exhibit laser-induced bond softening, which is a prerequisite for nonthermal melting.In agreement with our theory, we note that thermal phonon squeezing in bismuth [18] has been observed in the same general direction, where at higher fluences a lattice instability induces nonthermal melting [8].

IV. THERMAL SQUEEZING IN REAL SPACE
Further insight into the nature of the atomic motions can be obtained from the time-dependent variance hu 2 i À hui 2 of the atomic displacements u in the harmonic approximation (see the Appendix), from which we notice that N ¼ 640 already gives a sufficient sampling of all phonon modes in the first Brillouin zone [Fig.2(f)].By comparing the classical and quantum-mechanical time-resolved variances [Fig.2(f)], we conclude that, even though the meansquare atomic displacements vary between only 1.5 and 4.4 times the zero-point variance, quantum effects are not important for thermal phonon squeezing in silicon.
Figure 2(g) shows the results of our molecular-dynamics simulation with 640 atoms per supercell.Comparison with Fig. 2(f) demonstrates that anharmonicities and phononphonon interactions, which are fully included in the molecular-dynamics simulations and totally absent in the harmonic approximation, lead to an approximately 30% lengthening of the oscillation period and an equally large relative increase of the height of the first variance maximum.In Fig. 2(h), we draw phonon densities of states before and after laser excitation (see the Appendix).Three peaks, 1 ¼ 2:3, 2 ¼ 7:9, and 3 ¼ 12:6 THz, are discernible.Using the orthonormal eigenvectors fe i g of the 3N Â 3N dynamical matrix, we define the projection operators as e i e T i and P 2; The sum runs over all eigenpairs up to the first minimum in the phonon density of states, so that hðP 1 uÞ 2 i gives the mean-square atomic displacements projected onto the directions of the phonon modes with frequencies that approximately equal 1 and so that hðP 2;3 uÞ 2 i gives the contributions from all other lattice vibrations.In Fig. 2(g), we plot these spectrally projected variances along with the total hu 2 i, showing that acoustic phonon modes, which already swing with the largest amplitudes before laser excitation, dominate squeezing in real space.

V. THERMAL SQUEEZING IN MOMENTUM SPACE
A more even spectral weighting is obtained by considering the conjugated variables, i.e., the atomic momenta, whose distribution has a variance proportional to the instantaneous atomic temperature k B T atomic ¼ mhv 2 i=3 that we plot together with properly normalized spectral projections in Fig. 3.Although initially the phonons share a common temperature, the laser excitation is seen to lead to an average cooling of the atoms that is most pronounced for the phonon modes with frequencies around 1 .We further note that acoustic and optical phonons have not reached a common temperature after 600 fs, which is a new finding that suggests a line along which semiempirical theories, like the two-temperature model, can be extended.On top of these time-averaged effects, we notice at least two oscillations with periods that approximately equal 1=2 1 and 1=2 3 , demonstrating that phonon squeezing is not limited to the acoustic lattice vibrations that are responsible for nonthermal melting but that it also occurs in the optical part of the spectrum.In analogy, we infer that thermal squeezing is likely to occur in almost any material after femtosecond-laser excitation, as long as there is some modification of the interatomic force constants.

VI. FURTHER CONSIDERATIONS
In addition to the atomic forces that arise from laserexcited potential energy surfaces, which are fully incor-porated into the present simulations, incoherent electronphonon interactions will eventually lead to a thermalization of the electrons and the atoms.The essential step in this equilibration process is the emission of high-frequency phonons when hot carriers, i.e., electrons and holes, cool down [19].At very low excitation densities, emission of phonons happens on a time scale of 0 ¼ 200-260 fs (see [20]), but due to screening by the hot carriers this process has theoretically been predicted to slow down quadratically with the excitation density according to where is the excitation density, i.e., the number of electron-hole pairs per unit volume, and c is the so-called critical density [21].A recent experiment on silicon performed at ¼ 2:2 Â 10 À21 cm À3 has shown exponential atomic heating with a time constant ¼ 2 ps and has been explained by the above theory using 0 ¼ 230 fs and c ¼ 8 Â 10 20 cm À3 [20].In our simulations performed at the electronic temperature of 50 mhartree, ¼ 1:5 Â 10 22 cm À3 , so that we can expect the screening to become considerably more effective, but the estimate of 81 ps from the above theory is probably too large.A second interaction that is important for the electron-atom thermalization process is Auger recombination, where an electronhole pair recombines, transferring its energy to a third carrier and thereby keeping the electrons and holes hot [21].Competing with the above heat-exchange processes is the diffusion of hot carriers away from the excited surface region.In a thin-film geometry, where diffusion into the bulk is irrelevant, a complete thermalization of the excited carriers initially at 50 mhartree with the lattice initially at 316 K would lead to a phonon temperature that equals approximately 5700 K, assuming the high-temperature heat capacity of the lattice, which is well above the melting point of 1687 K.The estimated increase of the atomic temperature during the first 600 fs after the laser pulse depends on the exponential time constant , which is currently unknown for our excitation density but should, according to [20,21], lie between 40 K ( ¼ 81 ps) and 1408 K ( ¼ 2 ps).As mentioned above, in a bulk geometry, these values can be considerably lower due to carrier diffusion, which is an important process because of the very short penetration depth of near-ultraviolet light that equals approximately 10 nm.Additionally, by tuning the laser wavelength, it might be possible to excite electrons just above the band gap, likewise reducing the energy transfer from electrons to atoms.In any case, a rapid loss of coherence in the electronic subsystem in Si causes the electron-phonon thermalization to occur through incoherent scattering events.In bismuth, this has been shown to lead to a roughly linear increase of the variance of the atomic displacements with time, which is superimposed on the phonon-squeezing effect [18].In silicon, we expect a similar additional linear increase of the variance in Fig. 2(g) and of the ionic temperature in Fig. 3.However, in view of the above discussion, the slopes are hard to FIG. 3. Total and spectrally projected instantaneous atomic temperatures as a function of time after femtosecond-laser excitation.The data used and the decomposition are the same as in Fig. 2(g).
determine and will definitely depend on experimental parameters, such as the sample thickness and the wavelength of the ultrashort laser pulse.

VII. RELATION TO EXPERIMENTS
To complete our story, we now briefly discuss previous occurrences of phonon squeezing in crystals in the literature in the context of our new findings.First, quantum or vacuum squeezing has been induced by Garrett and coworkers in a transparent medium, KTaO 3 , using secondorder Raman scattering, which does not induce bond softening [22].The goal of this study is to squeeze the phonon wave packet below its zero-point width, which requires very low temperatures.Thermal squeezing of phonons, i.e., squeezing at elevated temperatures, has been observed in an opaque material, namely, bismuth [18].This study [18] has been limited to the low-fluence regime, well below the threshold for nonthermal melting, but in Sec.III we have linked thermal squeezing in Bi at low fluences and the direction in which it has been observed to nonthermal melting at higher fluences, in agreement with our theory.It is, however, important to point out that the precise role of the atomic displacements perpendicular to the eigenvector of the coherently excited A 1g mode during nonthermal melting has not yet been confirmed.In another study that is worth mentioning, harmonic models that are used to compare with experimental data for InSb that is excited at fluences both above and just below the melting threshold show oscillations of the root-mean-square atomic displacements [Eq.( 2) and Fig. 2 in Ref. [23] and Eq. ( 5) and Fig. 7 in Ref. [24] ], which are a direct manifestation of phonon squeezing comparable to our results in Fig. 2(g).Unfortunately, this effect could not be observed experimentally in the time-dependent Debye-Waller factors [23,24].This discrepancy between harmonic theory and experiment has been attributed to the models, in particular, to the absence of anharmonicities therein [23,24], which are indeed far from negligible, as can, for example, be seen either by comparing Figs.2(f) and 2(g) or from Refs.[10,11].It is, however, important to note that harmonic models cannot describe melting.Only a computation of the full pathways, such as the one performed here, can establish a connection between complex confined and unbound atomic motions.We therefore believe that the absence of squeezing in the experimental study of InSb [23,24] should not be explained as an artifact of the models but that it is probably due to other, experimental reasons, such as insufficient time resolution or large statistical errors.

VIII. PRECURSOR TO NONTHERMAL MELTING
In summary, we have found that phonons in silicon become thermally squeezed after relatively low-fluence femtosecond-laser excitation and we have shown that the atoms follow the same pathways in this process and during the first stages of nonthermal melting at higher fluences, demonstrating a close relationship at the microscopic level between these two seemingly independent phenomena.We have attributed these two processes to a common origin, namely, the laser-induced softening of the interatomic force constants, where, however, anharmonicities and phonon-phonon interactions enhance the squeezing effect by 30%.Based on a comparison of the atomic pathways, we have established that phonon squeezing is the precursor to nonthermal melting as a function of fluence.In this capacity, it has some predictive powers.First, if we expect no squeezing after femtosecond-laser excitation because the phonons neither soften nor harden considerably, laserinduced melting is a thermal process.This is, for example, the case in aluminum [16].The time scale of thermal melting is governed by incoherent electron-lattice equilibration, which lies typically in the range of picoseconds [25,26] but might also be faster, leading possibly to ultrafast thermal melting.Second, systems exhibiting nonthermal melting will definitively show thermal phonon squeezing at fluences below the melting threshold.Thermal squeezing may thus be used to experimentally obtain information about the nonthermal or thermal microscopic nature of ultrafast melting.It is also interesting to point out that, if there is laser-induced bond hardening, e.g., in gold [16], from our physical picture [Figs.2(a)-2(e)] we expect the squeezing oscillations to have the opposite sign compared to silicon, meaning that the distribution of displacements initially becomes narrower.For gold, delayed thermal melting has been reported [27].In the case of laser-induced bond softening, the squeezing is in the reported direction [cf.Figs.2(b) and 2(c)] and its amplitude increases with fluence up to the melting fluence.Thermal phonon squeezing can thus be used to follow and characterize the nature of ultrafast melting up to the melting transition as a function of fluence and could also provide a control parameter for femtosecond-laser-induced materials processing, which can relatively easily be measured as oscillations in the Bragg intensities.

ACKNOWLEDGMENTS
We thank BMBF (Project No. 05K10SJA) and DFG (Project No. GA 465/15-1) for funding.Computations were performed at the ITS, University of Kassel.

Ab initio molecular-dynamics simulations after femtosecond-laser excitation
To describe the effect of a hot-electron plasma, we used electronic-temperature-dependent density-functional theory [28], where the occupancies of the Kohn-Sham orbitals [29] are described by a Fermi function.This model, which is usually used to describe silicon after intense femtosecond-laser excitation, implicitly assumes very fast electron-hole equilibration.An interesting alternative model completely neglects electron-hole recombination [30].However, in the current work, we have restricted ourselves to the first more usual model.In order to be able to treat sufficiently large supercells, we have developed our own code for highly excited valence electron systems (CHIVES) [31], which uses norm-conserving pseudopotentials [32]; a primitive Gaussian basis set with exponents a 1 ¼ 1:177 69a À2 0 (s and p orbitals), a 2 ¼ 0:403 48a À2 0 (s, p, and d orbitals), and a 3 ¼ 0:129 89a À2 0 (s and p orbitals); and a regular three-dimensional grid with a cutoff energy of 1330 eV to describe the Hartree and the exchange and correlation potentials in the local density approximation [33].For our 96-atom supercell, we use a 1 Â 2 Â 2 k grid.All other cells are treated without noticeable loss of accuracy in the À-point approximation.
With these choices, we reproduce the laser-induced softening of the Brillouin-zone-boundary phonons responsible for nonthermal melting [10].In particular, for the transverse acoustic phonons at the X point, we find a groundstate frequency of 4.30 THz that decreases to 0 THz at an electronic temperature that equals approximately 65 mhartree, in excellent agreement with previous density-functional-theory calculations [34].Using the computed forces and the velocity Verlet scheme, we perform molecular-dynamics simulations with a time step of 2 fs, which is approximately 1=40 3 [Fig.2(h)].We initialize the atomic velocities and displacements before laser excitation using inverse transform sampling and 6N true random numbers [35] r i lying uniformly on [0,1] as

Phonons
We calculate all phonon modes that are compatible with a supercell by displacing an arbitrary atom from its equilibrium position by 0:001a 0 , first in the x direction and then in the y direction.Using the computed forces on all the atoms and the symmetry of the lattice, we construct the 3N Â 3N dynamical matrix, whose orthonormal eigenvectors e i are phonon directions that are related to the usual phonon vectors by " , where m is the atomic mass of silicon.We obtain phonon frequencies i from the eigenvalues i ¼ m! 2 i with !i ¼ 2 i .For the initial eigenvectors and eigenfrequencies at low electronic temperature before laser excitation, we use the notation e ðiÞ i and ðiÞ i .We denote the phonon directions and frequencies after laser excitation (electronic temperature of 50 mhartree) by e i and i .We compute phonon densities of states before and after laser excitation by, respectively, convolving the frequencies ðiÞ i and i , obtained for the N ¼ 1200 supercell, with Gaussians with full widths at half maximum of 1 THz [Fig.2(h)].

Harmonic theory of thermal phonon squeezing
In order to derive an expression for the time-dependent variance of the atomic displacements after laser excitation [Fig.2(f)], we first consider the stationary distributions of the atomic velocities and displacements in thermodynamic equilibrium at a temperature k B T ðiÞ atomic ¼ 1 mhartree, where k B is the Boltzmann constant, which applies to the situation before laser excitation.We call the vectors of length 3N with the velocities and displacements of all atoms v and u, respectively, and define the displacements and velocities in the phonon directions before laser excitation as which is our final result per atom for the classical timedependent variance of the atomic displacements in the harmonic approximation.A comparison with quantummechanical theory [18] shows that quantum effects modify the above equation by the substitution k B T ðiÞ atomic !@! ðiÞ The quantum-mechanical zero-point variance in the laserexcited potential equals

FIG. 1 .
FIG. 1. Atomic pathways induced by a femtosecond-laser pulse for four different fluences.(a) Averaged root-mean-square displacements from equilibrium as a function of time.The vertical linewidths are twice the standard deviations of the averages shown.The dashed line indicates the Lindemann stability limit.(b) Snapshot after 146 fs of a particular run for each excitation strength.The atomic coordinates are projected in the z direction.The lines indicate ð1; À1; 0Þ lattice planes.

FIG. 2 .
FIG. 2. Time-dependent variance of the atomic displacements at the electronic temperature of 50 mhartree.(a)-(e) Schematic illustration of thermal phonon squeezing.The black curves represent laser-induced changes in the potential.The blue curves represent the distribution of atomic displacements.(f) Ab initio results in the harmonic approximation.The blue curve is for N ¼ 640.The red curves are for N ¼ 1200.The solid curves represent classical results.The dashed curve represents the quantummechanical result.The hatched region represents the zero-point motion in the laser-excited state.(g) Ab initio moleculardynamics results: total and spectrally decomposed variances averaged over nine runs (N ¼ 640).The light shaded areas represent standard deviations of the averages.(h) Phonon density of states before (purple curve) and after (black curve and shaded areas) femtosecond-laser excitation in the range from 0 to 17.5 THz.
After the laser excitation, we obtain new phonon directions e i that are related to the initial directions by with coefficients C ij ¼ e T i Á e ðiÞ j that follow immediately from our computed phonon directions.The atomic displacement in the direction e i at time t ¼ 0 after laser excitation is given by u i ðt ¼ 0Þ e T i Á u ¼