# Strain-controlled high harmonic generation in two-dimensional materials with Dirac fermions

###### Abstract

Two-dimensional (2D) materials with zero band gap exhibit remarkable electronic properties with wide tunability. High harmonic generation (HHG) in such materials offers unique platforms to develop electronic-engineered novel optoelectronic devices at nanoscale, as well as to investigate strong-field and ultrafast nonlinear behaviour of massless Dirac fermions. Here we report controllable HHG from monolayer silicene by tuning the electronic properties via strain engineering using an ab initio approach based on time-dependent density-functional theory (TDDFT). We show that under small biaxial and uniaxial strains, large anisotropy and anharmonicity can be induced in the band structure of silicene while preserving the Dirac cone, which can lead to significantly enhanced harmonic intensity and broadened spectral width. With the advantage of silicene in compatibility and integration into the current silicon-based electronic industry, this study may open a new avenue to develop tailored extreme-ultraviolet and attosecond pulse sources at nanoscale, and provide a valuable tool to understand the strong-field and mechanically induced ultrafast nonlinear response of Dirac carriers in 2D materials.

## Introduction

The recent active investigations on high harmonic generation (HHG) from solid-state materials has opened up exciting opportunities in both fundamental physics and potential applications[1, 2, 3, 4, 5]. It provides an unique platform for the study of strong-field and ultrafast electron dynamics in the condensed phase[6]. HHG has been demonstrated as an attracting tool to explore the electronic structure of bulk crystals[7, 8, 9, 10, 11]. Experimental reports such as all-optical reconstruction of band structure of ZnO[12], retrieving energy dispersion profile of the lowest conduction band of SiO[13], and measurement of the non-vanishing Berry curvature in symmetry-broken -quartz[14], etc., mark some of the impressive advancement. HHG in solids also offers a promising approach for the development of novel compact coherent light sources in the extreme ultraviolet and soft x-ray spectral region[15], as well as for multi-petahertz-frequency optoelectronic[16] and attosecond photonics[17]. High harmonics and/or attosecond pulses can be shaped in terms of polarization and carrier-envelop phase[18, 19], and enhanced by means of tailored semiconductors[20] and metallic nanostructures[21].

In addition to bulk solids, two-dimensional (2D) materials have also attracted considerable attention for HHG which exhibit distinctive electronic properties. Monolayer h-BN driven by out-of-plane laser fields has been calculated to exhibit atomic-like HHG while with a more favorable wavelength scaling[22]. Driven by in-plane fields, enhanced HHG efficiency from isolated monolayer MoS compared to the bulk has been measured[23]. Besides, there are a few theoretical and experimental investigations of HHG from 2D materials with massless Dirac fermions (MDF), though so far limited to graphene[24, 25, 26, 27, 28, 29, 30]. The study of nonperturbative HHG from such zero-gap 2D materials offers new possibilities to access strong-field and ultrafast nonlinear dynamics of MDF[31].

Silicene, as the 2D form of silicon and the silicon analogue of graphene, also exhibits exceptional electronic properties such as MDF behavior, linear energy dispersion, and a high Fermi velocity. In addition, as silicon being the element of IV family close to carbon, silicene has apparent advantage over graphene in compatibility and integration into the current silicon-based electronic industry. Thus, HHG from silicene, which has not been studied before, is highly relevant not only for investigating the remarkable properties of massless Dirac quasiparticles, but also for developing future integrated all-optical solid-state devices at nanoscale.

