High-Harmonic generation from spin-polarised defects in solids

Generation of high-order harmonics in gases enabled to probe the attosecond electron dynamics in atoms and molecules with unprecedented resolution. Extending the techniques developed originally for atomic and molecular gases to solid state materials requires a fundamental understanding of the physics that has been partially addressed theoretically. Here we employ time-dependent density-functional theory to investigate how the electron dynamics resulting in high-harmonic emission in monolayer hexagonal boron nitride is affected by the presence of vacancies. We show how these realistic spin-polarised defects modify the harmonic emission, and demonstrate that important differences exist between harmonics from a pristine solid and a defected-solid. In particular, we found that the different spin channels are affected differently because of the presence of the spin-polarized point defect, and that localisation of the wavefunction, the geometry of the defect and the electron-electron interaction are all important ingredients to describe high-harmonic generation in defected-solids. We show that different vacancies lead to qualitatively different effects, thus opening the door to the high-harmonic imaging of spin-polarised defects in solids.

Due to the growth processes, defects are inevitable in real solids [24]. Defects in materials can appear in the form of vacancies, impurities, interstitials (all of these can be neutral as well as charged), dislocations, etc. Defect-induced microscopic modifications in a material significantly affect on its macroscopic properties [25]. Electronic, optical, vibrational, structural and diffusion properties of solids with defects have been thoroughly reviewed over the past century [26][27][28][29][30][31][32]. Defect engineering is used to achieve desirable characteristics for materials, e.g., doping has revolutionised the field of electronics [33]. Defects can also be highly controlled, and it is possible to create isolated defects such as nitrogen-vacancy defects in diamond [34,35] or single photon emitters in two-dimensional materials [36].
The influence of defects in solids is not well explored in strong-field physics. In this work, we aim to address the following questions: Is it possible to observe defect-specific fingerprints in strong-field driven electron dynamics? Or does the electron-electron interaction play a different role for defects and bulk materials? Moreover, some defects are also spin polarised in nature. This raises the question if it is possible to control the electron dynamics for different spin channels independently in a non-magnetic host material. As we will show below, HHG is a unique probe that helps us to explore these interesting questions.
We aim to model strong-field driven electron dynamics in a defected-solid with the least approximations. In this work, we employ ab initio time-dependent density functional theory (TDDFT) to simulate the strong-field driven electron dynamics in defected-solid [37]. This allows us to come up with theoretical predictions relevant for real experiments on defected-2 solid without empirical input. There are various theoretical predictions for the HHG from doped semiconductors by using simple one-dimensional model Hamiltonians [38][39][40]. In Ref. [38], by using TDDFT, it is predicted that there is an enhancement of several orders of magnitude in the efficiency of HHG in a donor-doped semiconductor. Using an independent particle model, Huang et al. contrastingly found that the efficiency of the second plateau from the doped semiconductor is enhanced [39]. Similar calculations within a tight-binding Anderson model indicates that the disorder may lead to well-resolved peaks in HHG [41].
The conceptual idea for the tomographic imaging of shallow impurities in solids, within a one-dimensional hydrogenic model, has been developed by Corkum and co-workers [42].
Even though these pioneering works have shown that defects can influence HHG, many points remain elusive. So far, only model systems in one dimension have been considered, and no investigation of realistic defects (through geometry optimisation and relaxation of atomic forces) has been carried out. Beyond the structural aspect, several other essential aspects need to be investigated in order to obtain a better understanding of HHG in defectedsolids such as the importance of electron-electron interaction (that goes beyond single-active electron and independent-particle approximations), the role of the electron's spin, the effect of the symmetry breaking due to the defects, etc. Our present work aims to shed some light on some of these crucial questions. For that, we need to go beyond the one-dimensional model Hamiltonians.
In order to investigate how the presence of defects modifies HHG in periodic materials, we need to select some systems of interest. Nowadays, two-dimensional (2D) materials are at the centre of tremendous research activities as they reveal different electronic and optical properties compared to the bulk solids. Monolayer 2D materials such as transitionmetal dichalcogenides [14,43], graphene [44,45] and hexagonal boron nitride (h-BN) [46][47][48] among others have been used to generate strong-field driven high-order harmonics. Many studies have examined HHG in h-BN when the polarisation of the laser pulse is either inplane [47,48] or out-of-plane [46] of the material. Using an out-of-plane driving laser pulse, Tancogne-Dejean et al. have shown that atomic-like harmonics can be generated from h-BN [46]. Also, h-BN is used to explore the competition between atomic-like and bulk-like characteristics of HHG [46,47]. In MoS 2 , it has been demonstrated experimentally that the generation of high-order harmonics is more efficient in a monolayer in comparison to its bulk counterpart [43]. Moreover, HHG from graphene exhibits an unusual dependence on the laser ellipticity [44]. Light-driven control over the valley pseudospin in WSe 2 is demonstrated by Langer et al. [14]. These works have shed light on the fact that 2D materials are promising for studying light-driven electron dynamics and for more technological applications in petahertz electronics [49] and valleytronics [50].
Different kinds of defects in h-BN can be classified as vacancy (mono-vacancies to cluster of vacancies), antisite, and impurities. In particular, defects like monovacancies of boron and nitrogen atoms in h-BN are among the most commonly observed defects. Recently, Zettili and co-workers have shown the possibility of engineering a cluster of vacancies and characterising them using ultra-high-resolution transmission electron microscopy [56,61,63].
Signatures of defects in h-BN are identified by analysing cathodoluminescence and photoluminescence spectra [51,64]. It is shown that the emission band around 4 eV originates from the transitions involving deep defect levels [65]. Ultra-bright single-photon emission from a single layer of h-BN with nitrogen vacancy is achieved experimentally, for example in large-scale nano-photonics and quantum information processing [36,53]. In Refs. [54,[65][66][67]69], different kinds of defects in h-BN are modeled, and their effects on the electronic and optical properties are thoroughly investigated.
In the present work, using the well-established supercell approach to model the defects in h-BN, we analyse their influence on HHG. Due to the partial ionic nature of its bonds, h-BN is a wide band-gap semiconductor with an experimental band-gap of 6 eV. This makes h-BN an interesting candidate for generating high-order harmonics without damaging the material. Moreover, 2D material like h-BN enable us to easily visualise the induced electron density and localised defect states easily. The presence of boron or nitrogen vacancies in h-BN acts as a spin-polarised defect. These factors make h-BN an ideal candidate for exploring spin-resolved HHG in non-magnetic defected solids. Below, we will discuss HHG in h-BN with and without spin-polarised monovacancies.

