Energetics and Dynamics of Laser-Assisted Field Emission from Silicene Nanoribbons : Time-Dependent First-Principles Study

We investigate laser-assisted field emission (LAFE) from a silicene nanoribbon (SiNR) using time-dependent density functional theory (TDDFT) simulation. The emission mechanism in the present study is considered to be over-barrier photoemission and is found to be governed by electronic dipole transitions, the characteristics of the excited states, and the energy levels of the excited states relative to the potential hump. The qualitative features of emission from the SiNR are similar to those from graphene nanoribbons (GNRs). The emission currents from SiNR, however, are found to be much larger than those from GNR for the same laser parameters in spite of a larger work function for SiNR. We reveal the emission currents in real time and space on an atomic scale, and observed current being driven back and forth in the early stage of emission. We further elucidate the dynamical correlation among the laser pulse, the Kohn-Sham potential and emission currents wiggling under the ponderomotive force. [DOI: 10.1380/ejssnt.2015.115]


I. INTRODUCTION
Studies of ultrashort-pulse field emission (FE) from nanoscale metal tips using femtosecond laser pulses [1,2] have attracted great interest because of the possibility of spatiotemporal resolution down to an atomic scale with an attosecond range utilizing a coherent electron source.Recently, an ultrashort-pulse FE source with site selectivity on the scale of a few tens of nanometers was reported [3,4].Emission mechanisms have been investigated and categorized phenomenologically by several processes such as photofield emission, multi-photon photoemission (or above-threshold photoemission (ATP)) and optical field emission.Which of these processes is dominant depends on the value of the Keldysh parameter γ, which is determined by laser parameters and the work function of the emitter [5].Optical field emission begins to play an important role in attosecond physics with an increasing laser field [6][7][8][9][10].The electron rescattering with surface atoms [9,11] that results in a plateau and a cutoff in the photoemission spectrum [9,10] and the carrier envelope phase (CEP) [10,12,13] of few-cycle near-infrared laser pulses are typical phenomena in attosecond emission dynamics.In order to further understand the emission processes, different calculation methods have been presented in order to interpret the emission mechanisms.A modified Fowler-Nordheim formula has been utilized for optical field emission [1,2], and the Boltzmann equation has been solved for photofield emission [3,4,14].More recently, the finite-difference time domain method has provided a field strength distribution around the tip apex in 100 nm dimensions with the presence of the CEP [9].The authors of Ref. 9 also applied a 1D time-dependent density functional theory (TDDFT) [15] simulation of jellium substrate combined with 3D-quasiclassical Monte Carlo simulations to the time-dependence analysis of emitted electrons and their energy spectra in the ATP processes.
The objectives of our computational study are to understand the origin of the electronic states responsible for laser-assisted excitation and emission from nanostructures and demonstrate ultrafast electron dynamics in real time.We first perform density-functional theory (DFT) calculations of 3D-atomic nanostructures to determine the electronic band structures.Next, time evolution of electron emission from the atomic surface is demonstrated by TDDFT simulations to elucidate the dynamical correlations among the femtosecond laser pulse, Kohn-Sham (KS) potential and emitted electrons on a femtosecond time scale.The TDDFT method has become a standard tool for studying the electronic states of excitations and non-equilibrium processes for many-electron systems.Since the KS Hamiltonian includes electron-electron interaction (Hartree and exchange-correlation), electron-ion interaction, and the interaction of electrons with external fields (DC field and laser field) automatically, the time evolution of electron wavefunctions obtained from the solution of the time-dependent KS equation gives detailed information on excitation and emission dynamics.
In our previous study on this topic with TDDFT simulations, we clarified the microscopic mechanism of electron emission from graphene nanoribbons (GNRs) by an analysis of the electronic dipole transition matrix element and the electronic properties of the emitting states [16].We observed photoemission for the laser parameters used in the simulation and found that the emission properties strongly depend on the electronic structures of GNRs.
In the present study, we perform a simulation of laserassisted field emission (LAFE) from silicene nanoribbons (SiNRs) [17], which are similar to GNRs in electronic state and atomic structure.Silicene, which is a one atom thick silicon sheet arranged in a honeycomb lattice like graphene, has been shown to have properties different from graphene.Although the electrons in silicene were considered to behave as Dirac fermions [18], many studies have shown that silicene synthesized on Ag(111) surfaces does not have Dirac fermions due to the strong coupling with the substrate [19].Further, even freestanding silicene does not have a flat layer but a buckled one.(See a review article 20 for the recent studies on silicene).Therefore, it is of interest to compare the electron emission properties of SiNRs with those of GNRs, as done in this study.