Being amenable to extensive nanoengineering technologies, the electronic properties of 2D materials can display very rich and complex phenomenons and wide tunability, which can have significant influence on HHG process. Control of HHG from MDF dynamics in 2D materials by tuning the electronic properties may thus open new perspectives for the realization of nanoscale devices and shed light on the ultrafast dynamics of MDF. To date, however, HHG control from 2D materials by tuning the electronic properties has remained largely unexplored. Moreover, the reported theoretical investigations so far are based on simplified models without properly considering the full electronic and real crystal structures. In this work, we report, to our knowledge, the first controlled HHG from monolayer silicene tuned by strain engineering using an ab initio approach based on the framework of time-dependent density-functional theory (TDDFT)[32, 33]. We use the methods established recently by Tancogne-Dejean et al[11, 18, 22]. The first-principles simulation results show that the electronic band structures of silicene and the HHG are sensitive to strains. Strain engineering thus can be a promising strategy for controlling HHG and studying MDF dynamics in 2D materials.

## Results

With a 2D honeycomb lattice structure like graphene, silicene is predicted to favor a low-buckled structure[34] (Fig. 1a, b). In this study, the geometric structures of silicene are obtained by using the Octopus package. The calculated lattice constant of the hexagonal lattice, Si-Si bond length, and the buckling distance for silicene without strain are 3.81, 2.24, and 0.42 Å, respectively, which are in good agreement with previous calculations[35, 36]. We also calculate the electronic structure of silicene without strain. The energy dispersion surface clearly shows the Dirac cones at the Brillouin zone (BZ) boundary at the Fermi level (Fig. 1c), which induces the massless Dirac fermion behavior near the Fermi level like graphene. Uniform biaxial and uniaxial strain are applied to silicene to investigate the strain effect on HHG. Two special directions of silicene, the armchair (AC) and zigzag (ZZ) directions are marked in Fig. 1a. In our calculations, strain is defined as , where and are the lattice constants with and without strain, respectively.

Figure 1d shows the HHG spectra from a monolayer silicene without strain with a driving laser intensity of W/cm. Only odd harmonic orders are present, reflecting the centrosymmetric nature of the crystal lattice. Up to the seventh harmonic order are obtained. It is worthy to note that the experimentally measured HHG from monolayer graphene is up to the fifth order[30] with a driving intensity of W/cm and up to the ninth order[29] with a driving intensity of W/cm. Considering the similarity of electronic structures between silicene and graphene, our simulation results are in good agreement with the experimental observations, showing the reliability of our simulations. Besides, both experiments have confirmed the HHG is in the nonperturbative regime[29, 30]. The same characteristic is expected for HHG in the parameter regime of our simulations.

We then investigate the effect of laser polarization orientation on HHG from the 2D material with Dirac cones. Previous studies state that HHG is isotropic in graphene due to rotational symmetry in the band structure of the Dirac cone[29, 30]. Here we consider laser polarization orienting along two specific directions, i.e., the AC and ZZ directions. Simulation results show evidently different HHG spectra, especially for the high-order harmonics. Therefore, even 2D materials with symmetric Dirac cone exhibit anisotropic emission of high harmonics. On the other hand, the anisotropic HHG spectra reflect the symmetries of the BZ of the material. We note that the band structure is isotropic only in the vicinity of the Dirac points. It is significantly anisotropic away from the Dirac points (Fig. 1b). Further increasing anisotropy of the band structure by strain engineering is thus expected to result in enhanced anisotropic HHG emission for different laser polarization. More importantly, modulated band structure by strain engineering should also lead to significant change in the HHG spectra compared to the case where no strain is applied.

The computed HHG spectra of silicene controlled by strain engineering are presented in Fig. 2. For the case of silicene under uniform biaxial strain, the silicene structure maintains the hexagonal symmetry while with a elongate lattice constant. Comparing to the no strain case, although the strain is as small as 6%, the intensity of the fifth harmonics is greatly enhanced, which is more than 6 and 4 time stronger for the AC and ZZ laser polarization configurations, respectively. Meanwhile, the intensity of other harmonic orders is nearly unchanged for both laser polarization configurations. In addition to the intensity enhancement, the full-width at half-maximum (FWHM) of the fifth harmonic peak increases by a factor of about 1.7 and 1.6 for the AC and ZZ laser polarization configurations, respectively.