Computational approach
h-BN has a two-atoms primitive cell. We model the vacancy in a 5×5 (7×7) supercell with 50 atoms (98 atoms) (See methods section for the details of structural optimisation) following the methods in Ref [54]. The size of the 5×5 supercell is large enough to separate the nearest defects with a distance greater than 12Å. This minimises the interaction between the nearby defects [54,69], and the defect wavefunctions are found to be well-localised within the supercell [66]. In the present work, we are not considering more than one point defect.
Both of our vacancy configurations (single boron as well as single nitrogen vacancy) have a total magnetic moment of +1µ B , which is consistent with the earlier reported results [54,66,69]. For the low defect-concentration considered here (2% and ∼ 1%), the strength of the total magnetic moment is independent of the defect concentrations [69]. We have found that the 7×7 and 5×5 supercells converged to similar ground-states [54].
In all the present calculations, we consider a laser pulse of 15-femtosecond duration at fullwidth half-maximum with sine-squared envelope and a peak intensity of 13.25 TW/cm 2 in matter (for an experimental in-plane optical index n of ∼ 2.65 [70]). The carrier wavelength of the pulse is 1600 nm, which corresponds to a photon energy of 0.

High-harmonic generation in hexagonal boron nitride with a boron vacancy
We start our analysis by comparing the HHG from pristine h-BN and from h-BN with a boron vacancy. Removal of a boron atom makes the system spin-polarised [54,69]. The high-harmonic spectrum of h-BN with a boron vacancy (V B ) and its comparison with the spectrum of pristine h-BN is presented in Fig. 1a. The spectrum of V B is different from the pristine h-BN as evident from the figure. There are two distinct differences: First, the below band-gap harmonics corresponding to V B are significantly enhanced (see Fig. 1c). Second, the harmonics, closer to the energy cutoff, have a much lower yield for V B in comparison to the pristine material.
An interesting aspect for the spin-polarised defects in a non-magnetic host is that the effect of the defect states, compared to bulk states, can be identified clearly by examining the spin- Unlike the pristine h-BN, spin-up and spin-down electrons in V B see a different bandstructure near the band-gap, since the spin-resolved in-gap states are different. Therefore, spin-up and spin-down electrons evolve differently in the presence of the laser pulse. It means that interband transitions and ionisation involving the defect states will contribute differently to the spectrum. Hence, the spectral enhancement of the third harmonic can be understood as follows: As visible from the spin-down band-structure, there is an additional defect state near the valence band, which allows spin-down electrons to be ionised or to recombine through multiple channels and contribute more to the third harmonic (see Fig. 1).
It is evident from the band structure that additional defect-states effectively reduce the minimal band gap needed to reach the conduction bands. Due to the relaxation of atoms next to the vacancy, the pristine bands are also slightly modified. However, this modification is found to be negligible compared to the photon energy of the laser and is not further discussed. In order to understand how the presence of the defect states influences the interband-tunneling, the number of excited electrons during the laser pulse is evaluated (see  Fig. 3a). In contrast to this, the pulse is not able to promote a significant portion of the valence electrons to the conduction bands permanently in the case of pristine h-BN as the band-gap is significantly large (see the black curve in Fig. 3a). More precisely, for the 5 × 5 supercell, we found that 1.6 electrons are ionised, compared to 0.25 in the case of the same cell in the pristine material.
To achieve a better understanding of the possible ionisation mechanisms, we also consider the induced electron density (n ind ) at two different times near the peak of the vector potential show that electrons remain in the defect states even after the laser pulse is over (see Figs. 3f and 3g). Considering that the three defect states are originally unoccupied, we cannot conclude that more electrons are ionised to the conduction bands. It is most likely that in-gap defect states are filled during the laser excitation as 1.6 electrons are excited and V B is a triple acceptor.
Overall, these results show that the electron dynamics in acceptor-doped solids imply a net transfer of population to the originally unoccupied gap states, but for a wide bandgap host no more electrons are promoted to conduction bands. This explains why only the loworder third harmonic is directly affected by the presence of spin-polarised defect states (as evidenced by the spin dependence of the spectrum). The low density of these defect states means that fewer photons are absorbed. We note that the irreversible population change, assisted by the defect states, ultimately implies that more energy is absorbed by the defected solid, which leads to a lower damage threshold in comparison to the pristine. However, the intensity considered here is low enough to see such an effect.

