Ultra-fast carriers relaxation in bulk silicon following photo-excitation with a short and polarized laser pulse

Silicon (Si) is a fundamental building block of semiconductors physics and microelectronics industry [1]. The miniaturization of Si-based devices to the nano-scale regime and the never-ending search for faster devices call for a deep understanding of the fundamental quantum-mechanical process that governs the ultra-short time dynamics of electrons and holes [2,3]. Most of the knowledge of the electronic and optical properties of Si remain, however, limited to the equilibrium regime. Only recently the development of ultra-fast laser pulses [4,5] has opened the opportunity to directly investigate the real-time dynamics in the non-equilibrium (NEQ) regime [6].

In real-time experiments the system is initially perturbed with a short laser pulse (the pump) followed by a second weaker pulse (the probe) that measures a specific physical observable like, for example, the absorption [7,8] or the photo-emission [9,10] spectra. The dynamics induced by the pump is, then, monitored by observing and analyzing the modifications induced in these observable by the presence of photo-excited carriers.

Despite the enormous experimental interest and the continuous development of more refined experimental techniques, the simulation methods are still based either on equilibrium first-principles approaches or on NEQ model Hamiltonians.

In the case of model Hamiltonians the relaxation paths can be calculated by using the NEQ Green's function [2,11,12] or the Monte Carlo [13] methods. The advantage of these approaches is that the modifications induced by the presence of photo-excited charges are correctly taken into account in the evaluation of the scattering transitions. However ad hoc parameters must be introduced to describe both the photo-excitation process and the specific material properties.

First-principles simulations are commonly performed by using time-dependent density-functional theory (DFT) [14] or equilibrium many-body perturbation theory [15]. In the first case the coupling with the laser pulse is described but the dissipative processes are neglected [16] or described in an empirical way [17]. In the second case the laser pulse is replaced by some ad hoc initial guess of the carriers distribution and the scattering rates are derived from the equilibrium quasi-particle (QP) lifetimes [18].

In this letter we demonstrate that an ab initio description of both the photo-excitation process and of the time dependence of the carrier scattering rates is essential for the correct interpretation of the experimental results. In particular we will reproduce the time evolution carriers in bulk Si, observed in a recent two-photons photo-emission (2PPE) experiment [9,10], without relying on any parameter. We will highlight and discuss the different scattering channels created by the pump excitation showing the existence of two different decay regimes: an ultra-fast regime due to transitions between energetically degenerated states, made possible by the symmetry breaking caused by the pump pulse; and a slower regime, where the carriers are taken to the minimum of the conduction band. In addition we will investigate the very fundamental problem of defining the lifetime of a photo-excited carrier. We will show that this definition differs from the equilibrium case, also at a very low carriers concentration. Moreover this lifetime time dependence will clearly mark the different regimes (fast and slow) of the carrier relaxation.

Our theoretical framework is based on a non-linear equation for the time-dependent occupation factors [19], f_i(t) ($i\equiv \lbrace n\textbf{k}\rbrace$ represents a general band, n, and reciprocal space point, k):

$$\begin{equation} \partial_t f_{i} (t) = \big |\partial_t f_{i} (t)\big|_{pump}+ \partial_t f_{i} (t)\big|_{\textit{relax}}. \end{equation} \tag{ 1 }$$

All ingredients of eq. (1) are calculated ab initio [20,21] by using DFT [22]. $\partial_t f_{i}(t)\big|_{\textit{relax}}$ could be derived by using a standard semi-classical Boltzmann approach [2]. Instead $\partial_t f_{i}(t)\big|_{\textit{pump}}$ requires an approach based on Green's functions. This is, indeed, the why we use NEQ Green's function theory, reducing the complex quantum kinetic Baym-Kadanoff equations to the full set of eqs. (1)–(5) we are going to present. The detailed derivation is described in refs. [19,23]. It is indeed possible to show that the quantum kinetic equations reduce, within a completed collision approximation [19,24], to a Boltzmann-like equation, with a smeared out energy conservation, plus a coherent term which describes the creation of carriers in presence of an external field.