II. METHOD AND MODEL
Throughout this section, we adopt atomic units in which ℏ = 1, m = 1 and e 2 = 1, unless stated otherwise.The time-dependent KS equation is given below as where ψ i (r, t) is the KS wavefunction of the valence electron, and i denotes the level index that runs over all occupied valence states.The Hamiltonian is constructed by where R α (t) is the position of the αth ion, V L (r − R α ) is the local pseudopotential and V l,m NL (r − R α ) is the nonlocal pseudopotential [21] for the αth ion at R α , formulated by Troullier and Martins [22] based on the Kleinman-Bylander separable form in real space [23].Pl is the projection operator on the angular momentum (l) and m is a magnetic quantum number that takes integral values from −l to l. µ xc is the exchange-correlation potential, and v ext (r, t) is the external field potential.
From the time-dependent KS wavefunctions, the electronic density ρ and emission current density J are respectively calculated as and where f i is the occupation number of the i-th KS level.In the present study, the time evolution operator is approximated by a fourth-order Taylor expansion method [24].The Ceperley-Alder functional [25] for local-density approximation with Perdew-Zunger parameterization [26] is used for the exchange-correlation interaction.The KS Hamiltonian includes the potentials of an electrostatic field and a laser field, which are simultaneously applied to the target nanostructure.
The calculation model is shown in Fig. 1.The target is chosen to be a hydrogen (H)-terminated zigzag-edged SiNR.Its atomic structure is similar to GNRs, as used in the previous simulation in Ref. 16, except for the buckled structure in Fig. 1(a).The optimized buckling distance of the Si atoms in the z-direction is 0.44 Å, which agrees well with previous calculations [18].The directions of the external fields, as well as the unit cell with the periodic boundary condition, are given in (b).Both of the electrostatic and laser fields are applied on the x axis.The unit cell size is set at 65.5×3.97×15.0Å3 .Twelve k y points were chosen along the ribbon axis in the first Brillouin zone.A complex absorbing potential (CAP) [27] is introduced in the vacuum to absorb emitted electrons, preventing them from being reflected at the boundary of the unit cell.
The emission current is measured in an observation plane placed at 10 Å above the H atoms, while three positions labeled by A, B, and C are chosen for monitoring the spatiotemporal evolution of the current densities.The distance from the H atom to position A is 2.7 Å, while the A-B and B-C distances are both 1.2 Å. Figure 1(c) shows the time evolution of the electric field applied to the SiNR, which consists of an electrostatic field of −0.075 V/ Å and a laser field with a strength of 0.05 V/ Å (power of 3.3 × 10 10 W/cm 2 ) and an energy of 6 eV.The laser field is spatially uniform and has a Gaussian envelope with a FWHM of 6 fs.