Effect of electron-electron interaction
We have established so far that the increase of the low-order harmonic yield is compatible with the presence of the defect states in the band gap of h-BN. This is an independentparticle vision, in which we used the ground-state band-structure of V B to explain the observed effect on the HHG spectrum. We now turn our attention to the higher energy harmonics, for which the harmonic yield is decreased. This seems not to be compatible with a simple vision in terms of the single-particle band structure, especially in view of the fact that more electrons are excited by the laser pulse, as shown in Fig. 3a.
To understand this, let us investigate the effect of the electron-electron interaction on the electron dynamics in V B and pristine h-BN. Within the dipole approximation, the HHG spectrum can be expressed as [37,71] Here, n ind (r, t) is the induced electron density, v 0 is the electron-nuclei interaction potential, N e is the total number of electrons, and E is the applied electric field. n ind (r, t) is the difference of the total time-dependent electron density n(r, t) and the ground-state electron density n 0 (r), i.e., n ind (r, t) = n(r, t) − n 0 (r). Also, n(r, t) is decomposed in spin-polarised 8 fashion as n(r, t) = n ↑ (r, t) + n ↓ (r, t). If one analyses this expression, it is straightforward to understand how the introduction of vacancies can change the harmonic spectrum, through a change in the local potential structure near the defect. This results in the change in gradient of the electron-nuclei interaction potential v 0 , which is independent of electrons interacting with each other or not. Apart from this explicit source of change, it is clear that the dynamics of the induced density, evolving from a different ground-state also leads to the modifications in the harmonic spectra. The fact that the ground states are different, and hence n 0 (r) is different, does not affect the harmonic spectrum because of the absence of time dependence in both n 0 and v 0 . The possible difference between the HHG spectra can, therefore, be understood in terms of independent particle effects (originating from the structural change of the nuclear potential v 0 ) and interaction effects from the induced density n ind .
In order to disentangle these two sources of differences between pristine and defect h-BN, we simulate the harmonic spectra within an independent particle (IP) model by freezing the Hartree and exchange-correlation potentials to the ground-state value. The harmonic spectra for pristine h-BN and V B within TDDFT and IP are compared in Fig. 4. In the case of pristine h-BN, the HHG-spectra obtained by the IP model and TDDFT are similar.
Hence, there is no significant many-body effect in HHG from pristine h-BN with an in-plane laser polarisation, at least as described by the PBE functional used here. A similar finding has been reported for Si [37] and MgO [72], within the local-density approximation. Only the high-order harmonics display small differences and there the electron-electron interaction reduces the harmonic yield, as found in Ref. [38]. In contrast, the HHG-spectra obtained by IP and TDDFT are significantly different for V B as reflected in Fig. 4b. This indicates that the electron-electron interaction is essential for HHG in defected solids. The total currents corresponding to the harmonic spectra for the pristine h-BN and V B are shown in Supplementary Fig. 2.
For the case of V B , the induced electron density is displayed in Fig. 3. This helps us to understand how the spatial structure of the defect wavefunctions influences the spectrum.
Similarly, Fig. 3c indicates that the spatial density oscillations are responsible for HHG in pristine h-BN. It is clear that the induced density is different in the two systems. The substantial difference in the harmonic spectra of V B obtained by two approaches, TDDFT and IP, is due to the so-called local field effects. This is explored in detail in the following paragraphs. It will be shown that this difference is responsible for the decrease of the harmonic yield for V B .
In the presence of an external electric-field, the localised induced-charge acts as an oscillating dipole near the vacancy. The dipole induces a local electric field, which screens the effect of the external electric field. This is usually referred to as local field effects. The same mechanism is responsible for the appearance of a depolarisation field at the surface of a material driven by an out-of-plane electric field [46,47]. It is important to stress that this induced dipole is expected to play a significant role here, due to the σ-character of the vacancy wavefunctions. The induced electric field is directly related to the electron-electron interaction term as clearly shown in Ref. [47]. As shown in Fig. 4b, the harmonic yield at higher energies is increased if we neglect local field effects, i.e., we treat electrons at the IP level. Moreover, the electron-electron interaction also affects the third harmonic, as shown in Fig. 4, but lead to an increase of the harmonic yield. We, therefore, attribute this effect not to local field effects but to correlation effects. It is important to stress that in Maxwell equations, the source term of the induced electric field is the induced density (summed over spin). Therefore, the induced electric field acts equally on the HHG from both spin channels, which is why the decrease of the yield occurs equally on both spin channels, see Fig. 1b.
We conclude that the modifications in the HHG spectrum due to a point defect originate from a complex interplay of two important factors: the electronic transitions including the in-gap defect states, and the electron-electron interaction. This indicates that HHG in defected-solids cannot be fully addressed via IP approximation as this can lead to wrong qualitative predictions in certain cases.
Finally, we discuss the dependence of the defect concentration on the HHG spectrum by computing the HHG for a 7×7 supercell with a boron vacancy, which corresponds to a ∼1% doping concentration. The harmonic spectrum is presented in Supplementary Fig. 4.
The third harmonic enhancement persists with comparable intensity even where there is a lower defect concentration, whereas the higher energy region of the harmonic spectrum matches well with the pristine spectrum. This is consistent with the observation that the higher energy spectrum is dominated by the bulk bands. This indicates that some of the present findings depend on the defect concentration and may not be observed below a certain concentration threshold.

