## Abstract

The generation of high-order harmonics from atomic and molecular gases enables the production of high-energy photons and ultrashort isolated pulses. Obtaining efficiently similar photon energy from solid-state systems could lead, for instance, to more compact extreme ultraviolet and soft x-ray sources. We demonstrate from ab initio simulations that it is possible to generate high-order harmonics from free-standing monolayer materials, with an energy cutoff similar to that of atomic and molecular gases. In the limit in which electrons are driven by the pump laser perpendicularly to the monolayer, they behave qualitatively the same as the electrons responsible for high-harmonic generation (HHG) in atoms, where their trajectories are described by the widely used semiclassical model, and exhibit real-space trajectories similar to those of the atomic case. Despite the similarities, the first and last steps of the well-established three-step model for atomic HHG are remarkably different in the two-dimensional materials from gases. Moreover, we show that the electron-electron interaction plays an important role in harmonic generation from monolayer materials because of strong local-field effects, which modify how the material is ionized. The recombination of the accelerated electron wave packet is also found to be modified because of the infinite extension of the material in the monolayer plane, thus leading to a more favorable wavelength scaling of the harmonic yield than in atomic HHG. Our results establish a novel and efficient way of generating high-order harmonics based on a solid-state device, with an energy cutoff and a more favorable wavelength scaling of the harmonic yield similar to those of atomic and molecular gases. Two-dimensional materials offer a unique platform where both bulk and atomic HHG can be investigated, depending on the angle of incidence. Devices based on two-dimensional materials can extend the limit of existing sources.

## INTRODUCTION

The possibility of realizing an efficient high-harmonic generation (HHG) from bulk materials has recently been attracting much attention (*1*–*9*). Beyond this very active field of research lies the possibility of extending the present techniques used in the atto-science community for atomic and molecular gases to solid-state systems. For instance, solid-state–based HHG not only could lead to more compact and brighter extreme ultraviolet sources but also enables on-chip ultrafast photonics (*10*), all-optical processors (*11*), and multi-petahertz electronics (*12*). It was reported that, for the same laser intensity, HHG from van der Waals solids can extend beyond the atomic limit for the same species (*6*). However, the maximum intensity that could be applied to a bulk material is intrinsically limited by its damage threshold. Because bulk materials also exhibit a less favorable scaling of the energy cutoff in terms of electric field strength (*1*) or wavelength (*13*), solid-state materials are not yet at the point of replacing the well-established HHG sources based on atomic and molecular gases. Yet, due to their infinite extension, solids can be a good alternative to gases for specific tasks, such as in producing circularly polarized high-order harmonics. HHG from circularly polarized driving field is possible in solids (*7*, *14*), and even the generation of circularly polarized harmonics from a single-color circularly polarized driving field is now possible (*14*).

In the context of finding new and more efficient solid-state materials for HHG, quasi–two-dimensional (2D) or monolayer materials have recently attracted considerable attention. Graphene (*9*, *15*, *16*) and molybdenum disulphide (MoS_{2}) (*9*, *17*) have been investigated both theoretically and experimentally. In particular, it has been shown that an in-plane driving electric field generates harmonics both along and perpendicular to the laser electric field polarization direction (*17*), and that harmonics can even be obtained from elliptically polarized pumping lasers (*9*). In a previous work, we have demonstrated that monolayer and few-layer materials are predicted to be efficient materials for HHG because of their strong inhomogeneity of the electron-nuclei potential (*13*). This improved efficiency, with respect to their bulk counterpart, has already been verified experimentally for MoS_{2} (*17*). Also, plasmon-enhanced harmonic generation has recently been proposed from finite size nanostructures such as graphene nanotriangles (*18*).