Thus, in our approach, eq. (1) originates from the diagonal terms of the equation of motion for time-diagonal NEQ Green's function, i.e. the time-dependent density matrix: $n(t)=-i G^<(t,t)$. This is a matrix in the electronic basis ${n_{ij}(t)}$, whose diagonal elements can be interpreted [2] as NEQ occupation factors, ${n_{ii}(t)=f_i(t)}$. The general equation has the form

$$\begin{eqnarray} \partial_t n_{ij}(t)&= i \left[H^{QP}+\Delta\Sigma[n](t)+U(t),n(t) \right]_{ij} \nonumber\\ &\quad-\,\gamma^1_{ij}(t) n_{ij}(t) +\gamma^2_{ij}(t). \end{eqnarray} \tag{ 2 }$$

Here U describes, in the length gauge, the coupling with external pump field: $U_{ij}(t) = -e \mathbf{r}_{ij} \cdot \textbf{E}(t)$. It is chosen to reproduce the exact experimental setup. The dynamics is governed by the DFT Kohn-Sham [22] Hamiltonian, H^KS, corrected by equilibrium many-body effects, $\Delta E^{QP}$, with ${H_{ij}^{QP}=H_{ij}^{KS}+{\Delta E}_{ij}^{QP}=\delta_{ij}\epsilon_i^{QP}}$. The deviation from equilibrium of the "Hartree plus exchange correlation effects" are described by the ${\Delta\Sigma}[n]$ operator. The equation for $i\neq j$ exactly reduces, within the linear regime, to the well-known Bethe-Salpeter equation [15,23] if: i) Σ is evaluated within a "Hartree plus statically screened Fock" approximation and ii) we use ${\gamma^2_{ij}(t)=0}$ with ${\gamma^1_{ij}(t)=\eta}$, a constant de-phasing of the polarization. such equation accounts for excitonic effects and, thus, ensures a very accurate description of the interaction with the laser pulse.

The term $\partial_t f_{i}(t)\big|_{\textit{pump}}$ originates from the commutator on the r.h.s. of eq. (2) and describes the creation of carriers by the interaction between the polarization and the external field. Thus, the equation for $f_i(t)$ is coupled to the one for n_ij, i.e. the polarization ${\textbf{P}(t)=-e \sum_{ij} \mathbf{r}_{ij} n_{ij}(t)}$, and the two must be propagated together.

The term $\partial_t f_{i}(t)\big|_{\textit{relax}}$ originates from the last two terms on the r.h.s. of eq. (2). It describes the relaxation and dissipation of the photo-excited carriers. It ensures that the carriers dissipate energy and scatter with phonons and electrons in such a way to relax towards the lowest energy states and has the form

$$\begin{equation} \partial_t f_{i}(t)\big|_{\textit{relax}} =-\gamma^{(e)}_{i}(t) f^{(e)}_{i}(t)+\gamma^{(h)}_{i}(t) f^{(h)}_{i}(t), \end{equation} \tag{ 3 }$$

with ${f^{(e)}_{i}=f_{i}}$ and ${f^{(h)}_{i}=1-f_{i}}$ the electron and the hole occupations. In this case the dynamics is fully dictated by the electron and hole lifetimes, ${\gamma^{(e)}}$ and ${\gamma^{(h)}}$, include both an electron-phonon (e-p) and an electron-electron (e-e) contribution: ${\gamma^{(e/h)}=\gamma^{(e/h)}\big|_{e-e}+ \gamma^{(e/h)}\big|_{e-p}}$.

In the e-p channel we have