High-harmonic generation in hexagonal boron nitride with a nitrogen vacancy
Until now, we have discussed HHG in h-BN with a boron vacancy. Now we will explore HHG in h-BN with a nitrogen vacancy. Similar to the V B case, h-BN with a nitrogen vacancy (V N ) is also spin polarised. Fig. 5a presents the high-order harmonic spectrum of V N and its The weaker local field effects in V N are fully consistent with the HHG spectra of V N , obtained through TDDFT and IP models (see Supplementary Fig. 3). Here, the third harmonic enhancement is well captured within the IP model, but the increase in the yield between 30 to 35 eV is found to originate from the electron-electron interaction (see Supplementary Fig. 3). Finally, in the HHG spectrum of V N , the defect state plays a role even at higher orders (see the spin-resolved spectra in Supplementary Fig. 5), though these effects are feeble.
In essence, the total HHG spectra corresponding to V B and V N include the contributions from the energy-bands of pristine h-BN as well as from the transitions including gap states.
Moreover, electron-electron interaction plays significant and different roles in both cases.
The harmonic spectrum of the defected solid preserves some piece of information of the pristine structure in the higher energy regime along with the characteristic signatures of the 12 defect in the near band-gap regime, or close to the cutoff energy for V N .

DISCUSSION
In summary, we have investigated the role of vacancy-defect in solid-state HHG. For this purpose, we considered h-BN with a boron or a nitrogen atom vacancy. In simple terms, one may assume that h-BN with a boron atom vacancy or with a nitrogen atom vacancy would exhibit similar HHG spectra since a single atom from h-BN has been removed. However, this is not the case as boron and nitrogen vacancies lead to qualitatively different electronic structures, and this is apparent from their corresponding gap states. It has been found that once an atom is removed from h-BN, either boron or nitrogen, the system becomes spin-polarised with a non-zero magnetic moment near the vacancy. As a consequence, the defect-induced gap states are found to be different for each spin channel and for each vacancy, which we found to be strongly reflected in the low-order harmonics. These contributions are strongly spin-dependent, according to the ordering and occupancy of the defect states.
Altogether, the role of the defect states can be understood by analysing the spin-polarised spectra, and the findings are in accordance with the spin polarised band-structure. This establishes one aspect of the role of defect states in strong-field dynamics in solids.
In addition, the vacancy wavefunctions of V B and V N show σ-and π-characters respectively, which lead to different qualitative changes in the harmonic spectra of vacancies, due to the local-field effects and electron correlations. These different behaviours are caused by the creation or absence of an induced dipole, which may counteract the driving electric field, and directly depends on the spatial shape of the defect state wavefunctions. Moreover, the electron-electron interaction also manifests itself in the decrease of the harmonic yield close to the energy cutoff in the V B case, whereas this effect is completely absent in the V N or pristine cases. This implies that the nature of the vacancies in V B and V N is entirely different, as reflected in their HHG spectra, even at a defect concentration as low as 2%. The HHG spectrum of V N is similar to the pristine h-BN, whereas the spectrum of V B differs significantly. These effects essentially imply that some defects are more suited to the modification of the HHG spectra of the bulk materials. This in turn opens the door for tuning HHG by engineering defects in solids.
From our work, we can also estimate of other known defects in h-BN. If one considers 13 the doping impurity instead of vacancy, e.g., carbon impurity, the band-structure of h-BN remains spin-polarised in nature near the band gap [54,65]. If a boron atom is replaced by a carbon atom, one occupied spin-up and one unoccupied spin-down defect levels appear.
Both the defect levels are near the conduction band with the wavefunctions contributed from the p z orbitals of carbon and nearby nitrogen atoms. So, the carbon doping defect is expected to show qualitatively similar behaviour in the HHG spectra as that observed in the case of V N . On the other hand, one occupied spin-up and one unoccupied spin-down states appear near the conduction bands when nitrogen is replaced with a carbon atom.
The wavefunction in this case is contributed mostly by the p z orbitals of the carbon as well as nearby boron atoms. In this case, the enhancement in the below band-gap harmonics is expected due to the defect states near the valence bands similar to V B . However, the effect of screening is expected to be lower compared to V B as the nature of wavefunctions here is similar to that of V N .
In the case of bi-vacancy in h-BN, one occupied defect state exists near the valence band and two unoccupied defect states near the conduction band (see Ref. [65]). In this situation, if the separation between defect states and the nearby bands is small compared to the photon energy, no significant changes in the below band-gap harmonics are expected.
Let us finally comment on the appealing possibility of performing imaging of spin-polarised defects in solids using HHG. As spin channels are not equivalent in the studied defects, one might think about using circularly polarised pulses to probe each spin channel independently. However, in the context of dilute magnetic impurities, as studied here, a crystal will host as many defects with positive magnetic moments as the negative ones, and the signals for up or down spin channel will appear as identical after macroscopic averaging.
Our work opens up interesting perspectives for further studies on strong-field electron dynamics in two-dimensional and extended systems, especially involving isolated defects. Further work may address the possibility of monitoring electron-impurity scattering using HHG, more complex defects such as bi-vacancies, and a practical scheme for imaging buried defects in solids.