Thus far, the investigation of HHG from monolayer materials has been limited to in-plane driving electric field, and the possible harmonic generation from an out-of-plane driving electric field has remained unexplored. Here, we theoretically investigate the HHG from a free-standing monolayer of hexagonal boron nitride (hBN) driven by an intense (*I*_{0} = 10^{14} W cm^{–2}) and short (15 fs) laser pulse polarized purely out of plane. This corresponds to the limit of grazing incidence. We estimate from our simulations the minimum angle of incidence required to avoid damage by the in-plane component of the field in the more realistic grazing incidence geometry, and found a minimum angle of incidence of 77°, keeping the out-of-plane intensity at *I*_{0} = 10^{14} W · cm^{–2}. Our findings are general and should apply to all quasi-2D materials, such as transition metal dichalcogenides or graphene, and to both transmission and reflection HHG measurements. It is important to note that we do not take any assumption as the single-active electron approximation or the independent-particle approximation, which we will show below is not valid for monolayer materials. We included a full description of electron-electron and electron-ion interactions in our time-dependent density functional ab initio simulations. We focus the discussion below to the similarities and discrepancies between the atomic case, as predicted by the widely used three-step model (*19*, *20*), and the results of our first-principles simulations for the quasi-2D material.

## RESULTS

### Atomic-like harmonic emission

We consider a monolayer of hBN as a prototypical material to illustrate the relevance of the electron-electron correlations, electronic band structure, and spatial anisotropy effects in HHG. This material has a calculated work function *E*_{w} of 6.02 eV within the local density approximation (LDA), obtained from the highest occupied state, located at the *K* point in the Brillouin zone. This work function is the minimum energy that must be transferred to an electron to promote it to the continuum. This energy plays the equivalent role of the ionization energy in the three-step model of HHG in gases (*19*, *20*). Because of the high energy required to promote electrons to the continuum, we found that monolayer hBN can be exposed to a laser field as intense as *I*_{0} = 10^{14} W · cm^{–2} in vacuum, for an out-of-plane electric field, without breaking the material. As will be discussed later, this is in sharp contrast to the case of an in-plane driving electric field.

If we assume the validity of the three-step model for monolayer materials driven by an out-of-plane electric field, we would expect to obtain a cutoff energy of *E*_{c} ≈ *E*_{w} + 3.17 *U*_{p} for the HHG spectrum (*21*), where *U*_{p} is the ponderomotive energy of the electron in the applied laser field. Taking a laser intensity in vacuum of *I*_{0} = 10^{14} W · cm^{–2} and a pumping wavelength of λ = 1600 nm, the cutoff energy should be about *E*_{c} ≈ 81.8 eV. We computed the HHG spectra of monolayer hBN under these conditions in the framework of time-dependent density functional theory (TDDFT) (*13*, *14*). Numerical details are provided in Methods.

As shown in Fig. 1A, harmonics up to the order 110 are obtained, which correspond to an energy cutoff of 85 eV. It is noteworthy that, in comparison, the energy cutoff for the in-plane polarized driving laser field is much lower, extending up to the harmonic order 25 (see discussion below and also fig. S2). The extension of the harmonic plateau, and thus of the energy cutoff, is in very good agreement with the one predicted by the three-step model, indicating an atomic-like behavior of the free-standing monolayer hBN when driven by an out-of-plane electric field. This is not the case for the in-plane polarization (see, for instance, fig. S2) for which the energy cutoff behaves as in bulk crystals (*1*, *13*). To confirm the semiclassical behavior of the electrons, we performed a time-frequency analysis of the harmonic emission (see Fig. 1B). Like in gases, we observe well-defined trajectories that give rise to the different harmonics.

Our time-frequency analysis also reveals that the harmonic emission extends in energy up to half-cycle energy cutoffs defined by , where *E*_{hcm} is the half-cycle maximum field strength, as expected from atomic-like harmonic generation from a few-cycle laser pulse (*22*).

We note that even if the intensity considered here is quite high, it was recently demonstrated experimentally (*23*) that fused silica is not irreversibly damaged by a short driving field of the same intensity as considered here (10^{14} W · cm^{–2}), with a similar wavelength (1700 nm).

In Fig. 2, a color plot of the evolution of the induced electronic density, averaged in the plane of the monolayer, is shown. The wave packets evolve in bundles and follow well visible real-space trajectories. The trajectories that return to the monolayer are those that result in harmonic generation, as in the case of gases.

Our results show that we can obtain high-order harmonic emission from monolayer materials, as compared to bulk crystals, based on the same mechanism as the one responsible for atomic and molecular HHG. The electrons follow well-defined trajectories that can be matched by the semiclassical model of electron dynamics in the presence of an electric field. This similarity with the atomic and molecular case opens the path toward developing robust, compact, and efficient experimental devices to generate controlled attosecond pulses from low-dimensional solids.