III. RESULTS
We first calculated the SiNR band structures in the absence of external fields together with the averaged KS potentials along the x axis under an electrostatic field of E DC = −0.075V/Å, as shown in Fig. 2(a) and (b), respectively.To compare, we also made the same calculations for the counterpart GNR with hydrogen-terminated zigzag edges and the results are shown in (c) and (d), which are essentially the same as Fig. 2(a) and (b) in Ref. 16, respectively except that the Kohn-Sham potential in Fig. 2(b) in Ref. 16 was put downward slightly due to a different energy origin.In all figures, the energies are measured from the Fermi energies.In both (a) and (c), flat bands near the X point, which correspond to the edge states, and quasi-continuum states above an energy of about 4 eV, which correspond to the vacuum levels, can be observed.Quantitatively, the energy-level spacings are smaller for the SiNR than for the GNR, which means that there might be more contributing channels to http://www.sssj.org/ejssnt(J-Stage: http://www.jstage.jst.go.jp/browse/ejssnt/) e-Journal of Surface Science and Nanotechnology the LAFE from the SiNR.The KS potentials in (b) and (d) are similar to each other, showing a linear decrease from greater than x ∼25 Å with a slope corresponding to the applied field E DC .The work function is larger for the SiNR (4.6 eV) than for the GNR (4.0 eV), which means that the FE should be more difficult from the SiNR than that from the GNR in the absence of laser fields.
To further clarify the band contributions to the FE currents from the SiNR, we calculated the dipole transition matrix elements, V d = ⟨i|x|f ⟩.These are the off-diagonal matrix elements of the dipole operator x between the initial states |i⟩ and the final states |f ⟩ corresponding to the energy bands in Fig. 2(a).The calculated matrix elements at the Γ point are shown in Fig. 3(a), since the elements at Γ are significantly larger than those at other k points.The red color shows the largest V d , and blue shows the smallest.The levels marked red do not make an important contribution to emission because the 26th level is below and 20th level is far below the potential hump, as seen in Fig. 2(b).In contrast, the 28th and 30th states are found to be emission states, although a V d (10 → 28) (solid circle, the energy difference between the two states is 6.8 eV) and V d (14→30) (dashed circle, the energy difference between the two states is 6.0 eV) do not seem large.The energy levels of the final states, i.e., the 28th and 30th levels at the Γ point, are sufficiently higher than the potential hump.
The corresponding KS wavefunctions for the final states are shown in Fig. 3(b).It can be seen that these levels have electronic structures with a σ-like character having electronic distribution along the direction parallel to the external fields and protruding toward the vacuum.According to our previous study, the relative angle between the direction of electronic distribution and that of applied fields was found to be a critical factor for determining the emission probability; electrons can be emitted efficiently from σ-like states (parallel to the electric field) but not from π-like states (perpendicular to the electric field).This has been verified both in conventional FE from graphitic nanostructures and LAFE from a GNR [16,28].Therefore, the orbital characteristics in Fig. 3(b) are important for the discussion below.In determination of the electronic states responsible for electron emission, we took into account the value of V d and the orbital characteristics of the emitting states.It would be necessary to evaluate the the time-dependent decomposition of Kohn-Sham orbitals that contribute electron emission for con- firmation of the emitting states.This will be the next study.
Having discussed the origin of the electronic states responsible for emission from the SiNR, we present the results of the emission currents, compared with those from the counterpart GNR in Fig. 4. The results for laser energies of 6.0 and 6.8 eV are shown; for energies of 5.0 and 5.5 eV, we observed emission currents two orders of magnitude smaller (not shown here).In spite of the larger work function of the SiNR (4.6 eV) compared with that of the GNR (4.0 eV), the emission current from the SiNR under the 6.8 eV laser is found to be about six times as large as that of the GNR.The present results can be interpreted by a combination of the factors discussed above: the V d , the energy level of the emitting states relative to the potential hump and the orbital characteristics of the emitting states.In addition, we observed that there are more emission channels (higher density of states) in the SiNR.
For the 6.8 eV laser, V d (10→28) makes a dominant contribution to the FE current for the SiNR (see Fig. 3(a), solid circle), while the transition from the 16th to the 20th level dominates the current from the GNR [16].The orbital characteristics for the 28th level of the SiNR enable the electrons to emit efficiently due to its protruding structure near the H atom.As shown in Fig. 2(a), the energy of the 28th level is sufficiently high such that the emitted electrons can efficiently go beyond the potential hump.Therefore, excitation from the 10th to the 28th level is a favored channel for SiNR electron emission.The number of emission channels also makes important contributions to the larger current of the SiNR.Under the 6.8 eV laser, electrons in the 13th to 17th levels can be excited to the quasi-continuum above the 31st level, where the transition probabilities are greatened by the high-density of continuum states, even with a small V d , and then the excited electrons can be emitted efficiently towards the vacuum.These are the reasons for the large current from the SiNR under the 6.8 eV laser with an emission mechanism considered to be over-barrier photoemission [10].
For the 6.0 eV laser, V d (14 → 30) (dashed circle in Fig. 3(a)), which looks small, is considered to be a possible channel of emission from the SiNR because the final state is sufficiently high relative to the potential hump.Compared with emission from the 28th level, that from the 30th level could be less efficient because KS wavefunction projection of the level near the H atom is smaller than that for the 28th level (see Fig. 3(b)), although the energy of the state is above the potential hump.In addition to the 30th level emission, excitation from the 17th level to the quasi-continuum also contributes to the current under the 6.0 eV laser.This explains the finite FE current from the SiNR under the 6.0 eV laser.In contrast, the GNR has no channels for emission under the 6.0 eV laser and thus the FE current is negligible.
To understand the spatiotemporal FE dynamics in more detail, we examine the time evolution of emitted current density at the points labeled A, B and C under the 6 eV laser as in Fig. 1.The corresponding results are shown in Fig. 5.At point B, the KS potential under the electrostatic field without the laser has a spatially maximum value (hump); this corresponds to x ∼28 Å in Fig. 2(b).The emission current density in Fig. 5 cillates and increases in amplitude initially, it reaches a maximum at ∼12.5 fs, and then decreases.The oscillation period corresponds to that of the laser.Similarly, the current density in (b) also shows an oscillating feature and reaches a maximum at ∼14 fs; however, the oscillation amplitudes in (c) are greatly reduced with a maximum at ∼17 fs.The oscillations are the strongest in (a) because position A is located closest to the ribbon edge, while the other two points are located farther away, where FE currents from different channels can overlap and cancel out the oscillations.This can also explain the delay in time at which the maximum is reached in (b) and (c).Another interesting feature observed in Fig. 5 is that the oscillating current density in (a) and (b) changes polarity at the initial time; the emitted electrons are driven back and forth every 0.69 fs.The average of the current density is positive, which indicates emission of electrons toward the vacuum.This shows the wiggle motion of the emitted electrons is caused by the laser field [29], and the polarity change in (b) implies that the electrons go back into the ribbon.by the solid and broken lines, respectively.The interval between the two times is close to a half period of the laser oscillation (0.34 fs); the slopes of the broken and solid lines beyond the region x = 30 Å are -0.039 and -0.124 V/ Å respectively.Here, it is noted in Fig. 6(a) that the KS potential hump does not decrease noticeably with increasing the time-dependent field strength.The feature is understood by the fact that the charge polarizations at both ends of the ribbon that is induced by the laser field hinder the electron emission by increasing the amount of the electric-dipole layers at ribbon edges.Snapshots of the FE current density distribution at the two times are shown in Fig. 6(b).In this short timescale, we cannot clearly see current propagation toward the vacuum.However, in Fig. 6(c), which are the snapshots over a longer time scale, propagation of the emission current toward the vacuum region is clearly seen, especially for t = 14.03 and 18.87 fs.As the laser has its peak at t = 12 fs (Fig. 1(c)), the electrons begin to emit after the laser peak (see the center panel of Fig. 6(c)) and are then driven by the electrostatic field to the vacuum.