$$\begin{eqnarray} &\gamma_{i}^{(e/h)}\Bigl|_{e\text{-}p}\propto\sum_{j\lambda}\sum_{I=\pm} N^{\lambda}_I\left| g^{\lambda}_{i\rightarrow j}\right|^2 \nonumber \\ &\times\,P\left(\omega_{\lambda}\pm \Delta_{ij}\right) f^{(h/e)}_j(t). \end{eqnarray} \tag{ 4 }$$

while in the e-e channel

$$\begin{eqnarray} &\gamma_{i}^{(e/h)}\Bigl|_{e\text{-}e}\propto\sum_{j k l}\left|W^{l\rightarrow k}_{i\rightarrow j}\right|^2P(\Delta_{ij}-\Delta_{kl})\nonumber\\ &\times\,f^{(h/e)}_j(t) f^{(h/e)}_k(t) f^{(e/h)}_l(t). \end{eqnarray} \tag{ 5 }$$

Here $\Delta_{ij}=\epsilon^{QP}_i -\epsilon^{QP}_j$ and λ represents a generic phonon mode with momentum q and branch η. $N^{\lambda}_{+}=(1+N_\lambda(T)/2)$ and $N^{\lambda}_{-}=N_\lambda(T)/2$, with $N_\lambda(T)$ the Bose distribution function at energy $\omega_\lambda$ and temperature T. W is the statically screened Couloumb interaction. The P functions represent a smeared energy conservation condition. Finally, $g^\lambda$ are the screened ionic potential derivatives, calculated within density-functional perturbation theory [25,26]. We have verified that, in the present work, the e-e scattering channel have only a negligible effect on the carrier dynamics due to the very low density of the photo-excited electrons. This is in agreement with the results of recent simulations [18,27,28].

Equations (3)–(5) make clear the different role played by ${\gamma^{(e)}}$ and ${\gamma^{(h)}}$. They describe the elemental process where an initial electron (hole) is scattered in another electron (hole) emitting or absorbing an electron-hole pair (e-e channel) [3] or a phonon [3,29,30]. In the e-p case the energy is transferred back and forth from the electronic to the phonon sub-systems until a thermal equilibrium is reached [2]. Thus, eq. (3) describes both relaxation and dissipation. In the $(e)$ channel the term ${\gamma^{(e)}f^{(e)}}$ describes the removal of electrons from the state $(i)$ and gives a negative contribution to $\partial_t f^{(e)}_{i}(t)$, while ${\gamma^{(h)}f^{(h)}}$ describes removal of holes, thus the filling of the state $(i)$, and gives a positive contribution.

In the 2PPE experiment we aim at describing [9,10], a Si wafer, oriented both along the [111] and the [100] surface directions, is excited with a laser pulse at room temperature. The photo-excited sample is, then, probed with a second laser pulse that photo-emits in the continuum the excited carriers. The photo-emitted current of electrons is measured as a function of the time delay between the pump and the probe. The final result is a measure of the time-dependent occupation of the valence bands (represented by the dots in the main frame of fig. 1). More specifically we consider the population of carriers near the point L₁, i.e. at $E\approx 1.6\ \text{eV}$ above the Fermi level¹ .

Zoom In Zoom Out Reset image size

We thus consider the electronic real-time dynamics in bulk Si, under the action of a laser pulse, whose parameters are taken directly from the 2PPE experiment [9,10]. The pulse is centered at around $3.4\ \text{eV}$, with duration (the full width at half-maximum) of $\sigma=110\ \text{fs}$ and intensity of $10^4\ \text{kW/cm}^2$ corresponding to an electric-field intensity of $6\times 10^6\ \text{V/m}$. The total fluence is $4\times 10^{-8}\ \text{KJ/cm}^3$, which means that the pump field creates a carrier density of about $4.5\times 10^{18}\ \text{el/cm}^3$, $\approx 1.8\ 10^{-4}\ \text{el}/\Omega$ where Ω is the unit cell. As the laser transferred momentum is negligible, on the scale of the solid unit cell size, all pumped carriers are excited vertically from the valence to the conduction bands along the $\Gamma\text{-}L$ line.