### Anisotropy effects

Most of the simulations of HHG in atomic and molecular gases rely on the validity of the single-active electron approximation. However, in polyatomic molecules, this approximation fails (*24*, *25*). In the case of bulk materials, we have shown that a related approximation, namely, the independent-particle approximation, performs extremely well for the HHG of not strongly correlated materials, such as bulk silicon (*13*) or MgO (*14*). In the case of monolayer materials, it is natural to ask whether this approximation is still valid or not.

To assess the validity of the independent-particle approximation in the case of a free-standing monolayer driven by an electric field polarized along the out-of-plane direction, we followed the approach used by Tancogne-Dejean *et al*. (*13*). More precisely, we computed the time-dependent electronic current by performing either the full time-dependent electronic evolution of the self-consistent Kohn-Sham Hamiltonian or the time evolution of the electrons in a frozen potential corresponding to the ground-state Kohn-Sham potential. The results are shown in Fig. 3, and they reveal that assuming the independent-particle approximation leads to a severe overestimation of the electric field acting on electrons and therefore predicts incorrect harmonic spectra (see Fig. 3B). The line shape of the harmonic spectra is incorrect in the low-energy region, and the overall intensity is incorrect by a factor of 10. This can be understood easily in terms of the anisotropic screening of electrons in and out of plane [sometimes referred to as “local-field effects” in solid-state physics (*26*)]. As shown in Fig. 4, the induced electronic density of the monolayer becomes polarized during the interaction with the external electric field, like a dipole aligned along the direction of the applied laser field, that is, perpendicular to the monolayer’s plane. These fluctuations of the electronic density locally induce an electric field that strongly screens the external applied laser field. Therefore, the total electric field felt by electrons, which is responsible for harmonic generation, is strongly reduced. Because the electric field strength dictates the ionization rate, this local-field effects tend to considerably reduce the probability of ionization and thus of promoting electrons in the continuum. Therefore, harmonic generation is reduced.

The results indicate that the first step (ionization) of the commonly accepted three-step model of HHG in gases is strongly affected by the electron-electron interaction (that is, dynamical electronic screening) in the case of monolayer materials. The independent-particle approximation, which is valid for solids such as bulk silicon, is no longer valid for monolayer materials. Because of the local nature of the induced electric field, the energy cutoff is correctly predicted within the independent-particle approximation. However, this approximation leads to qualitatively incorrect (i) harmonic spectral shape at low energy, (ii) harmonic yield on the entire spectrum (as the ionization step is modified by the local induced field), and (iii) time profile of the electronic current. This conclusion applies to all monolayer materials and might also have important consequences for the HHG of heterostructures made by stacking few-layer materials. Note that similar effects will be present at the surface of a bulk material, and it would be important to carefully investigate them to distinguish between bulk and surface contributions to HHG from solids.

### Wavelength scaling of the HHG yield

Because of the λ^{2} dependence of the ponderomotive energy (*U*_{p} ∝ *I*_{0}λ^{2}), it is always tempting to increase the driver wavelength to obtain higher-order harmonics. However, in the case of single-atom HHG, it was found that the harmonic yield decreases drastically with increasing wavelength (*24*, *27*). The decrease of the harmonic yield with increasing wavelength is partially understood in terms of a semiclassical picture based on a three-step model in which the electron is first promoted to the continuum, then evolves as a free particle driven by an electric field, and finally recombines with the parent ion. The quantum nature of the process is taken into account by considering that the electron is not a point-charged particle but is described by a wave packet that is initially localized on the parent ion. During the evolution of the electron as a free particle in vacuum, this wave packet spreads in time and space. The efficiency of the harmonic emission in this semiclassical picture is dictated by the overlap of the wave packets of the electron and parent ion (*28*). Because of the time spreading of the electron wave packet, a λ^{− 3} scaling is expected (*21*, *27*). However, a theoretical investigation (*29*) of the evolution of the harmonic yield with respect to the driver wavelength found, from the solution of the time-dependent Schrödinger equation, that the integrated harmonic yield of the isolated atom scales between λ^{−5} and λ^{−6} (*29*). Later, it has been shown experimentally, for wavelengths ranging from 800 to 1850 nm, that the single-atom HHG emission scales as λ^{−6.3} and λ^{−6.5} for Xe and Kr, respectively (*30*). The wavelength dependence of the yield of atomic HHG has been intensively investigated [see, for instance, previous studies (*27*, *29*–*33*)] and is, by now, well understood in terms of quantum diffusion of the electron wave packet, which corresponds to a time and space spreading of the electron wave packet during its free evolution in vacuum. Because of this spreading, which increases with driver wavelength, the recombination of the electron wave packet with the parent ion becomes less and less efficient, reducing the HHG yield.

