Thermal modeling and numerical simulations have been performed to describe the ultrafast thermal response of band gap materials upon optical excitation. A model was established by extending the conventional two-temperature model that is adequate for metals, but not for semiconductors. It considers the time- and space-dependent density of electrons photoexcited to the conduction band and accordingly allows a more accurate description of the transient thermal equilibration between the hot electrons and lattice. Ultrafast thermal behaviors of bismuth, as a model system, were demonstrated using the extended two-temperature model with a view to elucidating the thermal effects of excitation laser pulse fluence, electron diffusivity, electron-hole recombination kinetics, and electron-phonon interactions, focusing on high-density excitation.

## INTRODUCTION

With increasing fundamental interest in ultrafast far-from-equilibrium dynamics and increasing practical demand for laser machining, phase-change memory devices such as GST (Ge_{2}Sb_{2}Te_{5}), and solar thermoelectric materials, it is of fundamental importance to understand light-matter interactions, thermal equilibration between non-equilibrium electrons and phonons, and their heat transport at different excitation densities and time scales.^{1–3} An ultrashort laser pulse can rapidly heat the electrons in metals to an initial temperature much higher than that of the lattice because the heat capacity of electrons is significantly lower than the heat capacity of the lattice.^{4–7} Thermalization between the electrons and the lattice is usually achieved via electron-phonon scattering within several picoseconds of photoexcitation. The conventional two-temperature model has been used successfully to describe the time evolution of the electron and lattice temperatures in photoexcited metals.^{4–9} These models, however, cannot describe the electron-lattice thermalization in bandgap materials accurately for several reasons. First, the density of photoexcited electrons in the conduction band is not constant in bandgap materials, due to electron-hole recombination, in contrast to the high carrier density Fermi gas in metals. Radiative and nonradiative recombination kinetics influence the rate and amount of energy flow from excited carriers to the lattice. The carrier heat capacity is proportional to the electron temperature in metals, but it must also be treated differently when the density of photoexcited electrons is low in band gap materials.

In semiconductor physics, diffusion models have described carrier diffusion, cooling, and recombination of photoexcited carriers.^{10–12} In semiconductors, excitation near the band edge deposits most of the pulse energy into creation of electron-hole pairs, so the carrier temperature contains a relatively small fraction of the pulse energy. Therefore, in semiconductors, carrier-lattice thermalization is less important than nonradiative carrier recombination as a source of lattice heating when exciting near the band edge. Semimetals lie somewhere in between and require a model that can bridge the regimes of low and high equilibrium carrier densities and small and large band gaps. Such models have been presented before^{13,14} to describe x-ray and electron diffraction results in photoexcited bismuth.

In this article, an extended thermal model is presented to investigate the ultrafast thermal behavior of band gap materials by extending the conventional two-temperature model to include materials with non-constant carrier density. Using photoexcited bismuth as an example, it describes the time evolution of the thermalization process between the photoexcited electrons and the lattice, highlighting the effects of laser pulse fluence, electron diffusion, electron-hole recombination kinetics, and electron-phonon coupling, including their roles in phase transitions of photoexcited materials.

## HEAT AND MASS TRANSPORT EQUATIONS

We derive a heat equation from the one-dimensional heat transfer process in which particles diffuse and carry their thermal energy from *z* to *z* + d*z*. *j* (particles/cm^{2} s^{−1}) indicates particle flux, which measures how many particles cross per unit area *A* per unit time Δ*t*. The heat accumulated in volume *A*Δ*z* per unit time is given by

where *Q* (J/cm^{3}) indicates the heat density and *C* represents the heat capacity of the particles, which may be dependent on the local temperature. The 3-D spatial problem can be reduced to a 1-D problem because the skin depth (∼15 nm) of bismuth at our 800 nm photoexcitation wavelength is several orders of magnitude smaller than the diameter of the excitation laser pulse (∼100-200 *μ*m). Thus, heat flow and electron diffusion occur predominantly perpendicular to the surface, where there are large temperature and density gradients. Hereafter, we will be concerned with the heat and mass transport along the perpendicular direction (z) only.

As Δ*t* and Δ*z* approach zero, Eq. (1) is transformed into a linearized differential equation,

Using Fick’s first law, $j=\u2212D d n d z $, Eq. (2) is rewritten as

where *n* denotes the density of particles (cm^{−3}) and *D* the diffusivity (cm^{2}/s), respectively.

Equation (3) describes the heat flow that is driven by the diffusion of particles. Depending on the electron kinetic energy, electron-phonon coupling time, and mean free path, diffusive transport can break down and Eq. (3) needs to be modified to include ballistic transport. In our thermal modeling, however, we assume that carrier transport is diffusive and not ballistic. Hereafter, particles represent the hot electrons photogenerated in the conduction band, indicated by the subscript *e*,