IV. CONCLUSION
In conclusion, the energetics and dynamics of LAFE from nanoribbons were clarified by TDDFT simulations.The emission mechanism, which is common to the SiNR and GNR, is found to be governed by three factors: the electronic dipole matrix element, the orbital characteristics of the excited states, and the energy levels of the excited states relative to the potential hump.However, the emission currents from the SiNR are found to be much larger than from the GNR for the same laser parameters in spite of the larger work function of the SiNR.This can be interpreted by the higher density of states, or more channels for dipole transition in SiNRs than that of GNRs.We demonstrated the dynamical correlation among the laser pulse, the Kohn-Sham potential and emission currents wiggling under the laser field in real time and space on an atomic scale, and observed emission current being driven back and forth in the early stage of emission.

FIG. 1 .
FIG. 1.(a) Side view of hydrogen (H)-terminated SiNR.(b) Top view of H-SiNR under an electrostatic field (EDC) and a laser pulse (ELaser).The Si and H atoms are indicated by large and small balls, respectively.The currents are emitted toward the right vacuum region.The unit cell used in the calculations is shown by the dashed box.(c) Time evolution of sum of electrostatic field and laser field.

FIG. 2 .
FIG. 2. Energy band structures of (a) SiNR and (c) GNR.Integer numbers indicate the band indices.KS potentials under the electrostatic field for the SiNR and GNR are shown in (b) and (d), respectively.The dotted lines in figures (b) and (d) indicate the work function (upper line) in the absence of an electric field and the Fermi energy (lower line), respectively.Figures (c) and (d) are essentially the same as Fig. 2(a) and (b) in Ref. 16, respectively except for the details.

FIG. 3 .
FIG. 3. (a)Map of dipole matrix elements for SiNR at Γ point.The red color shows the largest V d , and blue shows the smallest.(b) Side view of the SiNR KS wavefunctions for the 28th and 30th levels at Γ are shown.

Figure 6 (FIG. 6 .
FIG. 6.(a) Averaged KS potentials for SiNR along x axis at t = 14.15 and 14.51 fs, denoted by solid and broken lines, respectively.(b) Corresponding FE current density distribution snapshots along x-axis for (b) short and (c) long time scales.