However, in the case of a quasi-2D material, the quantum diffusion of the wave packet should play a minor role in the electronic recombination; hence, high harmonics could be generated efficiently. To asses this issue, we performed HHG calculations for monolayer hBN driven by different wavelengths (from 800 to 1600 nm), keeping the intensity and pulse duration fixed. The corresponding calculated HHG spectra are shown in fig. S1. The comparison with the cutoff energy predicted by the three-step model (vertical lines in fig. S1) confirms that the harmonic cutoff scales as λ^{2} for monolayer materials when driven in the out-of-plane direction.

We found that the integrated harmonic yield from 15 to 35 eV (see Fig. 5) exhibits a more favorable scaling of λ^{−4.73} than the experimental or theoretical ones reported before for a single atom (*27*, *29*–*33*) or molecules (*34*) driven by similar wavelengths. Phase matching plays a fundamental role in HHG from gases. One important consequence of phase matching is that the wavelength scaling observed experimentally deviates from the one predicted by single-atom HHG simulations. However, it was recently shown theoretically that when HHG from gases occurs only at the nanometer scale, using plasmonic field enhancement, most of the sources of phase mismatch do not contribute (*35*). Because of the 2D nature of a free-standing monolayer, the interaction length between light and matter is of subnanometer size even at grazing incidence. We therefore expect that the harmonic emission from monolayer materials will not suffer from phase-matching and propagating effects that are inherent to gases. The present findings demonstrate unambiguously that monolayer materials have a better wavelength scaling than single atoms and consequently gases.

### Grazing incidence

Thus far, we only considered a perfectly out-of-plane driving electric field. However, experimentally, the electric field should be p-polarized and incident with a certain angle. This could drastically limit the possibility of using free-standing monolayer materials, as discussed above. Because of the bulk nature of the material, any monolayer is very sensible to the strength of the in-plane component of the driving electric field, and a high in-plane electric field could lead to an irreversible damage of the material, thus prohibiting its use as a continuous HHG source. It is therefore important to quantify what is the effect of the in-plane component of the electric field in the setting of an experiment performed at grazing incidence.

Because, to the best of our knowledge, there are no available experimental values for the damage threshold of monolayer hBN, we estimated them here from our simulations by assuming that the damage threshold is reached when contributions to the power spectra of the electronic current from the harmonics become more important than those of the fundamental. From these criteria, we estimated that an in-plane electric field up to 0.47 V Å^{−1} can be used for our excitation conditions (λ = 1600nm, pulse duration of 15-fs FWHM), corresponding to an intensity in matter of 7.95 × 10^{12} W · cm^{–2}, taking the experimental in-plane refractive index of *n* = 2.65 (*36*) of bulk hBN. The calculated spectra are shown in fig. S2. We note that our estimated damage threshold value is in good agreement with the previously reported experimental value of the single-shot damage threshold of graphene [3 × 10^{12} W · cm^{–2} in vacuum (*37*)], showing the validity of our estimation method. We also note that the use of the static dielectric constant and refractive index is well justified because our excitation wavelength is very small compared to the band gap of the material.

Because of the large difference between the critical in-plane and out-of-plane fields mentioned above, only grazing incidence experimental excitation can be envisioned. If we want to drive electrons out of plane with an electric field strength of 2.74 V Å^{−1} (corresponding to *I*_{0} = 10^{14}W · cm^{–2} in vacuum), then we should consider at least an incidence angle of 77.7° (or 80.5° if we only want a lower in-plane electric field strength of 0.3 V Å^{−1} in matter), in which we again use the in-plane bulk static dielectric constant (ϵ_{||,0} = 7.04) (*36*) to evaluate the transmission of the free-standing hBN monolayer.