If there is a heat source ($ Q \u0307 e , + $) and a heat sink ($ Q \u0307 e , \u2212 $), these are included as follows:

The heat source is the additional energy introduced to the electronic subsystem, including the energy deposited by the ultrashort optical laser pulse. The heat sink represents heat loss from the electronic subsystem by electron-phonon coupling or electron-hole recombination. Equation (5) is the generalized electronic heat equation.

The heat flow in the lattice subsystem is simpler than that in the electronic subsystem because the diffusion of atoms is negligible. The conventional heat equation describes the lattice subsystem,

where *n _{i}* and

*C*represent the density of the lattice atoms and the heat capacity of the lattice, respectively. Both are independent of space and time coordinates.

_{i}*K*is the thermal conductivity of the lattice. The heat source and heat sink in the lattice subsystem are included as $ Q \u0307 i , + $ and $ Q \u0307 i , \u2212 $, respectively. Energy flow from the electronic subsystem to the lattice subsystem via electron-phonon coupling and nonradiative recombination is described by the heat source term.

_{i}The spatial and temporal distribution of electrons is determined by diffusion and recombination processes,

where the first term represents the diffusion process and the second term describes the first-order recombination process with rate constant *k* = 1/*τ*, where *τ* is the photoexcited electron (or hole) lifetime.

Equation (5) can be simplified through several steps. Since *Q* = *nCT*,

Equation (7) was used when propagating Eqs. (9) and (10). Comparison between Eqs. (5) and (10) yields

In the case of bismuth, the fraction of valence electrons that are photoexcited into the conduction band are on the order of 1% for a laser fluence of ∼5 mJ/cm^{2} at 1.55 eV photon energy^{15,16} and each photoexcited electron has high initial kinetic energy since the photon energy greatly exceeds the bandgap. Under this situation, every photoexcited electron can be considered as a free electron and they have the kinetic energy of $ 3 2 k B T e $. The corresponding heat capacity *C*_{e} is $ 3 2 k B $, independent of *T*_{e} as electron cooling occurs. In contrast to the electronic specific heat capacity of a metal, which is proportional to the temperature due to the high Fermi temperature, *C*_{e} will be assumed to be constant in Eq. (11) which becomes

As mentioned earlier, the term $ Q \u0307 e , \u2212 $ describes energy loss via electron-phonon coupling and electron-hole recombination,

where *g* is the electron-phonon coupling constant and *ξ* is the fraction of electron-hole recombination that occurs nonradiatively, thereby contributing to lattice heating, while (1 − *ξ*) is the quantum yield. Inserting Eq. (13) into Eq. (12) leads to

The energy gain in the lattice subsystem is

The energy loss term $ Q \u0307 i , \u2212 $ which appears in Eq. (6) was ignored in Eq. (16) because nonradiative phonon emission (which is a main channel for lattice energy loss) is negligible. Equations (7), (14), and (16) are coupled equations whose solution yields the time and space dependent carrier density, carrier temperature, and lattice temperature.

## PHYSICAL CONSTANTS, INITIAL CONDITIONS, AND BOUNDARY CONDITIONS

Electron transport plays a critical role in heat transfer from the hot electrons to the lattice. It has been reported that the electron diffusivity is as high as 100 cm^{2}/s for the ground state of bismuth at room temperature.^{17} A smaller diffusivity of 2.3-5 cm^{2}/s was also reported.^{18–20} Recent measurements yielded an intermediate value of 30-70 cm^{2}/s (3000-7000 nm^{2}/ps),^{13} which is used in this paper. These values are large enough that a considerable fraction of the photogenerated hot electrons move out of the excitation skin depth within 1 ps after photoexcitation, resulting in a far lower initial lattice temperature than would be reached in the absence of electron diffusion. It is also clear that the time evolution of the lattice temperature will be significantly different in bismuth films whose thickness is comparable to the skin depth, with an insulating substrate that suppresses carrier diffusion, than in a bulk bismuth sample whose surface is irradiated.

Inelastic scattering of hot electrons by phonons is responsible for lattice heating. This is the major mechanism for lattice heating in metals after impulsive photoexcitation. In the conventional two-temperature model, the electron-phonon coupling constant, *g*, is given by

where *m _{e}* is the electron effective mass,

*n*is the electronic density,

_{e}*v*is the speed of sound in the material, and

_{s}*τ*is the electron-phonon collision time.

_{ep}^{4,5}Since

*τ*∼ 1/

_{ep}*T*,

_{i}*g*is a constant independent of temperature. It is assumed for the present simulation that the electron-phonon collision time is

*τ*= 50 fs at 290 K and the free electron mass is used for

_{ep}*m*. The speed of sound is

_{e}*v*= 1.79 nm/ps in bismuth. Substitution of these values in Eq. (17) yields

_{s} where *n _{e}* is the time-dependent density of electrons excited in the conduction band in #/nm

^{3}. Eq. (17) is appropriate for estimating the constant of metals, but not for photoexcited band gap materials. Arnaud and Giret

^{14}measured the electron-phonon coupling constant in bismuth using the Debye-Waller factors and found it to vary between 10

^{−27}and 10

^{−23}J/(ps K nm

^{3}) for electron temperatures from 0 to 4000 K. Bismuth’s indirect bandgap and strongly varying density of states near the Fermi level lead to a strong change in the electron-phonon coupling as the carriers collect in the valleys at the T and L points.

^{13,14}In our model, we scale the electron-phonon coupling constant with the carrier density rather than temperature using Eqs. (17) and (18). Because of the uncertainties in the parameters such as the electron-phonon collision time and the electron effective mass in Eq. (17), we vary the proportionality constant,

*g*

_{0}, widely. Using an electron-phonon coupling constant of

*g*= 3.3 × 10

^{−24}

*n*J/(ps K nm

_{e}^{3}) produces the same heating rate as calculated by Arnaud and Giret,

^{14}consistent with experimental determination of the lattice heating rate. In this model, we consider a single temperature of the lattice that is coupled with carriers through the carrier-density-dependent coupling constant. This is a good approximation at low excitation fluence but can be treated in more sophisticated ways for high fluence. Large-amplitude lattice motion dramatically increases the effective electron-phonon coupling. Displacive excitation accelerates a bismuth atom toward the center of the unit cell along the A

_{1g}mode coordinate. At high amplitudes, this lattice motion results in strong coupling to the E

_{g}lattice mode ultimately culminating in nonthermal melting,

^{21,22}where the kinetic energy deposited into the lattice via displacive excitation is sufficient to melt the lattice.

^{23}This type of coupling can be treated by a three-temperature model

^{24,25}where the A

_{1g}mode is treated separately from the rest of the phonon bath, and the coupling to the bath is represented by a phonon-phonon coupling constant.

The specific heat capacity of bismuth is 25.52 J/mol K^{−1} at room temperature, which is close to the classical limit of the heat capacity for solids, 3*R*. This results from the low Debye temperature of 86.5–112 K for this material.^{26} In this work, the value for the lattice heat capacity of bismuth is taken as *C _{i}* = 3

*R*.

The thermal conductivity of bismuth is 8 W/m K. Since the electron density in the conduction band is low, the electronic contribution to the thermal conductivity is negligible. The value 8 W/m K is used for the lattice conductivity, *K _{i}*, in this paper. The corresponding thermal diffusivity,

*D*=

_{i}*K*/