For the case of silicene under uniaxial strain, the symmetry of silicene is broken, which is expected to bring new features to HHG. In this work, we only consider the effect of tensile strain along the ZZ direction for demonstration. Due to Poisson’s effect, the lattice constant along the AC direction decreases when the lattice is stretched along the ZZ direction. Consequently, the geometric structure and BZ no longer maintain the regular hexagonal symmetry. For the HHG spectra under uniaxial strain of 6%, in contrast to the biaxial strain case under the same strain, not only the fifth harmonic order, but also the other harmonic orders are all significantly enhanced in intensity. The enhancement for the fifth harmonics reaches an order of magnitude, which is approximately 12 and 10 times compared to the no strain case for the AC and ZZ laser configurations respectively, larger than that of the biaxial strain case. Moreover, the spectral width exhibits a strong anisotropic modulation. The FWHM of the fifth harmonic peak is broadened by 1.4 times for the ZZ laser configuration, while nearly unchanged for the AC laser configuration.

To understand the above results, we next study the electronic structure of siliene under different strain profiles, since the in-plane HHG is closely related to the band structure for 2D materials. First we explore the band structure over the whole BZ. In order to get a simple physical picture, we assume the HHG is mainly contributed by the highest valence band and lowest conduction band. Due to the zero bandgap feature, Dirac fermions can easily respond to external laser field. The carriers can give rise to HHG through coupled contribution from nonlinear intraband current, i.e., laser-driven dynamic Bloch oscillation within the non-parabolic band, and interband polarization, i.e., direct electron-hole recombination between the conduction and valence bands. The carrier dynamics and thus the HHG can be strongly influenced by the modulation of band structures. From Fig. 3a-c, we observe that the highest valence band does not change much under various strain profiles comparing to the lowest conduction band counterparts. Therefore we focus on the the lowest conduction band, especially the states near the Fermi level due to their large population of carriers. Without the presence of strain, the states with low energies all locate at the vicinity of the six K points (i.e., the Dirac cones). The area in which Dirac cones keep the rotational symmetry is also the largest among the three strain profiles (Fig. 3d-f). For the biaxial strain case, although the BZ has the same symmetry as the no strain case, the area where the Diac cones have rotational symmetry shrinks. Furthermore, besides the states around the Dirac cones, the states near the point also locate closely above the Fermi level (Fig. 3e). For the uniaxial case, the BZ and band structure have a lower symmetry than those of the other two strain cases. Similar with the biaxial case, new states emerge closely above the Fermi level, but they locate near the S point and show high directivity in -S direction (Fig. 3f). The lower the states in the conduction bands are, the more easily the free carriers can be excited by the laser field from the highest valence band to these lowered states, which will also contribute the HHG process.

To get further insights, we analyse the paths in the BZ which correspond to the MDF dynamics in the reciprocal space. For the sake of better illustration, we choose four representative paths where the states are close to the Fermi level, which are expected to have major contributions to the HHG process. Path I and II correspond to the ZZ laser polarization configuration (Fig. 4a-b), while path III and IV to the AC laser polarization configuration (Fig. 4c-d). Notice the direction in reciprocal and real space are opposite. Under the uniaxial strain, the Dirac points will shift parallel to the K2 direction[36]. However, they only shift very slightly in the considered small strain of 6% , thus we use the same paths for the uniaxial strain case as the other two strain cases for consistency. We first investigate MDF dynamics in the ZZ laser polarization configuration. Considering path I, when there is no strain, most of the states near the Fermi level are from the -band. Besides, there is only a small overlap between the - and -band near the point in CBM (Fig. 4e). With the presence of strain, the overlap between the - and -band becomes larger, and the states near the points shift towards the Fermi level (Fig. 4f-g). The hybridization of - and -band leads to a larger anharmonicity of CBM comparing to the no strain case, which may result in an enhanced contribution to the HHG via the intraband mechanism. In path II, we also observe this hybridization and enhanced anharmonicity of the CBM under strains (Fig. 4h-j). Thus both path I and II may have a contribution to the HHG enhancement in the ZZ laser polarization configurations. We also notice that in path I the states near the points are closer to the Fermi level in the biaxial strain case (Fig. 4f) than those in the uniaxial strain case (Fig. 4g). However, the states in path II near the S1 point in the uniaxial strain case (Fig. 4j) are even closer to the Fermi level than those near the point in the biaxial case in path I (Fig. 4f). Thus the CBM in the uniaxial strain case can have a higher carrier population, which may lead to the enhanced HHG intensity of all harmonic orders in the uniaxial strain case.