Therefore, our ab initio estimation indicates that it could be possible to expose a free-standing monolayer of hBN to an out-of-plane intense electric field, still being below the in-plane damage threshold of the material, at an experimentally realistic grazing incidence angle. Note that similar conclusions should be valid for other quasi-2D materials.

## DISCUSSION

In conclusion, we investigated the high-harmonic emission from a free-standing monolayer, taking the hBN as a prototypical case. We demonstrated atomic-like HHG when electrons are driven out of plane by the electric field of the pump laser, together with the normal bulk HHG response (*13*, *14*) for electron driven in the plane of the monolayer. By replacing the ionization potential with the work function of the material, we found that cutoff energy predicted by the three-step model of atomic HHG agrees well with the results of our simulations. The time-frequency analysis of the harmonic emission revealed that the HHG from the monolayer originates from clear out-of-plane trajectories, similar to HHG from gases. This allows a direct extension of present experimental techniques to produce isolated attosecond pulses from monolayer materials. We found that the electron-electron interaction plays a major role in the HHG from quasi-2D materials, due to strong anisotropic screening and local-field effects. We further investigated the wavelength scaling and found that it is more favorable than in the case of single-atom HHG, making monolayer materials an attractive alternative to generate high-energy photons. Finally, we estimated the in-plane damage threshold and showed the possibility of performing a grazing incidence experiment, using the intense out-of-plane electric field, while staying below the in-plane damage threshold at an experimentally feasible incidence angle. Further investigations should address excitonic effects and the role of the electron-phonon coupling, defects, and substrates in the HHG from monolayer materials. Another important question is the effect of the confinement field from the substrate on the HHG spectra in transmission and reflection of monolayer materials. Because of their unique nature, 2D materials offer the opportunity to investigate both bulk and atomic HHG, using the angle of incidence as the control parameter between the two limiting cases of in-plane and out-of-plane electric field. This work opens the door to a new type of solid-state HHG sources, which could become a viable alternative to the current gas-based HHG sources.

## METHODS

The time evolution of the wave functions and the evaluation of the time-dependent electronic current were computed by propagating the Kohn-Sham equations in real space and real time within TDDFT, as implemented in the Octopus code (*38*), in the adiabatic LDA (*39*) and with semiperiodic boundary conditions. All calculations were performed using fully relativistic Hartwigsen, Goedecker, and Hutter (HGH) pseudopotentials (*40*).

We used a simulation box of 360 bohr along the nonperiodic dimension, including 40 bohr of absorbing regions on each side of the monolayer. The absorbing boundaries are treated using the complex absorbing potential method, and the cap height η is taken as η = −1 atomic units (a.u.) to avoid the reflection error in the spectral region of interest (*41*). The real-space cell was sampled with a grid spacing of 0.4 bohr. We found that sampling the 2D Brillouin zone with a 24 × 24 *k*-point grid for purely out-of-plane laser field simulations yielded converged results, whereas a 42 × 42 *k*-point grid was required in the case of grazing incidence to reach the convergence. The boron nitride bond length is taken here as 1.445 Å. The laser was described in the velocity gauge (that is, by its vector potential) to treat both out-of-plane and in-plane perturbations, and we used a sin-square envelope. In all our calculations, we used a carrier-envelope phase of φ = 0.

The HHG spectrum was directly obtained from the total time-dependent electronic current **j**(**r**, *t*) as(1)where FT denotes the Fourier transform.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/2/eaao5207/DC1

fig. S1. HHG spectra of monolayer hBN for various wavelengths.

fig. S2. In-plane HHG spectra of monolayer hBN.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgment:**We would like to acknowledge O. D. Mücke and M. Altarelli for very interesting and fruitful discussions.

**Funding:**We acknowledge financial support from the European Research Council (ERC-2015-AdG-694097).

**Author contributions:**N.T.-D. and A.R. conceived and designed the project. N.T.-D. carried out the code implementation and the numerical calculations. N.T.-D. and A.R. participated in the discussion of the results and contributed to the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The data that support the findings of this study are available from the corresponding authors upon request and will be deposited on the NOMAD repository. The Octopus code is available from www.octopus-code.org.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).