_{i}*n*, is 6.7 nm

_{i}C_{i}^{2}/ps. This is very small compared to the carrier diffusivity

*D*∼ 7000 nm

_{e}^{2}/ps, discussed in the section on parameter selection.

The pump laser pulse is treated as a delta function instead of a Gaussian function with a finite width. Equivalently, it can be assumed that electrons are excited to the conduction band. This is a good approximation because our main concern is to examine the lattice thermalization and cooling, which occur after the end of the pump pulse. Lattice thermalization occurs in several picoseconds and lattice cooling proceeds slowly afterward.

The initial condition for the electron density in the conduction band is given by

where *F* is the pump laser fluence, *σ* is the skin depth, R is the reflectivity, and *n*_{g} is the equilibrium conduction electron density under ambient conditions (*n*_{g} = 0.001% of valence electrons at 300 K in bismuth)^{27} and *hν* is the photon energy. The initial density of electrons is assumed to vary linearly with the pump laser fluence up to our highest value of 60 mJ/cm^{2}.

The photoexcited electrons have initial kinetic energy given by

where *n*_{0} is the excited electron density. In principle, this is a function of depth *z* into the sample, but for the excitation densities, we are concerned with *n*_{0}(*z*) ≫ *n*_{g} even for *z* = 5*σ*, so we can use the approximate result everywhere.

In order to solve three (coupled) second-order differential Eqs. (7), (14), and (16), six boundary conditions are needed, in addition to the three initial conditions. The boundary conditions are given by

where *z* = 0 corresponds to the surface of bismuth and *z _{L}* indicates the thickness of bismuth used for numerical calculations (600 nm). The thickness of 600 nm is large enough to simulate bulk bismuth. Three Neumann boundary conditions are introduced at