We now study the case of AC laser polarization configuration. From the band structures in path III, we observe there is a small overlap between the - and -band in the biaxial strain case (Fig. 4l) while these two bands remain untouched in the no strain (Fig. 4k) and uniaxial strain (Fig. 4m) cases. But for all strain cases most of the low-lying states still locate near the Dirac points (Fig. 4k-m). Therefore this path contribute little to the change of HHG intensity under different strain profiles. In contrast, new states emerge close to the Fermi level under biaxial and uniaxial strains in path IV, which locate near the and S points for the biaxial (Fig. 4o) and uniaxial (Fig. 4p) strains, respectively. The states in the uniaixal strain case (Fig. 4p) are even closer to the Fermi level than those in the biaxial case (Fig. 4o). Moreover, these states in the uniaxial strain case form a rather flat band near the S1 and S2 points (Fig. 4p) and thus induce a high density of states (DOS). The low-lying states near the S points with a high DOS can lead to a high carrier population. Consequently, this may result in a higher HHG intensity in the uniaxial strain case than in the other two strain cases.

The presence of strain changes the lattice structures of silicene, which changes the charge density. Moreover, strain also brings profound modulation to the electronic band structures. Consequently, the carrier density is further varied. These strain effects can lead to modification of the nonlinear refraction index of the material and thus phase changes of the light waves. As a result, spectral broadening can be induced via nonlinear processes such as self-phase modulation and cross phase modulation. It should be noted that the HHG in 2D materials is a complicated process originating from the interplay between the intraband and interband dynamics. To elucidate more details underlying the HHG physical mechanism, it is helpful to develop appropriate theoretical models to assist interpretation of the simulation results, which is beyond the scope of the present work.

## Discussion

In summary, we have investigated tunable high-harmonic emission from a monolayer 2D material with Dirac fermions driven by in-plane electric fields using first-principles TDDFT calculations. We show that HHG in silicene can be actively controlled in terms of intensity and spectral width by modulating the electronic properties of the 2D material via strain engineering. The HHG is sensitive to the electronic band structure changes. The HHG intensity can be enhanced by an order of magnitude under a small uniaxial strain of only 6%. The strain engineering of HHG should be applicable also in other similar 2D materials with MDF, e.g., graphene. Many other approaches, such as using substrate, doping, and electrostatic gating, can be introduced to tune the electronic properties of 2D materials. The study thus opens a new avenue to the development of compact solid-state optoelectronic nano-devices for tailored extreme-ultraviolet and attosecond pulses. It can also provide important insights for the investigation of strong-field and mechanically induced nonlinear response of Dirac carriers in 2D materials. By measuring the time evolution of HHG, it should be possible to obtain novel ultrafast dynamical information of the Dirac fermions. Besides, HHG offers a valuable tool to understand the strain effect on silicene monolayer, which is also of great importance for 2D material growth.

## Methods