In fig. 1 the experimental occupation of the L₁ state (dots) is compared with the solution of eq. (1) (continuous line). The agreement between theory and experiment is excellent. Both the gradual filling and emptying of the L₁ state follows quite nicely the experimental curve. The theoretical results correctly describe the ultra-fast decay time-scale $(\sim 180\ \text{fs})$ and the $40\ \text{fs}$ shift of the population peak from the maximum of the pump pulse.

The $40\ \text{fs}$ delay reflects the delicate balance between the photo-excitation and the e-p scattering and can only be described by treating both processes on the same footing.

Experimentally [9,10] the ultra-fast decay of the L₁ state is interpreted as due to $L_1 \rightarrow X_1$ transitions. However, a deeper analysis of the theoretical result reveals a different scenario.

In fig. 1(a), the population of the levels at t = 0 is shown. Blue lines represent charges added and red lines charges removed. The band structure is computed along the $L \Gamma L^{\prime}$ high symmetry path in the Brillouin zone (BZ). In bulk Si L and $L^{\prime}$ are equivalent points but fig. 1(a) shows that the level $L_1^{\prime}$ is not populated and most of the carriers are injected in the L₁ level.

This symmetry breaking mechanism is made possible by the external field (U(t) operator in eq. (1)) which, in the 2PPE experiment, is polarized along the crystallographic [111] direction. This breaks the $L \leftrightarrow L^{\prime}$ symmetry as the operation that moves L in $L^{\prime}$, although being a symmetry of the unperturbed system, does not leave the [111] direction unchanged (see also fig. 2). In practice this means that eq. (1) does not respect this symmetry anymore and k-points connected by a rotation that does not leave the pumping field unchanged are populated in a different way [31]. Electrons are injected in the conduction band along the $\Gamma\text{-}L$ line but not, for symmetry reasons, along the $\Gamma\text{-}L^{\prime}$ line. This is clearly shown in fig. 1 where the population of the $L_1^{\prime}$ state is represented with a dashed line. The $L_1^{\prime}$ state is gradually filled while L₁ is depleted revealing that the real source of the ultra-fast decay observed experimentally is the $L_1\rightarrow L_1^{\prime}$ scattering.

Zoom In Zoom Out Reset image size

This scattering is faster than any other scattering as it involves states with the same energy. Indeed fig. 1 shows that initially the dynamics equilibrates the populations of the L₁ and $L_1^{\prime}$ states which reach the same value at $t\sim 220\ \text{fs}$. After this point both states decay simultaneously by using the slower $L \rightarrow X$ channel towards the conduction band minimum (CBM)² . After $t\sim 500\ \text{fs}$ the electrons (holes) can be already described by two Fermi distributions around the CBM VBM (=valence band maximum) with very high temperatures ($T_e\approx 9000\ \text{K}$ for electrons and $T_h\approx 2350\ \text{K}$ for hole). Once the Fermi distributions are created, the relaxation process is mostly dissipative and phonons are emitted in order to cool the carriers temperatures (after $t=1\ \text{ps}$, for example we obtain, $T_h\approx 550\ \text{K}$ and $T_e\approx 1900\ \text{K}$).

To better disentangle the $L_1\rightarrow L_1^{\prime}$ process from the slower $L \rightarrow X$ channel we consider a shorter laser pulse with $\sigma=50\ \text{fs}$. We also consider a gedanken experiment where electrons are manually excited. One of the approximations most widely used in the literature is to mimic the effect of the laser pulse with some, ad hoc, initial population of carriers in the valence bands. This approximation corresponds to put U = 0 in eq. (2) defining some initial arbitrary population $f_{i}(t=t_0)$. Here we have chosen an initial population around the $\Gamma_{15}$ state, with a carrier density equal to the one measured experimentally.