*z*= 0. In the case of photoexcitation of a thin film, the left side boundary conditions can be modified to be the same as the right side boundary conditions, confining all carriers to the film.

## NUMERICAL RESULTS AND DISCUSSION

Coupled heat and diffusion Eqs. (7), (14), and (16) were solved numerically. Figure 1(a) shows the lattice temperature as a function of depth to the surface normal. The maximum temperature at *z* = 0 is reached at 18 ps and deeper regions show time delays for reaching a maximum temperature. At deeper regions, lattice heating by hot electrons through electron-phonon coupling and recombination is less significant because hot electrons lose most of their energy in the skin depth and disappear due to recombination within a characteristic time of 29 ps. Lattice heating at deeper regions is caused mainly by slower, phonon-mediated heat conduction. This is supported by the observation that the time for reaching the maximum temperature increases rapidly with increasing depth. The fluence of 10 mJ/cm^{2} raises the bismuth lattice temperature by 150 K, approximately halfway to the melting point. The present lattice temperature change at the surface coincides with transient reflectivity of bismuth films in recent ultrafast measurements,^{28} in which reflectivity signal grows rapidly until approximately 20 ps and then it decays very slowly.

Figure 1(b) shows the time evolution of the lattice and electron energy densities, integrated over the sample thickness (600 nm). The sum of these energies, which constitutes the total energy deposited into the system, is also shown. The lattice energy increases with time while the electron energy decreases. The characteristic time for the energy exchange is about 15 ps, which is determined by the electron-phonon coupling constant and the recombination rate. It was assumed that electron-hole recombination is entirely nonradiative (*ξ* = 1.0), as bismuth is an indirect bandgap material and does not luminesce. Therefore, there is no decay of energy as recombination occurs.

Figure 2(a) represents the density of photoexcited electrons integrated over the sample thickness (600 nm). The electron density decreases with time because of electron-hole recombination with a lifetime assumed to be 29 ps.^{13} Figure 2(b) shows the time evolution of the electron density at several sample depths, revealing the time scale for electron diffusion into the crystal.

Figure 2(c) shows the time-dependent electron temperature, which as discussed above is essentially independent of depth. If the band gap of bismuth is assumed to be 0.2 eV, 1.35 eV can be imparted to increase the kinetic energy of the electrons upon photoexcitation at 800 nm (1.55 eV). We can estimate the initial electron temperature to be approximately 16 000 K as shown in Fig. 2(c). It was observed from Fig. 1(b) that about 90% of the electron energy is lost to the lattice during the first 10 ps after photoexcitation. However, over the next 10 ps, a far smaller percentage of the remaining electron energy is lost. The reason is evident from Fig. 2(c). The electron-phonon coupling constant *g* is far higher at high electron density n_{e} than at moderate *n _{e}*. Thus, the electron temperature drops from about 16 000 K to 2000 K in the first 10 ps but is further reduced by less than 1000 K in the next 10 ps. Figure 2(d) shows the temporal local density of the electron energy which considers both the temperature and the electron density. Time-resolved ARPES (angle-resolved photoemission spectroscopy) experiments have found time scales of 5-6 ps for the time scale of electron-phonon coupling at comparable fluences,

^{29}which is comparable with our carrier cooling rate. Figure 2(e) shows the local density of carriers at several delay times where the density gradient reduces with time because of diffusion.

Figure 3(a) represents the effect of pump laser fluence on the lattice temperature. The maximum lattice temperature is proportional to the laser fluence. As discussed in the section on phase transitions, increasing the pump laser fluence increases the initial number of hot electrons excited in the conduction band without raising the electron temperature.

The effect of electron diffusivity on the lattice temperature is shown in Fig. 3(b). For the lowest diffusivity *D _{e}* value, 200 nm

^{2}/ps, considered in this work, the peak temperature change of the surface (

*Δ T*) is almost 370 K at 18 ps and is followed by a fast decay. In contrast, for the highest diffusivity value, 10 000 nm

_{l}^{2}/ps, the peak temperature change becomes about 90 K at 27 ps and is followed by a very slow decay. For the higher diffusivity, the initially deposited energy is moved far deeper into the crystal by the electrons before thermalization with the lattice occurs, reducing the peak lattice temperature and also reducing the thermal gradient that drives subsequent cooling.

The effects of the electronic diffusivity indicate substantially different thermal behaviors between bulk crystals and thin films of a material on an insulating substrate. In the bulk sample, hot electrons diffuse quickly to deeper regions, distributing their thermal energy over a large volume. In thin films, the diffusion of hot electrons is suppressed at the interface with an insulating substrate, resulting in a higher lattice temperature than that reached in the bulk sample.