Geometry relaxation
Geometry optimisation was performed using the DFT code Quantum ESPRESSO [73,74].
Both the atomic coordinate and the lattice constant relaxation were allowed. Forces were optimised to be below 10 −3 eV/Å. We used an energy cutoff of 150 Ry., and a k-point grid of 10×10×1 k-points. We used a vacuum region of more than 20Å to isolate the monolayer from its periodic copies. The relaxed lattice constants of 12.60 and 12.57Å for V B , and 12.48 A for V N were obtained while for pristine h-BN, it was found to be 12.56Å. The structure of h-BN with a boron vacancy was relaxed by lowering the local threefold symmetry [54]. The lowering of the symmetry with the boron vacancy is attributed to the Jahn-Teller distortion, which is found to be independent of the defect concentration [54,65,69].

TDDFT simulations
By propagating the Kohn-Sham equations within TDDFT, the evolution of the timedependent current is computed by using the Octopus package [75,76]. TDDFT within generalised gradient approximation (GGA) with exchange and correlations of Perdew-Burke-Ernzerhof (PBE) [77] is used for all the simulations presented here. The adiabatic approximation is used for all the time-dependent simulations. We used norm-conserving pseudopotentials. The real-space cell was sampled with a grid spacing of 0.18Å, and a 6×6 k-points grid was used to sample the 2D Brillouin zone. The semi-periodic boundary conditions are employed. A simulation box of 74.08Å along the nonperiodic dimension, which includes 21.17Å of absorbing regions on each side of the monolayer, is used. The absorbing boundaries are treated using the complex absorbing potential method and the cap height h is taken as h = -1 atomic units (a.u.) to avoid the reflection error in the spectral region of interest [78].
Effective single-particle band-structure or spectral function for the vacancy structures is visualised by unfolding the band-structure of 5×5 supercell with 50 atoms [79][80][81]. The spin-polarised calculations are used to address the spin-polarised vacancies.
The Fourier-transform of the total spin-polarised time-dependent electronic current j σ (r, t) is used to simulate the HHG spectrum as where FT and σ stand for the Fourier-transform and the spin-index, respectively.
The total number of excited electrons is obtained by projecting the time-evolved wavefunctions (|ψ n (t) ) on the basis of the ground-state wavefunctions (|ψ GS n ) where N e is the total number of electrons in the system and N k is the total number of kpoints used to sample the BZ. The sum over the band indices n and n run over all occupied states.

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

DATA AVAILABILITY
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.

COMPETING FINANCIAL INTERESTS
The Authors declare no Competing Financial or Non-Financial Interests.