We then show (fig. 1(b)) the population of the L₁ state in the gedanken experiment (dot-dashed line), and when the photo-excitation is performed with the shorter pulse for both L₁ (continuous line) and $L_1^{\prime}$ (dashed line). By comparing the results we notice that the decay of the L₁ state, when the carriers are manually excited is much slower compared to the case when the carriers are photo-excited. This is because the dynamics following the ad hoc population is symmetric and the L₁ and $L_1^{\prime}$ are equally populated. The ultra-fast $L_1\rightarrow L_1^{\prime}$ channel is switched off. Instead, when the shorter pump pulse is considered, the two decay time-scales (the ultra-fast $L_1\rightarrow L_1^{\prime}$ and the slower $L \rightarrow X$) are clearly visible. Immediately after the population peak we notice an almost vertical drop of the L₁ population. The characteristic time scale is even faster $(90\ \text{fs})$ than the one measured experimentally $(180\ \text{fs})$, and still dictated by the length of the laser pulse. After the vertical drop we notice a smoother decay of both the L₁ and the $L_1^{\prime}$ occupations, induced by the much slower $L \rightarrow X$ channel.

The most time-consuming part of solving eq. (1) is the update of the $\gamma^{(e,h)}(t)$ functions whose dependence on the occupations must be re-calculated at each time step. A very tempting possibility would be to keep $\gamma_{i}^{(e,h)}(t)$ constant. This is indeed the main ingredient of the relaxation time approximation (RTA) that is based on the assumption that $\gamma_{i}^{(e,h)}(t)\sim\gamma_{i,eq}^{(e,h)}$, with $\gamma_{i,eq}^{(e,h)}$ the equilibrium lifetimes, calculated without the presence of any external field³ . This approach has been recently used in ref. [18] to describe the carrier relaxation in Si excited by the weak sunlight.

In order to investigate further the meaning of the $\gamma_{i}^{(e,h)}(t)$ time dependence and their crucial role in describing the experimental results, let us introduce the totally relaxed occupations, $f^{\infty}_{i}$, defined as the occupations at the time t_rel such that ${\partial_t f_{i}(t)|_{t=t_{rel}}=0}$. From eq. (1) and eq. (3) it follows that

$$\begin{equation} f^{\infty}_{i}\equiv f_{i}(t_{rel})=\frac{\gamma_{i}^{(h)}(t_{rel})} {\gamma_{i}^{(h)}(t_{rel})+\gamma_{i}^{(e)}(t_{rel})}, \end{equation} \tag{ 6 }$$

where we have used the fact that at $t=t_{rel}$ the pump field is switched off and $\partial_t f_{i}(t_{rel})\big|_{\textit{pump}}=0$.

In the equilibrium regime the dependence of the $\gamma_{i,eq}^{(e/h)}$ lifetimes on the electronic energies is well known. Indeed for conduction bands $\gamma_{i,eq}^{(h)}=0$ while for valence states $\gamma_{i,eq}^{(e)}=0$. This means that, from eq. (6), the totally relaxed occupations would be zero for any conduction band.

This is the reason why in the RTA [18] $f^{\infty}_{i}$ is added as an adjustable parameter in the simulation. It is commonly parametrized as a Fermi distribution with a given temperature and chemical potential. But our simulations reveal that, in general, electrons and holes are distributed with two different Fermi distributions. This means two different chemical potentials and temperatures. In our approach $f^{\infty}_{i}$ is a by-product of the simulation and it must not be provided at the beginning.

Moreover, the present scheme allows to go well beyond the RTA by formally defining a NEQ carrier lifetime, ${\overline{\gamma}^{(e/h)}_{i}(t)}$, such that occupation functions satisfy the simple equation ${\partial_t f^{(i)}=-\overline{\gamma}^{(i)} f^{(i)}}\ ({i=e/h})$ with