Silicene structures are studied by using the semiperiodic supercell model. Silicene structures without strain and under biaxial strain are modelled by the hexagonal primitive cell containing two silicon atoms, while silicene structures under uniaxial strain are modelled by a rectangular supercell containing four silicon atoms. A vacuum of 30 Bohr, which includes 3 Bohr of absorbing regions on each side of the monolayer, is chosen to eliminate the interactions between adjacent silicene sheets and avoid the reflection error in the spectral region of interest. All calculations are performed by the Octopus package. The ground state electronic structure properties and geometric structure relaxation are performed within the density functional theory (DFT) framework in the local density approximation (LDA). Time evolution of the wave functions and time-dependent electronic current are studied by propagating the Kohn-Sham equations in real time and real space within the time-dependent DFT (TDDFT) framework in the adiabatic LDA (ALDA). The real-space spacing is 0.4 Bohr. A Monkhorst-Pack k-point mesh for the BZ sampling is used for silicene without strain, and the sampling is scaled according to the size of the supercells. The fully relativistic Hartwigsen, Goedecker, and Hutter (HGH) pseudopotentials are used in all our calculations.

The laser is described in the velocity gauge. A sin-squared envelope and a carrier-envelope phase of is used for the laser profile. The central laser wavelength is 1600 nm, corresponding to a photon energy of 0.77 eV. The pulse duration is =15 fs. The peak laser intensity in vacuum is taken to be W/cm. The linearly polarized laser field is normally incident onto the silicene sample so that the driving electric field is in the plane of the monolayer.

The HHG spectrum was calculated from the total time-dependent electronic current as:

(1) |

where FT denotes the Fourier transform.

### Data availability.

The data that support the findings of this study are available from the corresponding authors upon request.

The OCTOPUS code is available from http://www.octopus-code.org.

## References