The electron-phonon coupling constant, *g*, is given by Eqs. (17) and (18). As indicated earlier, there is considerable variation in the reported values and it is believed that the parameter cannot be considered a “constant” but rather has a value that varies considerably as a function of the electron density *n _{e}* and temperature

*T*

_{e}. In order to gauge the impact of this variation, we show in Fig. 3(c) the lattice temperature change for different values of the proportionality constant,

*g*

_{0}. As the value increases, the peak lattice temperature increases and is reached more quickly. Therefore, faster and more localized lattice heating is observed with increasing electron-phonon coupling constant. An electron-phonon coupling constant near

*g*

_{0}= 3.3 × 10

^{−24}

*n*yields good agreement with experimental results.

_{e}^{13,14}This value shows that most of the lattice heating occurs through electron-phonon scattering rather than electron-hole recombination. By the time recombination occurs, most of the carrier kinetic energy has already been transferred to the lattice. Figure 3(d) emphasizes this result as well. The lattice heating does not depend strongly on the carrier lifetime because electron-phonon scattering is faster than any realistic lifetime values.

## HIGH FLUENCE AND PHASE TRANSITIONS

We can estimate the pump laser fluence needed to reach the melting temperature *T*_{m} = 545 K(Δ*T _{i}* = 245 K) of bismuth at the surface, at first neglecting the effects of electron diffusion in which case all of the initially deposited energy in the electron subsystem are transferred to the lattice in the initially excited sample depth, i.e., the skin depth. In this case, the lattice temperature would be given by

where *C _{i}* is the heat capacity per unit volume. The result is

*F*≃ 1.8 mJ/cm

^{2}. In order to transform the phase of the material from solid to liquid at its melting point, we need a higher pump fluence to provide the enthalpy of melting Δ

*H*

_{m}= 545 (J/cm

^{3}), yielding an additional fluence of $F\u22433.2 m J / cm 2 $ for a total of 5 mJ/cm

^{2}which is close to the experimental results obtained from free-standing thin film measurements by Sciaini

*et al.*

^{23}where diffusion is suppressed so spatial transport is not significant.

According to the numerical results shown in Fig. 3(a), a pump fluence of about 28 mJ/cm^{2} is necessary to reach the melting point of bismuth. The large discrepancy between this value and 1.8 mJ/cm^{2} obtained above originates from the assumption that the pump laser energy is deposited only into the lattice without transported into the sample depth. Similarly, in order to induce melting, a fluence of ∼90 mJ/cm^{2} would be required rather than 5 mJ/cm^{2} as indicated above. Several experiments have been carried out on thin films^{18,23} in which carrier diffusion is suppressed and they can be compared with the simulation results for a thin film as shown in Fig. 3(e). It shows the temporal evolution of the lattice temperature for a 50 nm thick film at several excitation fluences. Around 18 mJ/cm^{2}, the film exceeds the melting point and enthalpy of melting, and the lattice would melt. This value is close to the melting threshold reported by Sokolowski-Tinten *et al.*^{18} and is much lower than that of bulk bismuth. Figure 3(f) represents the threshold fluence required for the melting transition of thin films as a function of thickness. It increases linearly with increasing thickness for small thickness (<300 nm) and then increases slowly. It becomes almost constant as 90 mJ/cm^{2} for films thicker than 500 nm. On the other hand, it converges to 5 mJ/cm^{2} as the thickness approaches zero as discussed above. For simulations in Fig. 3(f), the diffusion coefficient of 7000 nm^{2}/ps was used. Depending on the coefficient, the thickness dependence of the threshold fluence would differ significantly. Clearly the effects of fast electron diffusion on photoinduced phase transition and associated photophysical processes can be dramatic.

## CONCLUSIONS

We have extended the two-temperature model for band gap materials by considering photoexcited transient electrons in the conduction band and examined the ultrafast thermal response of photoexcited bismuth thin films upon impulsive photoexcitation. Based on the model, we described the time- and depth-dependent electron density, electron temperature, and lattice temperature of the material elucidating the thermal effects of laser pulse fluence, electron diffusivity, electron-hole recombination kinetics, and electron-phonon interactions. The extended two-temperature model that we present here is adequate for describing ultrafast thermal behavior of photoexcited band gap materials, including solid-to-liquid phase transitions.

## Acknowledgments

This research was supported by a grant from the Office of Naval Research (No. N00014-12-1-0530) and the National Science Foundation (No. CHE-1111557). We would like to thank A.A. Maznev for helpful discussions.