$$\begin{eqnarray} &&\overline{\gamma}^{(e/h)}_{i}(t)=\pm\left[\frac{ \gamma_{i}^{(e)}(t) f^{(e)}_{i}(t)-\gamma_{i}^{(h)}(t) f^{(h)}_{i}(t) }{f^{(e/h)}_{i}(t)}\right] \text{.} \end{eqnarray} \tag{ 7 }$$

In eq. (7) the + (−) corresponds to electrons (holes). Equation (7) demonstrates that a true instantaneous NEQ carrier lifetime includes contributions from both the electron ($\gamma_{i}^{(e)}$) and the hole ($\gamma_{i}^{(h)}$) lifetimes.

The deviation of $\overline{\gamma}^{(e/h)}_{i}(t)$ from $\gamma^{(e/h)}_{i,eq}$ is, indeed, strictly connected with the symmetry breaking mechanism that explains the experimental result. In the $L_1$/$L_1^{\prime}$ case, indeed, $\gamma_{L_1,eq}^{(e)}=\gamma_{L_1^{\prime},eq}^{(e)}\approx 20\ \text{meV}$. From fig. 3 we see that, instead, both $\overline{\gamma}^{(e)}_{L_1}(t)$ and $\overline{\gamma}^{(e)}_{L_1^{\prime}}(t)$ changes, during the simulation, by an order of magnitude. The reason is that, in eq. (7), the $\gamma^{(e)}_{i}f^{(e)}_{i}$ and $\gamma^{(h)}_{i}f^{(h)}_{i}$ factors are of the same order and their balance measures the difference of population between the L₁ and $L_1^{\prime}$ states. We can, therefore, easily recognize in fig. 3 two well-defined regimes: when $\overline{\gamma}^{(e)}_L \gg \overline{\gamma}^{(e)}_{L^{\prime}}$ the ultra-fast $L_1 \rightarrow L_1^{\prime}$ scattering channel is active. When $\overline{\gamma}^{(e)}_L \approx \overline{\gamma}^{(e)}_{L^{\prime}}$, instead, the relative L₁ and $L_1^{\prime}$ populations are balanced and the dynamics is dictated by the slower $L\rightarrow X$ channel.

Zoom In Zoom Out Reset image size

To further illustrate the results of our simulations we also show, in fig. 4, an overview of the carriers dynamics by four snapshots. The carriers distribution is represented here on the band structure⁴ of Si. Added (removed) carriers are in blue (red). At $t=0\ \text{fs}$ (panel (a)), corresponding to the laser pulse maximum, at $t=200\ \text{fs}$ (panel (b)), when the L₁ and the $L_1^{\prime}$ populations are almost the same, at $t\approx 500\ \text{fs}$ (panel (c)) when both the L₁ and $L_1^{\prime}$ states are almost empty and finally after more than $1\ \text{ps}$ (panel (d)) when the electrons have reached the CBM. The dominant process at each time considered is schematically represented by the arrows.

Zoom In Zoom Out Reset image size

In conclusion we have presented a fully ab initio simulation of the carrier dynamics in Si. The present scheme, based on the merging of DFT with NEQ Green's function theory, successfully describes the ultra-fast decay of the L₁ carrier population measured in a recent 2PPE experiment. We have also highlighted that the microscopic mechanism that drives this ultra-fast decay is not a standard inter-valley scattering but it is due to an ultra-fast (as fast as $90 \ \text{fs}$) $L_1\rightarrow L_1^{\prime}$ scattering channel activated by the specific polarization the pump laser. This physical interpretation is, also, supported by introducing a novel definition of the non-equilibrium carrier lifetime that provides an intuitive picture of the physical processes activated by the initial photo-excitation.

Financial support was provided by the Futuro in Ricerca grant No. RBFR12SW0J of the Italian Ministry of Education, University and Research. DS thanks G. Onida for the access granted to the etsfmi cluster in Milano and C. Attaccalite for useful discussions.