- [1] Ghimire, S. et al. Observation of high-order harmonic generation in a bulk crystal. Nat. Phys. 7, 138–141 (2011).
- [2] Schubert, O. et al. Sub-cycle control of terahertz high-harmonic generation by dynamical Bloch oscillations. Nat. Photon. 8, 119–123 (2014).
- [3] Vampa, G. et al. Linking high harmonics from gases and solids. Nature 522, 462–464 (2015).
- [4] Hohenleutner, M. et al. Realtime observation of interfering crystal electrons in high-harmonic generation. Nature 523, 572–575 (2015).
- [5] Ndabashimiye, G. et al. Solid-state harmonics beyond the atomic limit. Nature 534, 520–523 (2016).
- [6] Kruchinin, S. Yu., Krausz, F. & Yakovlev, V. S. Strong-field phenomena in periodic systems. Rev. Mod. Phys. 90, 021002 (2018).
- [7] Vampa, G. et al. Theoretical analysis of high-harmonic generation in solids. Phys. Rev. Lett. 113, 073901 (2014).
- [8] Otobe, T. High-harmonic generation in -quartz by electron-hole recombination. Phys. Rev. B 94, 235152 (2016).
- [9] Osika, E. N. et al. Wannier-Bloch approach to localization in high-harmonics generation in solids. Phys. Rev. X 7, 021017 (2017).
- [10] You, Y. S., Reis, D. A. & Ghimire, S. Anisotropic high-harmonic generation in bulk crystals. Nat. Phys. 13, 345–349 (2017).
- [11] Tancogne-Dejean, N., Mücke, O. D., Kärtner, F. X. & Rubio, A. Impact of the electronic band structure in high-harmonic generation spectra of solids. Phys. Rev. Lett. 118, 087403 (2017).
- [12] Vampa, G. et al. All-optical reconstruction of crystal band structure. Phys. Rev. Lett. 115, 193603 (2015).
- [13] Luu, T. T. et al. Extreme ultraviolet high-harmonic spectroscopy of solids. Nature 521, 498–502 (2015).
- [14] Luu, T. T. & Wörner, H. J. Measurement of the Berry curvature of solids using high-harmonic spectroscopy. Nat. Commun. 9, 916 (2018).
- [15] Garg, M., Kim, H. Y. & Goulielmakis, E. Ultimate waveform reproducibility of extreme-ultraviolet pulses by high-harmonic generation in quartz. Nat. Photon. 12, 291–296 (2018).
- [16] Garg, M. et al. Multi-petahertz electronic metrology. Nature 538, 359–363 (2016).
- [17] Hammond, T. J. et al. Integrating solids and gases for attosecond pulse generation. Nat. Photon. 11, 594–599 (2017).
- [18] Tancogne-Dejean, N., Mücke, O. D., Kärtner, F. X. & Rubio, A. Ellipticity dependence of high-harmonic generation in solids originating from coupled intraband and interband dynamics. Nat. Commun. 8, 745 (2017).
- [19] Langer, F. et al. Symmetrycontrolled temporal structure of high-harmonic carrier fields from a bulk crystal. Nat. Photon. 11, 227–231 (2017).
- [20] Sivis, M. et al. Tailored semiconductors for high-harmonic optoelectronics. Science 357, 303–206 (2017).
- [21] Han, S. et al. High-harmonic generation by field enhanced femtosecond pulses in metal-sapphire nanostructure. Nat. Commun. 7, 13105 (2016).
- [22] Tancogne-Dejean, N. & Rubio, A. Atomic-like high-harmonic generation from two-dimensional materials. Sci. Adv. 4, eaao5207 (2018).
- [23] Liu, H. et al. High-harmonic generation from an atomically thin semiconductor. Nat. Phys. 13, 262–265 (2017).
- [24] Mikhailov, S. A. Non-linear electromagnetic response of graphene. Europhys. Lett. 79, 27002 (2007).
- [25] Bowlan, P. et al. Ultrafast terahertz response of multilayer graphene in the nonperturbative regime. Phys. Rev. B 89, 041408(R) (2014).
- [26] Al-Naib, I., Sipe, J. E. & Dignam, M. M. High harmonic generation in undoped graphene: interplay of inter- and intraband dynamics. Phys. Rev. B 90, 245423 (2014).
- [27] Al-Naib, I., Sipe, J. E. & Dignam, M. M. Nonperturbative model of harmonic generation in undoped graphene in the terahertz regime. New J. Phys. 17, 113018 (2015).
- [28] Cox, J. D., Marini, A. & García de Abajo, F. J. Plasmon-assisted high-harmonic generation in graphene. Nat. Commun. 8, 14380 (2017).
- [29] Yoshikawa, N., Tamaya, T. & Tanaka, K. High-harmonic generation in graphene enhanced by elliptically polarized light excitation. Science 356, 736–738 (2017).
- [30] Taucer, M. et al. Nonperturbative harmonic generation in graphene from intense midinfrared pulsed light. Phys. Rev. B 96, 195420 (2017).
- [31] Baudisch, M. et al. Ultrafast nonlinear optical response of Dirac fermions in graphene. Nat. Commun. 9, 1018 (2018).
- [32] Runge, E. & Gross, E. K. U. Density-functional theory for time-dependent systems. Phys. Rev. Lett. 52, 997–1000 (1984).
- [33] van Leeuwen, R. Causality and symmetry in time-dependent density-functional theory. Phys. Rev. Lett. 80, 1280–1283 (1998).
- [34] Cahangirov, S. et al. Two- and onedimensional honeycomb structures of silicon and germanium. Phys. Rev. Lett. 102, 236804 (2009).
- [35] Qin, R., Wang, C.-H., Zhu, W. & Zhang, Y. First-principles calculations of mechanical and electronic properties of silicene under strain. AIP Adv. 2, 022159 (2012).
- [36] Qin, R., Zhu, W., Zhang, Y. & Deng, X. Uniaxial strain-induced mechanical and electronic property modulation of silicene. Nanoscale Res. Lett. 9, 521 (2014).

## Acknowledgements

We acknowledge financial support from the National Natural Science Foundation of China (NSFC) (11705185) and the Presidential Fund of China Academy of Engineering Physics (CAEP) (YZJJLX2017002).

## Author contributions

Both authors contributed substantially to all parts of the work.

## Additional information

### Competing interests:

The authors declare no competing financial interests.