# Controlling high-harmonic generation from strain engineered monolayer
phosphorene

Saibabu Madas ID ELI ALPS, ELI-HU Non-Profit Ltd., Wolfgang Sandner utca 3.,
Szeged 6728, Hungary  Institute of Physics, University of Szeged, Dóm tér 9,
H-6720 Szeged, Hungary  Subhendu Kahaly ID ELI ALPS, ELI-HU Non-Profit Ltd.,
Wolfgang Sandner utca 3., Szeged 6728, Hungary  Mousumi Upadhyay Kahaly ID ELI
ALPS, ELI-HU Non-Profit Ltd., Wolfgang Sandner utca 3., Szeged 6728, Hungary

## Abstract

Phosphorene, a well-studied 2D allotrope of phosphorus, features unique
properties such as widely tunable bandgap, high carrier mobility, and
remarkable intrinsic in-plane anisotropy. Utilizing these structural and
electronic properties, we investigate ultrafast electron dynamics and high
harmonic generation (HHG) from phosphorene subject to band structure
engineering through external strain, based on ab initio time-dependent
density-functional theory approach. We show that strong field processes in
such systems can be optimized and controlled by biaxial tensile and
compressive strain engineering, that results in electronic structure
modification. While $-10\%$ strain resulted in closing of band gap, $2\%$
strain increased the gap by $22\%$ with respect to $0.9\leavevmode\nobreak\
eV$ in pristine phosphorene, consequently affecting the high harmonic yield.
With reduction of gap, by applying strain from $2\%$ to $-10\%$, the valence
band near $\Gamma-$point becomes more flat and discreet, resulting in large
electronic density of states and enhanced electronic excitation, which
reflects in their ultrafast sub-cycle dynamics under laser excitation.
Moreover, due to its intrinsic in-plane anisotropy, harmonic yield with laser
polarization along the armchair (AC) direction is found to be higher than that
of the zigzag (ZZ) direction for all the strain cases. Nearly, an order of
magnitude enhancement of harmonic intensity is achieved for $-10\%$ strain
along AC direction. The current study expands the research possibilities of
phosphorene into a previously unexplored domain, indicating its potential for
future utilization in extreme-ultraviolet and attosecond nanophotonics, and
also for efficient table-top HHG sources.

##  1 Introduction

The advent of high-harmonic generation (HHG) from atomic system [1] had
ignited massive interest in both experimental [2, 3] and theoretical [4]
attosecond science research community. It created the possibility of extending
the then existing techniques utilized in the attosecond science research field
for atomic systems to condensed phase systems. HHG from solids is an emerging
field in attosecond science and material physics [5, 6, 7, 8], demonstrating
great prospects as compact attosecond light-source technologies [9, 10], and
also as novel ultrafast spectroscopy techniques capable to probe electronic
structures [11]. Different solid structure exhibit their own distinct band
structure [12], thereby presenting diverse prospects for interband and
intraband dynamics [13], which have been identified as crucial factors at the
fundamental level for solid-state HHG [14]. However, the comprehensive
understanding of the coupling mechanism between interband and intraband
processes remains incomplete and difficult to be understood with conventional
atomic three-step model often used to address HHG processes from atomic
systems.

The mechanism of HHG in condensed phase systems has main contribution from
interband and intraband processes: (i) an electron-hole pair is generated due
to the excitation of an electron from the valence band to the conduction band,
usually by means of tunnelling near the band gap minimum, resulting in
interband transition. (ii) The solid crystal momentum will follow the time-
dependent vector potential (according to the acceleration theorem [15, 16])
resulting in carrier motion in the bands, which will lead to nonperturbative
intraband harmonics emission [5]. (iii) Simultaneously, interband harmonic
radiation can occur due to electron recombination, with frequencies
corresponding to the instantaneous band gap [17]. The complicated relationship
between intraband and interband processes can be controlled and utilized to
enhance the HHG emission [18, 19]. However, they are non-trivially coupled,
with the interplay between them being either co-operative or competitive [17,
20]; hence their contribution to the overall HHG yield should be investigated.
In this aspect, quantum confinement in solid-state systems can be used as an
additional useful means to control interband and intraband dynamics [21], in
parallel to optical tuning methods, such as tuning laser field polarization
direction [22] and ellipticity [23].

Quantum confined systems [24, 25], for instance, two dimensional (2D)
material’s optical and electronic properties are very sensitive to the
external perturbations, because of their atomic thickness [26]. Amongst wide
range of 2D materials synthesized in recent past [27], monolayer black
phosphorene (BP) displays band gap between that of graphene
($0\leavevmode\nobreak\ eV$) [28] and transition-metal dichalcogenides (TMDs)
($1\leavevmode\nobreak\ eV-2\leavevmode\nobreak\ eV$) [29], resulting in
efficient photo-response in the infrared (IR) range. Hence, due to its highly
anisotropic characteristics [30], together with its high carrier mobility,
significant band gap, and excellent mechanical properties that can withstand
significant amount of deformation before failure, monolayer black phosphorene
is emerging as a viable contender in the competitive field of 2D materials, in
particular for optoelectronics. Regardless of the fact that phosphorene is
highly unstable in ambient temperature, with several emerging ambient
stabilization techniques for phosphorene [31], this system offers a new test
bed for exploring its potential for IR optoelectronics and high harmonic
spectroscopy. One interesting feature of this 2D material is the possibility
to tune the band gap by structural engineering [24, 12]. The effects of
different types of strains (uni-axial or biaxial) along the two anisotropic
lattice vectors on the ground state properties of phosphorene using density
functional theory (DFT) [32] have been extensively studied so far [33, 34,
35]. Experimentally, biaxial tensile strain can be applied through
differential thermal expansion method that relies on a large difference in
thermal expansion coefficients between the 2D material and the substrate. When
a 2D material strongly bonds to the substrate and the system is heated or
cooled, the expansion difference of these materials result in the creation of
a uniform tensile biaxial strain [36]. This thermal strain will be evenly
distributed in two perpendicular directions if the substrate is either
amorphous, cubic, or possesses a (0001) surface in a hexagonal crystal, or a
(001) surface in a tetragonal crystal. In addition, piezoelectric materials
under the influence of an external electric field can be used for applying
biaxial compressive and tensile strains in 2D materials [37, 38, 39]. In this
method, a 2D material is placed on a hybrid substrate that includes a PMN-PT
(lead magnesium niobate-lead titanate) layer, which experiences thickness
variations when subjected to an external electric field. With the application
of an appropriate electrical bias, the substrate extends vertically and
contracts horizontally. This, in turn, induces a uniform biaxial compressive
strain onto the 2D materials. On the other hand, when the bias direction is
reversed, it can induce tensile strain as well. Phosphorene can withstand
tensile strains up to 30% and 27% along the armchair (AC) and zigzag (ZZ)
direction, respectively [40]. Hence, strain engineering can be a practicable
way to tune the electronic properties of phosphorene. Can the utilization of
strain be a potent means to manipulate ultrafast electron dynamics and
consequent optical properties in such materials?

In this study, using a full ab initio approach based on DFT [32] and real-time
time-dependent density functional theory (TDDFT) [41], we demonstrate the
strong-field nonlinear response and the pronounced sensitivity of the HHG
process to the modulation of band structures in monolayer black phosphorene,
that can be achieved through strain engineering. We find that substantial
enhancement of harmonic intensity by up to an order of magnitude is possible
for some compressive strained cases and explore the possible origins behind
the response/phenomena. We investigate the effects of intrinsic anisotropy of
phosphorene by applying electric field along the AC and ZZ directions on the
HHG and total harmonic yield. Our results indicating the tunability of HHG and
total yield through both strain engineering and in-plane anisotropy can be
useful for exploring phosphorene towards extreme-ultraviolet and attosecond
nanophotonics.

This article is organized in the following way. The simulation approaches and
parametric details are described in section 2. Then, we present our findings
for the pristine and strained phosphorene systems, their laser-induced
response and HHG properties in section 3. In section 4 conclusions are
summarized.

Figure 1: Atomic structure of monolayer black phosphorene with it’s
crystallographic directions x (Armchair), y (zigzag). Purple and red colored
spheres represent phosphorous atoms located in different planes within a
single puckered layer. (a) Top view with representative unit cell structure
(black dashed rectangle) and biaxial compressive (- ve) and tensile (+ ve)
strains denoted by red and dark green arrows, respectively. (b) Side view with
geometrical parameters: $d_{1}$ and $d_{2}$ are the bond lengths, $d_{b}$ is
the buckling length. $\alpha$ and $\beta$ are the bond angles. (c) Electronic
band structure along high symmetry directions (green vertical dashed lines),
and (d) corresponding DOS of relaxed pristine (black solid curve), $-10\%$
compressive (red dashed curve), $2\%$ tensile (blue curve) strained
phosphorene. Fermi level ($E_{F}$) is shifted to zero in (c) and (d). Inset in
(c) shows irreducible Brillouin zone (BZ) with high symmetry points of
phosphorene.

##  2 Method

In our work, structural optimization and electronic structure calculations of
pristine and strained phosphorene systems were carried out within DFT [32]
using a plane-wave basis set and ultrasoft pseudo-potential, as implemented in
quantum espresso (QE) [42] package. The generalized gradient approximation
(GGA) is used for the exchange and correlation energy as proposed by Perdew,
Burke, and Ernzerhof (PBE) [43]. To simulate a monolayer of black phosphorous,
a rectangular unit cell with periodic boundary conditions is created. To avoid
the interaction between system and the periodic images of slabs along the
z-direction, fairly large vacuum of 16 Å is created. The wave functions were
expanded in plane waves up to an energy cutoff of $60\leavevmode\nobreak\ eV$.
For Brillouin zone (BZ) integration, we used the k-point sets generated by the
$12\times 12\times 1$ Monkhorst-Pack [44] mesh. All atoms are relaxed till the
final force exerted on each atom is less than $0.001\leavevmode\nobreak\
eV/\text{\AA}$ and also until the electronic energy convergence of
$10^{-6}\leavevmode\nobreak\ eV$ is reached.

At room temperature, the most stable phosphorus allotrope is the black
phosphorus, which was synthesized by Bridgman in 1914 [45]. Black phosphorus
is a layered material in which each atomic layers are stacked together by the
vdW interactions. Unlike graphene, which is a planar layer of carbon atoms,
black phosphorene (monolayer of black phosphorus) is a non-planar layer of
phosphorous atoms. Four phosphorous atoms are present in the unit cell, among
them, two atoms lie in one plane and the rest lie in another parallel plane,
which results in a puckered honeycomb structure with each phosphorous atom
covalently bonding with three adjacent atoms, as a result it has an
anisotropic crystal structure, as shown in the Figure 1(a,b). In our
calculations, we have optimized the lattice constants of black phosphorene,
$a=4.62\leavevmode\nobreak\ \text{\AA}$ and $b=3.29\leavevmode\nobreak\
\text{\AA}$, which are well in agreement with the experiment [46] and the
existing literature [47, 48]. The relaxed geometrical parameters of
phosphorene (shown in Figure 1b), such as the in-plane ($d_{1}$), out-of-plane
($d_{2}$) bond lengths, inter-planar distance or the buckling length
($d_{b}$), angle among in-plane atoms ($\alpha$), and the angle among out-of-
plane atoms ($\beta$) are tabulated in Table 1.

The optimized pristine and strained phosphorene systems from QE calculations
are then used in Octopus package [49, 50] for further time-dependent
calculations, after ensuring force and energy convergence criteria are met. In
order to investigate the laser-driven time dependent electron dynamics in
pristine and the strained phosphorene systems, we employed semi-periodic
boundary conditions, with a simulation box size of 30 Bohr along the non-
periodic direction (along z axis), which includes 3 Bohr of absorbing regions
on both the sides of the phosphorene. This is a precautionary measure to
prevent unphysical reflection of field-accelerated electrons at the border of
the simulation box. The absorbing regions are treated using the complex
absorbing potential (CAP) method [51], with a CAP height of -1 atomic units
(a.u.). The real-space cell was sampled along all the three directions by a
uniform grid spacing of 0.46 Bohr and $38\times 46$ Monkhorst-Pack [44]
k-points grid was used to sample the 2D Brillouin zone. To study the nonlinear
processes, the time evolution of the wave functions and the estimation of
time-dependent electronic current were calculated by propagating the Kohn-Sham
equations in real time and real space within TDDFT implemented within Octopus
package [49, 50], within GGA PBE exchange and correlations [43].

##  3 Results and Discussion

In this work, we applied biaxial strain by tuning the in-plane lattice
constants a and b of pristine phosphorene [Figure 1(a)]. Two types of biaxial
strains are considered in this study: compressive (- ve) and tensile (+ ve)
along both x or armchair (AC) and y or zigzag (ZZ) directions, as shown in
Figure 1(a). Strain ($\varepsilon$) is defined as
$\varepsilon=(a_{\varepsilon}-a_{0})/a_{0}$, where $a_{0}$ and
$a_{\varepsilon}$ are the pristine and strained phosphorene lattice constants,
respectively. As it is evident from Table 1, applied strains affect the
geometrical parameters significantly: $d_{1}$, $d_{2}$, $d_{b}$, $\alpha$ and
$\beta$ are slightly changed due to the application of tensile strain
($\varepsilon=2\%$), where the lattice constants a and b increases due to the
elongation. On the other hand, for $\varepsilon=-10\%$ strain, the lattice
constants a and b decreases resulting in significant change of the geometrical
parameters. Note that such strains are well below the damage threshold as
reported in the existing literature [40]. As tabulated in Table 1, compared to
the pristine case, when $-10\%$ strain is applied, $d_{2}$ and $d_{b}$
increased by 0.8% and 4.7%, respectively, while $d_{1}$, $\alpha$ and $\beta$
decreased by 4%, 6.7% and 7.7%, respectively. It is evident that the
compressive strain has a strong influence on the phosphorene system.
Subsequently, the atomic structure is re-optimized and the electronic
properties of the elongated or compressed phosphorene lattice is compared to
that of pristine case.

Now we explore the differences in electronic band structure and highlight the
important changes in the electronic properties as a function of the applied
biaxial strain. For relaxed black phosphorene structure, we obtain a direct
band gap of $0.9\leavevmode\nobreak\ eV$ at $\Gamma-$point, while the reported
band gap from experiments is between $2.05-2.2\leavevmode\nobreak\ eV$ [52].
As shown in Figure 1(c) (in black) for pristine phosphorene, the maximum
contribution to the valence band maximum (VBM) and conduction band minimum
(CBM) near the $\Gamma-$point is from the $p_{z}$ orbital, as revealed through
projected bands and partial densities of states (PDOS) analysis [Figure SI
1(a,d)]. For $\varepsilon=2\%$, the bandgap increases to
$1.1\leavevmode\nobreak\ eV$, with slight downshift of both valence bands (VB)
and conduction bands (CB). Both the VBM and CBM originating still from $p_{z}$
orbital at $\Gamma-$point [see Figure SI 1(c,f)]. However, for
$\varepsilon=-10\%$, the compressive strain prominently shortens the in-plane
interatomic P-P bond length ($d_{1}$) and increases the in-plane $p_{x,y}$
orbital contributions. Consequently, superposition of the atomic orbitals
leads to shift in the energy of the states and an indirect band gap closure
due to the crossing of VBM and CBM at $\Gamma\longrightarrow X$ path as
denoted with red dashed lines in Figure 1(c). As projected bands and PDOS
analysis revealed that while CBM still mainly originates from the $p_{z}$
orbitals [see Figure SI 1(b,e)] at $\Gamma-$point, the in-plane orbitals
$p_{x}\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ p_{y}$ contributions
become dominant in the VBM along the $\Gamma\longrightarrow X$ path and
crosses the Fermi level [SI 1(b,e)]. Consequently, more electronic density of
states (DOS) near the VBM and CBM are obtained as compared to the other two
cases, as shown in Figure 1(d). These effects result in compressive strain-
induced semiconductor to metal transition in phosphorene. The calculated band
gap values of pristine and strained phosphorene systems (tabulated in Table 1)
are well in agreement with the existing literature [53].

max width=1  Strain _a_ (Å) _b_ (Å) _d 1_ (Å) _d 2_ (Å) _d b_ (Å) $\alpha$
(deg) $\beta$ (deg) $E_{g}$ (eV) 0% 4.62 3.29 2.22 2.26 2.11 104.08 95.6 0.9
-10% 4.16 2.96 2.13 2.28 2.21 97.07 88.15 - 2% 4.72 3.36 2.24 2.26 2.09 104.76
97 1.1

Table 1: Relaxed geometrical parameters and band gap ($E_{g}$) values of
pristine and strained phosphorene systems. Parameters such as lattice
constants _a_ , _b_ , bond lengths $d_{1}$, $d_{2}$, buckling length $d_{b}$,
bond angles $\alpha$, $\beta$ are depicted in Figure 1(b).

Such strain-induced local structural modifications [54] and resulting
electronic properties has direct impact on the time-dependent nonlinear
responses in phosphorene. To probe this, we consider a driving laser field
which is spatially uniform and is described in the velocity gauge, i.e. in
terms of a vector potential [$\boldsymbol{A}(t)$] to treat purely in-plane
perturbation, given as

$\boldsymbol{A}(t)=-\frac{cE_{0}}{\omega_{0}}f(t)\cos{(\omega_{0}t+\phi)}\leavevmode\nobreak\
,$ (1)

where c is the speed of light, $E{{}_{0}}$ is the electric field amplitude and
$\omega_{0}$ is the carrier frequency of the driving laser pulse. The envelop
of the vector potential $f(t)$ considered in this work is sine squared, that
is defined as $f(t)=sin^{2}(\pi t/2\tau)$ for $0<t<2\tau$, or else $f(t)=0$.
$\tau=20\leavevmode\nobreak\ fs$ is the pulse duration at full width half
maximum (FWHM). $E_{0}$ in atomic units is defined in terms of peak intensity
of the laser ($I_{0}$) and is given as $E_{0}=\sqrt{I_{0}/3.51\times
10^{16}}$, where $I_{0}=2\times 10^{11}\leavevmode\nobreak\ Wcm^{-2}$. All
calculations are performed with a central wavelength ($\lambda$) of
$1600\leavevmode\nobreak\ nm$ (corresponding to a photon energy of
$0.77\leavevmode\nobreak\ eV$), and with a zero carrier envelope phase
($\phi=0$). The corresponding electric field is defined as
$\boldsymbol{E}(t)=-\partial\boldsymbol{A}/(c\leavevmode\nobreak\ \partial
t)$. The applied laser waveform, i.e. the vector potential and the
corresponding electric field are depicted in Figure 2(a) and Figure 2(b),
respectively. The pulse parameters are chosen to represent Mid-IR laser regime
as available in ELI ALPS facility [55, 56] and is suitable to probe HHG from
solid-state materials. Note that we assume dipole approximation and the
contributions from magnetic component of the electromagnetic field is ignored.
Also any other relativistic terms, such as spin-orbit coupling are not
considered within our adopted framework. Our systems are excited by the laser
polarized either along AC or ZZ direction, where $\overline{\Gamma X}$ and
$\overline{\Gamma Y}$ shown in the BZ (Figure 1c), corresponds to the AC and
ZZ directions in real-space, respectively. Henceforth, $\overline{\Gamma X}$
and $\overline{\Gamma Y}$ are used as a substitute for AC and ZZ directions,
respectively.

Figure 2: (a) The applied vector potential $\boldsymbol{A}(t)$ and the
corresponding (b) electric field $\boldsymbol{E}(t)$ along with the induced
electronic currents in pristine phosphorene that are parallel ($\parallel$)
and perpendicular ($\perp$) to the incident laser polarization direction for
the AC (c) and ZZ (d) configurations. The arrangement of (e) armchair (AC) and
(f) zigzag (ZZ) configuration defined with respect to the laser polarization
or fundamental field direction and the harmonic fields. The pink and blue
arrows refers to the fundamental field and the harmonic fields, respectively.
The grey dashed line represents the mirror plane (M′) of the phosphorene. Note
that the current follows the vector potential; at peak electric field, the
current is zero for both AC and ZZ in (c,d).

In order to calculate the HHG spectra, we obtain total electronic current in
the system, $\boldsymbol{J}(t)=\partial/\partial
t\int_{\Omega}d^{3}\boldsymbol{r}\boldsymbol{j}(\boldsymbol{r},t)$, which is
obtained by integrating the microscopic electronic current density
[$\boldsymbol{j}(\boldsymbol{r},t)$] over the unit cell, where $\Omega$ is the
unit cell volume [57]. The microscopic electronic current
$\boldsymbol{j}(\boldsymbol{r},t)$ can be written as

$\frac{\partial}{\partial
t}\int_{\Omega}d^{3}\boldsymbol{r}\boldsymbol{j}(\boldsymbol{r},t)=-\int_{\Omega}d^{3}\boldsymbol{r}n(\boldsymbol{r},t)\nabla\upsilon(\boldsymbol{r},t)\leavevmode\nobreak\
,$ (2)

where $n(\boldsymbol{r},t)$ is the time-dependent electron density of the
material driven by laser field. External potential
$\upsilon(\boldsymbol{r},t)$ corresponds to both the electron-nuclei potential
and the applied laser field. The HHG spectrum [$I(\omega)$] is obtained by
applying a discrete Fourier-transform ($\mathcal{FT}$) to $\boldsymbol{J}(t)$,
given as

$I(\omega)=\bigg{|}\mathcal{FT}\Big{(}\boldsymbol{J}(t)\Big{)}\bigg{|}^{2}\leavevmode\nobreak\
.$ (3)

Time dependent dipole currents can capture the first signature of the
anisotropic nonlinear response. With respect to the crystal symmetry, we
analyze the polarization components of the harmonic radiation that are
perpendicular and parallel to the linearly polarized fundamental field. For
any particular laser polarization direction (along AC or ZZ direction),
prominent electronic currents (using Equation 2) and corresponding harmonic
radiation are obtained for both directions parallel and perpendicular to the
laser polarization direction, due to the anisotropy in phosphorene crystal
(elaborated in the next paragraph). The resulting electronic currents for
pristine phosphorene are shown in Figure 2(c) and Figure 2(d). In Figure 2(c),
the electronic current parallel to the electric field direction (along x) is
much higher than the perpendicular contribution (along y). Similarly, in
Figure 2(d), when the laser polarization is along ZZ direction, the parallel
current component (along y) is much greater than the perpendicular current
component (along x). From Figure 2(c,d), it is evident that the parallel
electronic current component is nearly three orders of magnitude higher than
the perpendicular component. The laser polarization direction with respect to
the crystal orientation is shown schematically in Figure 2(e,f). As shown, in
AC configuration, the laser is polarized along AC direction, and likewise for
ZZ. Each of these configurations consists of two different scenarios in which,
w.r.t. the laser polarization direction, the parallel ($\parallel$) and
perpendicular ($\perp$) polarization components of the harmonic radiation is
obtained. When the laser is polarized along ZZ direction, for
$\varepsilon=-10\%$, due to gross change in the crystal structure (as shown in
Table 1) and consequently the electronic properties, the two current
components are noted to have similar/comparable amplitude (see supplementary
Figure SI 2 ).

HHG signal in both parallel and perpendicular directions w.r.t. laser
polarization in phosphorene stems from its intrinsic anisotropy, which
manifests inversion and mirror reflection symmetries only along the ZZ
direction, as demonstrated in Figure 2(e,f). Due to this inversion symmetry,
only odd harmonics are obtained parallel to the laser polarization direction
[58, 59]. This symmetry conditions and symmetry-forbidden selection rules
result in electronic currents along both x and y, when incident laser
polarization is along AC; however the perpendicular current component is 3
orders of magnitude less than the parallel current component along x (Figure
2c).

###  3.1 Mechanism of HHG

HHG spectra for pristine (Figure 3(a,d)), -10% (Figure 3(b,e)), and 2% (Figure
3(c,f)) strained phosphorene systems driven by laser polarized along AC
(Figure 3(a-c)) and ZZ (Figure 3(d-f)) direction is investigated. We show the
harmonic resolved in parallel and perpendicular configurations corresponding
to the harmonics reported in Figure 4(a,b). Anisotropic HHG response is
clearly visible by the distinctive harmonic radiation obtained for AC and ZZ
cases in all the panels. In all the cases except, -10% ZZ, the parallel
component of the harmonics is much higher than the perpendicular component as
shown in Figure SI 2.

Figure 3: HHG spectra of pristine phosphorene (a,d), $-10\%$ strained (b,e),
and $2\%$ strained (c,f) phosphorene systems driven by a laser polarized along
AC (a,b,c) and along ZZ (d,e,f). Figure 4: HHG spectra of pristine (dashed
line), $-10\%$ (dash dotted line) and $2\%$ (dotted line) strained phosphorene
systems driven by the laser polarized along AC (a) and ZZ (b) direction.

Different features in HHG spectrum resolved in parallel and perpendicular
directions for AC and ZZ configurations [in Figure 4(a,b)] originates from the
intrinsic anisotropy of the phosphorene crystal. For the parallel harmonic
components, i.e., x in AC case [black lines in Figure 3(a-c)] and y in ZZ case
(blue lines in Figure 3(d-f)], for both pristine and strained structures, the
fundamental peak appears at $\omega$ ($=0.77\leavevmode\nobreak\ eV$),
followed by odd integer multiples of the driving laser frequency [see Figure
4(a,b)], suggesting that the inversion symmetry of the system is protected
even by applying biaxial compressive and tensile strains. In each case, the
HHG peak intensity is found to monotonically decrease, and is distinctly
visible until H13. In both AC and ZZ cases, compressive strain
($\varepsilon=-10\%$) leads to significant enhancement in the harmonic
intensity compared to the pristine case. Additionally, in AC case, for
$\varepsilon=-10\%$, HHG peaks become stronger and spectral width broadens, as
can be seen in Figure 4(a), while for ZZ there is no systematic enhancement.
In all the cases, except $\varepsilon=-10\%$ in ZZ case, the parallel
component of the harmonics is higher than the perpendicular component.
Interestingly, in -10% case, from the finite electronic population at the
Fermi level [as evident along $\overline{\Gamma X}$ and $\overline{\Gamma Y}$
directions in Figure 1(c,d)], some states are always available for field-
induced excitations. Thus, even when the field is along ZZ, prominent current
and HHG signal can be obtained along perpendicular direction. As found in
Figure SI 1(b), along $\overline{\Gamma X}$ (i.e. AC direction), electronic
states at and around the Fermi level ($E_{F}$) contribute significantly to
field induced transition, thereby resulting in higher strain-induced HHG
enhancement. For example, in the AC configuration for $-10\%$ case, the high
harmonic spectrum [from H5 in Figure 4(a)] clearly shows the harmonic
intensity enhancement by nearly an order of magnitude. In particular for -10%
case, due to significant structural changes (see Table 1), the effective
electron-nuclei potential is prominently modified (in Equation 2), resulting
in distinctive new features in the electronic band dispersion, thereby
affecting the microscopic electronic current, and eventually the HHG signal.
The lower the states are in the conduction band, the easier it is for the free
carriers to get excited by the laser field from the VBM to these lowered
states, which eventually leads to an enhanced HHG [60]. Furthermore, due to
the applied strain, the conduction bands and the valence bands shift
downwards. This downward shift of bands is more prominent in the case of
$\varepsilon=-10\%$ compared to $\varepsilon=2\%$. However, in
$\varepsilon=-10\%$ case, the available electronic states in the downward
shifted VBM along $\overline{\Gamma Y}$ are low compared to that in
$\varepsilon=2\%$ case. As a consequence, the harmonic yield for the ZZ
configuration is expected to be higher in the case of $\varepsilon=2\%$
compared with the rest of the cases, which is investigated and discussed
later. In order to clearly understand these results, we show the computed
number of electrons excited to the conduction band from the valence band
during the laser pulse in Figure 5 corresponding to Figure 4.

It is well known that in a semiconductor, recombination of an electron with an
existing hole in the valence band results in the emission of a photon with
energy exceeding the band gap energy. Consequently, it becomes evident that
harmonics below the band gap do not stem from the interband contribution.
Hence, below bandgap harmonics are only from intraband current. However, above
band gap harmonics can originate from both interband and intraband currents.
In case of phosphorene, the band gap is $0.9\leavevmode\nobreak\ eV$, and
induced photon energy is $0.77\leavevmode\nobreak\ eV$. Thus, in this
scenario, recombination of electrons and holes in different bands becomes
possible, and both interband and intraband processes are allowed to jointly
contribute to the generation of harmonics above the band gap. To estimate the
extent of interband transitions in our systems of study, we calculated the
temporal evolution of total number of excited electrons from valence band to
conduction band by projecting the time-evolved wave functions
($\ket{\psi_{i,\textbf{k}}(t)}$) on to the basis of ground state wave
functions ($\ket{\psi_{j,\textbf{k}}^{GS}}$), given as [61]

$N_{ex}(t)=N_{e}-\frac{1}{N_{\textbf{k}}}\sum_{i,j}^{occ}\sum_{\textbf{k}}^{BZ}|\langle\psi_{i,\textbf{k}}(t)|\psi_{j,\textbf{k}}^{GS}\rangle|^{2}\leavevmode\nobreak\
,$ (4)

where $N_{e}$, $N_{\textbf{k}}$, $\psi_{i,\textbf{k}}$, and
$\psi_{j,\textbf{k}}$ are the total number of electrons in the system, the
total number of k-points used to sample the BZ, Kohn-Sham state at the i-th
band at the k-point k, and the ground Kohn-Sham state at the j-th band at the
k-point k, respectively. The sum over the band indices $i$ and $j$ run over
all occupied states and the summation runs over the entire first Brillouin
zone.

Figure 5: Total number of excited electrons from valence bands to the
conduction bands during the laser pulse, which is polarized along AC (a) and
along ZZ (b) in pristine (in blue), $-10\%$ (in red), and $2\%$ (in black)
strained phosphorene.

Thus, $N_{ex}(t)$ indicates the amount of interband current. Temporal
evolution of $N_{ex}$ displays same sub-cycle periodicity as the induced laser
field (Figure 5(a) for AC , Figure 5(b) for ZZ). Figure 5 shows that the
number of excited electrons is larger in the case of AC-polarized laser, which
can hence lead to larger electronic current and higher harmonic yield along AC
(consistent with Figure 4). It is important to note that the anisotropic
electronic band dispersion shows flatter bands along $\overline{\Gamma Y}$
(ZZ), in comparison to $\overline{\Gamma X}$ (AC direction). This implies
lower effective mass and higher carrier mobility in AC direction, in
comparison to ZZ. Consequently interband current is prominently larger in AC
direction, as demonstrated by corresponding $N_{ex}$ (compare both panels in
Figure 5). Consequently, significant difference in the charge excitation along
different directions is expected. For -10% strain in AC case, the laser
induced excited electrons reaches a peak value of 0.12 $e^{-}$/atom (0.48
$e^{-}$/unit cell) at $\sim 21\leavevmode\nobreak\ fs$ and then oscillates at
around 0.035 $e^{-}$/atom (0.14 $e^{-}$/unit cell) [Figure 5(a)]. Similar
value of $e^{-}$/atom is obtained for 0% and 2% strain cases, with 0% slightly
higher than 2%. For ZZ configuration, the maximum excited $e^{-}$/atom of 0.10
is attained for -10% strain case and then oscillates at 0.02 by the end of the
laser pulse. However, for 0% and 2% strain cases, the excited electrons
(slightly higher in 2% case compared to 0% case) returns to almost zero after
each half cycle in the ZZ configuration, while for AC configuration for all
cases, a fraction of excited electrons remain effectively excited to
conduction bands at the end of the laser pulse ($\sim 40\leavevmode\nobreak\
fs$). These strain-induced number of excited electrons are well in agreement
with the corresponding harmonic spectra shown in Figure 4. Only for -10%
strain case, for both laser polarization along AC and ZZ direction, during
each half cycle, more electrons are excited, however fraction of excited
electrons do not return to ground state. This indicates that carriers have
been driven to higher energy levels (interband) and then remained there even
after the removal of the laser pulse. These excited electrons must occupy
broader areas of BZ, i.e., larger momentum space, thereby resulting in higher
intraband current under compressive strain, which contribute to the harmonic
generation collectively with interband current.

Figure 6: Time-dependent electron dynamics in strain-free and strained
phosphorene systems. (a) Vector potential of the driving laser pulse polarized
along AC (_x_ -direction). Snap shots of the sub-cycle dynamics of momentum-
resolved excited electrons near the peak and minimum of the vector potential
(marked as A, B, C in (a)) for pristine phosphorene (b,c,d), -10% (e,f,g), and
2% (h,i,j) strained phosphorene. The oscillation features of the excited
electrons in the BZ in panels (b-d) for pristine, (e-g) for -10%, and (h-j)
for +2% show evidence of intraband process significantly contributes to HHG.

Furthermore, to validate our interpretation and to obtain laser-driven
subcycle electron dynamics, we calculate the momentum-resolved excited
electrons distribution $N_{ex}(\textbf{k};t)$ in the following way [61]:

$N_{ex}(\textbf{k};t)=\frac{1}{N_{\textbf{k}}}\bigg{(}N_{e}-\sum_{i,j}^{occ}|\langle\psi_{i,\textbf{k}}(t)|\psi_{j,\textbf{k}}^{GS}\rangle|^{2}\bigg{)}\leavevmode\nobreak\
.$ (5)

We compare the momentum-resolved sub-cycle dynamics of excited electrons at
three time instances, which are near the maximum and minimum of the vector
potential, marked as A, B, C in Figure 6(a) for 0% (top), $-10\%$ (middle),
and $2\%$ (bottom) strained phosphorene systems [Figure 6(b-j)]. It is
interesting to observe that the excited electron density is concentrated at
the $\Gamma-$point for all the cases. For both 0% and 2% strained phosphorene,
both VBM and CBM are at $\Gamma-$point, hence electron transition probability
is maximum at $\Gamma$ (due to their similar band structure). On the other
hand, for -10% strained system, VBM and CBM overlap near $\Gamma-$point (see
$\overline{\Gamma X}$ band dispersion in Figure 1c). This results in the
higher concentration of excited electron density at and around the
$\Gamma-$point in this particular case. At time instant denoted by B in Figure
6(a), i.e. when the system has the lowest excitation, electrons wander to
broader areas in BZ under compressive strain [see Figure 6(f)]. Hence, the
number of excited electrons at each momentum ($\textbf{k}_{x},\textbf{k}_{y}$)
increases, which confirms that intraband and interband processes cooperatively
enhances the harmonic yield for -10% case. At the time instants indicated by A
and C in Figure 6(a), a large part of BZ has larger amount of excited
electrons, and it is clear that the excited electrons arrange themselves from
left(right) to right(left) boundaries in the BZ. These sub-cycle oscillations
of excited electrons driven by the laser field along with excited electrons
traversing to larger parts of BZ shows that intraband and interband processes
collectively contributes to HHG and enhanced harmonic yield compared to
pristine and strained phosphorene cases, which is highlighted later.

To reveal the nature of the HHG spectrum, we calculated the peak harmonic
intensity ($I_{n}^{max}$) of representative harmonics ($n$) as a function of
peak laser intensity ($I_{L}$). Figure 7 illustrates the scaling behavior of
harmonic peaks $n=$ 5, 7, 9 and 11 as a function of $I_{L}$, ranging from
$0.1-0.4\leavevmode\nobreak\ TWcm^{-2}$. In all the configurations shown in
Figure 7, for $n=5,\leavevmode\nobreak\ 7$, fit to the power law $I^{p}$
resulted in _p_ with a value that ranges between $2-5$. Whereas, for the high-
order harmonics, for instance H9 and H11, _p_ ranges between $5-7$. By fitting
the data to the power-law, the non-perturbative nature of the high-order
harmonics is clearly evident in the scaling of individual harmonic peaks with
$I_{L}$. Note that in the perturbative regime, the radiation is limited to
frequencies that are merely a few times higher than the frequency of the
driving field. Furthermore, the $n^{th}$ harmonic would scale proportionally
with $I^{n}_{L}$, where $I_{L}$ is the peak intensity of the laser and the
harmonic yield would reduce exponentially [14, 58].

Figure 7: Calculated peak harmonic intensity as a function of peak laser
intensity $I_{L}$ for 0% strain (a,d), -10% (b,e), and 2% strain (c,f) for
representative harmonics n $=5$ (blue circles), 7 (orange circles), 9 (green
circles), and 11 (red circles) obtained for laser polarized along AC (a-c) and
along ZZ (d-f). Solid lines are obtained by fitting data to power law,
yielding exponents showing the non-perturbative scaling of the harmonic
process. The resulting exponents are mentioned in the legends of individual
panel.  Figure 8: The relative change in harmonic intensity as a function of
applied strain for representative harmonics for (a) AC case and (b) ZZ case.
(c) Harmonic yield (integrated from H2 to H15) as a function of applied strain
for both AC and ZZ cases.

For each harmonic peak, the relative change in the harmonic intensity and
total harmonic yield are compared for AC and ZZ cases, as shown in Figure 8.
The relative change in the harmonic intensity is defined as
$I^{\varepsilon}_{n}/I^{\varepsilon=0}_{n}$, where $I^{\varepsilon}_{n}$ and
$I^{\varepsilon=0}_{n}$ stands for $n^{th}$ order harmonic peak intensity with
strain and without strain, respectively. Consistent with HHG signals (in
Figure 4), for laser polarized along AC (Figure 8a), the compressive strain of
-10% enhances HHG intensities by $100-400\%$ for all harmonic orders, H7 being
most prominently affected followed by H9, H13, H11 [see Figure 8(a)]. With
tensile strain, however, the relative harmonic intensity of all the harmonics,
except for $5^{th}$ order, decreased, as seen in Figure 8(a). With tensile
strain, however, the relative harmonic intensity of all the harmonics, except
for $5^{th}$ order, decreased, as seen in Figure 8(a). For laser polarization
along ZZ, each harmonic responded differently with respect to the applied
strain (Figure 8(b)), reflecting the lack of systematic variation of HHG
intensity. Different HHG peak intensities and harmonic yield for AC and ZZ
cases is consistent with anisotropic band dispersion of phosphorene and
distinctive features. Harmonic yield (integrated from 2nd order to 15th order)
as a function of applied strain (Figure 8(c)) shows almost 3 times enhancement
in total yield when $-10\%$ strain is applied for AC case. This demonstrates
the possibility of tailoring the HHG yield using both appropriate strain and
laser electric field direction.

##  4 Conclusions

In summary, using state-of-the-art TDDFT simulations, we have demonstrated the
effects of ultrashort laser induced non-linear response and yield of non-
perturbative HHG from a monolayer black phosphorene, extending to the 13th
order. In particular the role of biaxial compressive and tensile strain on the
high harmonic generation from monolayer phosphorene is investigated. Our
predictions show that tuning the crystal structure of phosphorene via biaxial
strain provides an efficient tool to tailor its electronic structure and HHG
in terms of harmonic intensity and spectral width. Due to the intrinsic in-
plane anisotropy of phosphorene, the harmonic yield for the AC configuration
is higher compared to that obtained for the ZZ configuration. The laser pulse
parameters as chosen in our study are in Mid-IR regime, that are suitable to
probe HHG from solid-state band-gap materials, and available in advanced laser
facilities like ELI ALPS facility [55, 56]. Our analysis shows that, the
harmonic yield can be enhanced by a factor of 3 by applying $-10\%$
compressive strain. The origin of this unique response can be attributed to
the electronic band structure engineering through external compressive strain,
that results in prominent increase in the number of excited electrons, which
eventually leads to an enhanced HHG. Our results and insights promotes the
possibility to tune the HHG process by band engineering through the external
application of strain and also opens the door to experimentally investigate
strain-induced strong-field nonlinear response of this promising 2D
nanostructure.

In light of the growing enthusiasm for developing compact table-top sources of
coherent extreme ultraviolet XUV radiation, which hold promise for practical
applications across various domains of science and technology, and given the
ever-expanding spectrum of potential applications for 2D materials, our study
showcases the promise of phosphorene and analogous 2D nanostructures. With
easily tunable band structures feasible through strain engineering, these 2D
materials exhibit significant potential for designing and optimizing
innovative devices in the field of attoscience, as well as for the creation of
stable table-top HHG sources.

## Declaration of Competing Interest

The authors confirm that they have no competing financial interests or
personal relationships that could have influenced the findings presented in
this paper.

## Acknowledgements

The ELI ALPS project (GINOP-2.3.6-15-2015-00001) is supported by the European
Union and it is co-financed by the European Regional Development Fund. This
research has been supported by the IMPULSE project which receives funding from
the European Union Framework Programme for Research and Innovation Horizon
2020 under grant agreement No 871161. SK, MUK and SM also acknowledges project
No. 2019-2.1.13-TÉT-IN-2020-00059, which has been implemented with support
provided by the National Research, Development and Innovation Fund of Hungary,
and financed under the 2019-2.1.13-TÉT-IN funding scheme. MUK and SM would
like to thank the ELI ALPS IT department and the HPC administration for their
support in providing computational resources. This research is supported by
the ÚNKP-23-4 - New National Excellence Program of the Ministry for Culture
and Innovation from the Source of the National Research, Development and
Innovation Fund.

## References

  * [1] A. McPherson et al.  “Studies of multiphoton production of vacuum-ultraviolet radiation in the rare gases”  In _Journal of the Optical Society of America B_ 4.4 The Optical Society, 1987, pp. 595  DOI: 10.1364/josab.4.000595
  * [2] P.. Corkum  “Plasma perspective on strong field multiphoton ionization”  In _Physical Review Letters_ 71.13 American Physical Society (APS), 1993, pp. 1994–1997  DOI: 10.1103/physrevlett.71.1994
  * [3] Jeffrey L. Krause, Kenneth J. Schafer and Kenneth C. Kulander  “High-order harmonic generation from atoms and ions in the high intensity regime”  In _Physical Review Letters_ 68.24 American Physical Society (APS), 1992, pp. 3535–3538  DOI: 10.1103/physrevlett.68.3535
  * [4] M. Lewenstein et al.  “Theory of high-harmonic generation by low-frequency laser fields”  In _Physical Review A_ 49.3 American Physical Society (APS), 1994, pp. 2117–2132  DOI: 10.1103/physreva.49.2117
  * [5] Shambhu Ghimire et al.  “Observation of high-order harmonic generation in a bulk crystal”  In _Nature Physics_ 7.2 Springer ScienceBusiness Media LLC, 2010, pp. 138–141  DOI: 10.1038/nphys1847
  * [6] T.. Luu et al.  “Extreme ultraviolet high-harmonic spectroscopy of solids”  In _Nature_ 521.7553 Springer ScienceBusiness Media LLC, 2015, pp. 498–502  DOI: 10.1038/nature14456
  * [7] Yong Sing You, David A. Reis and Shambhu Ghimire  “Anisotropic high-harmonic generation in bulk crystals”  In _Nature Physics_ 13.4 Springer ScienceBusiness Media LLC, 2016, pp. 345–349  DOI: 10.1038/nphys3955
  * [8] Seunghwoi Han et al.  “High-harmonic generation by field enhanced femtosecond pulses in metal-sapphire nanostructure”  In _Nature Communications_ 7.1 Springer ScienceBusiness Media LLC, 2016  DOI: 10.1038/ncomms13105
  * [9] Katsumi Midorikawa  “Progress on table-top isolated attosecond light sources”  In _Nature Photonics_ 16.4 Springer ScienceBusiness Media LLC, 2022, pp. 267–278  DOI: 10.1038/s41566-022-00961-9
  * [10] Arjun Nayak et al.  “Saddle point approaches in strong field physics and generation of attosecond pulses”  In _Physics Reports_ 833 Elsevier BV, 2019, pp. 1–52  DOI: 10.1016/j.physrep.2019.10.002
  * [11] Andrea Damascelli  “Probing the Electronic Structure of Complex Systems by ARPES”  In _Physica Scripta_ T109 IOP Publishing, 2004, pp. 61  DOI: 10.1238/physica.topical.109a00061
  * [12] Mousumi Upadhyay Kahaly and Umesh V. Waghmare  “Contrast in the Electronic and Magnetic Properties of Doped Carbon and Boron Nitride Nanotubes: A First-Principles Study”  In _The Journal of Physical Chemistry C_ 112.10 American Chemical Society (ACS), 2008, pp. 3464–3472  DOI: 10.1021/jp072340d
  * [13] T.. Kaloni, M. Kahaly, Y.. Cheng and U. Schwingenschlögl  “Ge-intercalated graphene: The origin of the p-type to n-type transition”  In _EPL (Europhysics Letters)_ 99.5 IOP Publishing, 2012, pp. 57002  DOI: 10.1209/0295-5075/99/57002
  * [14] Eleftherios Goulielmakis and Thomas Brabec  “High harmonic generation in condensed matter”  In _Nature Photonics_ 16.6 Springer ScienceBusiness Media LLC, 2022, pp. 411–421  DOI: 10.1038/s41566-022-00988-y
  * [15] Felix Bloch  “Über die Quantenmechanik der Elektronen in Kristallgittern”  In _Zeitschrift für Physik_ 52.7-8 Springer ScienceBusiness Media LLC, 1929, pp. 555–600  DOI: 10.1007/bf01339455
  * [16] W.. Houston  “Acceleration of Electrons in a Crystal Lattice”  In _Physical Review_ 57.3 American Physical Society (APS), 1940, pp. 184–186  DOI: 10.1103/physrev.57.184
  * [17] G. Vampa et al.  “Theoretical Analysis of High-Harmonic Generation in Solids”  In _Physical Review Letters_ 113.7 American Physical Society (APS), 2014  DOI: 10.1103/physrevlett.113.073901
  * [18] O. Schubert et al.  “Sub-cycle control of terahertz high-harmonic generation by dynamical Bloch oscillations”  In _Nature Photonics_ 8.2 Springer ScienceBusiness Media LLC, 2014, pp. 119–123  DOI: 10.1038/nphoton.2013.349
  * [19] Takuya Higuchi, Mark I. Stockman and Peter Hommelhoff  “Strong-Field Perspective on High-Harmonic Radiation from Bulk Solids”  In _Physical Review Letters_ 113.21 American Physical Society (APS), 2014  DOI: 10.1103/physrevlett.113.213901
  * [20] D. Golde, T. Meier and S.. Koch  “High harmonics generated in semiconductor nanostructures by the coupled dynamics of optical inter- and intraband excitations”  In _Physical Review B_ 77.7 American Physical Society (APS), 2008  DOI: 10.1103/physrevb.77.075330
  * [21] Kotaro Nakagawa et al.  “Size-controlled quantum dots reveal the impact of intraband transitions on high-order harmonic generation in solids”  arXiv, 2022  DOI: 10.48550/ARXIV.2212.11483
  * [22] Keisuke Kaneshima et al.  “Polarization-Resolved Study of High Harmonics from Bulk Semiconductors”  In _Physical Review Letters_ 120.24 American Physical Society (APS), 2018  DOI: 10.1103/physrevlett.120.243903
  * [23] T. Tamaya, A. Ishikawa, T. Ogawa and K. Tanaka  “Higher-order harmonic generation caused by elliptically polarized electric fields in solid-state materials”  In _Physical Review B_ 94.24 American Physical Society (APS), 2016  DOI: 10.1103/physrevb.94.241107
  * [24] M. Kahaly, Shikha Misra and S.. Mishra  “Photo-assisted electron emission from illuminated monolayer graphene”  In _Journal of Applied Physics_ 121.20 AIP Publishing, 2017  DOI: 10.1063/1.4984318
  * [25] Mousumi Upadhyay Kahaly and Umesh V. Waghmare  “Effect of curvature on structures and vibrations of zigzag carbon nanotubes: A first-principles study”  In _Bulletin of Materials Science_ 31.3 Springer ScienceBusiness Media LLC, 2008, pp. 335–341  DOI: 10.1007/s12034-008-0053-x
  * [26] Nicolas Tancogne-Dejean and Angel Rubio  “Atomic-like high-harmonic generation from two-dimensional materials”  In _Science Advances_ 4.2 American Association for the Advancement of Science (AAAS), 2018  DOI: 10.1126/sciadv.aao5207
  * [27] Beatriz Mendoza-Sánchez and Yury Gogotsi  “Synthesis of Two-Dimensional Materials for Capacitive Energy Storage”  In _Advanced Materials_ 28.29 Wiley, 2016, pp. 6104–6135  DOI: 10.1002/adma.201506133
  * [28] K.. Novoselov et al.  “Electric Field Effect in Atomically Thin Carbon Films”  In _Science_ 306.5696 American Association for the Advancement of Science (AAAS), 2004, pp. 666–669  DOI: 10.1126/science.1102896
  * [29] Qing Hua Wang et al.  “Electronics and optoelectronics of two-dimensional transition metal dichalcogenides”  In _Nature Nanotechnology_ 7.11 Springer ScienceBusiness Media LLC, 2012, pp. 699–712  DOI: 10.1038/nnano.2012.193
  * [30] S. Madas, S.. Mishra, S. Kahaly and M. Kahaly  “Superior Photo-thermionic electron Emission from Illuminated Phosphorene Surface”  In _Scientific Reports_ 9.1 Springer ScienceBusiness Media LLC, 2019  DOI: 10.1038/s41598-019-44823-x
  * [31] Yonghong Zeng and Zhinan Guo  “Synthesis and stabilization of black phosphorus and phosphorene: Recent progress and perspectives”  In _iScience_ 24.10 Elsevier BV, 2021, pp. 103116  DOI: 10.1016/j.isci.2021.103116
  * [32] W. Kohn and L.. Sham  “Self-Consistent Equations Including Exchange and Correlation Effects”  In _Physical Review_ 140.4A American Physical Society (APS), 1965, pp. A1133–A1138  DOI: 10.1103/physrev.140.a1133
  * [33] Baisheng Sa et al.  “Strain Engineering for Phosphorene: The Potential Application as a Photocatalyst”  In _The Journal of Physical Chemistry C_ 118.46 American Chemical Society (ACS), 2014, pp. 26560–26568  DOI: 10.1021/jp508618t
  * [34] Junbeom Seo, Sungwoo Jung and Mincheol Shin  “The Performance of Uniaxially Strained Phosphorene Tunneling Field- Effect Transistors”  In _IEEE Electron Device Letters_ 38.8 Institute of ElectricalElectronics Engineers (IEEE), 2017, pp. 1150–1152  DOI: 10.1109/led.2017.2712259
  * [35] Mohsen Yarmohammadi, Mohammad Mortezaei and Kavoos Mirabbaszadeh  “Anisotropic basic electronic properties of strained black phosphorene”  In _Physica E: Low-dimensional Systems and Nanostructures_ 124 Elsevier BV, 2020, pp. 114323  DOI: 10.1016/j.physe.2020.114323
  * [36] Geun Ho Ahn et al.  “Strain-engineered growth of two-dimensional materials”  In _Nature Communications_ 8.1 Springer ScienceBusiness Media LLC, 2017  DOI: 10.1038/s41467-017-00516-5
  * [37] Yeung Yu Hui et al.  “Exceptional Tunability of Band Energy in a Compressively Strained Trilayer MoS2 Sheet”  In _ACS Nano_ 7.8 American Chemical Society (ACS), 2013, pp. 7126–7131  DOI: 10.1021/nn4024834
  * [38] Fei Ding et al.  “Stretchable Graphene: A Close Look at Fundamental Parameters through Biaxial Straining”  In _Nano Letters_ 10.9 American Chemical Society (ACS), 2010, pp. 3453–3458  DOI: 10.1021/nl101533x
  * [39] Wenjing Jie et al.  “Effects of controllable biaxial strain on the Raman spectra of monolayer graphene prepared by chemical vapor deposition”  In _Applied Physics Letters_ 102.22 AIP Publishing, 2013, pp. 223112  DOI: 10.1063/1.4809922
  * [40] Qun Wei and Xihong Peng  “Superior mechanical flexibility of phosphorene and few-layer black phosphorus”  In _Applied Physics Letters_ 104.25 AIP Publishing, 2014, pp. 251915  DOI: 10.1063/1.4885215
  * [41] M.A.L. Marques and E.K.U. Gross  “TIME-DEPENDENT DENSITY FUNCTIONAL THEORY”  In _Annual Review of Physical Chemistry_ 55.1 Annual Reviews, 2004, pp. 427–455  DOI: 10.1146/annurev.physchem.55.091602.094449
  * [42] Paolo Giannozzi et al.  “QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials”  In _Journal of Physics: Condensed Matter_ 21.39 IOP Publishing, 2009, pp. 395502  DOI: 10.1088/0953-8984/21/39/395502
  * [43] John P. Perdew, Kieron Burke and Matthias Ernzerhof  “Generalized Gradient Approximation Made Simple”  In _Physical Review Letters_ 77.18 American Physical Society (APS), 1996, pp. 3865–3868  DOI: 10.1103/physrevlett.77.3865
  * [44] Hendrik J. Monkhorst and James D. Pack  “Special points for Brillouin-zone integrations”  In _Physical Review B_ 13.12 American Physical Society (APS), 1976, pp. 5188–5192  DOI: 10.1103/physrevb.13.5188
  * [45] P.. Bridgman  “TWO NEW MODIFICATIONS OF PHOSPHORUS.”  In _Journal of the American Chemical Society_ 36.7 American Chemical Society (ACS), 1914, pp. 1344–1363  DOI: 10.1021/ja02184a002
  * [46] L. Cartz et al.  “Effect of pressure on bonding in black phosphorus”  In _The Journal of Chemical Physics_ 71.4 AIP Publishing, 1979, pp. 1718–1721  DOI: 10.1063/1.438523
  * [47] Jingsi Qiao et al.  “High-mobility transport anisotropy and linear dichroism in few-layer black phosphorus”  In _Nature Communications_ 5.1 Springer ScienceBusiness Media LLC, 2014  DOI: 10.1038/ncomms5475
  * [48] Poya Yasaei et al.  “High-Quality Black Phosphorus Atomic Layers by Liquid-Phase Exfoliation”  In _Advanced Materials_ 27.11 Wiley, 2015, pp. 1887–1892  DOI: 10.1002/adma.201405150
  * [49] Xavier Andrade et al.  “Real-space grids and the Octopus code as tools for the development of new simulation approaches for electronic systems”  In _Physical Chemistry Chemical Physics_ 17.47 Royal Society of Chemistry (RSC), 2015, pp. 31371–31396  DOI: 10.1039/c5cp00351b
  * [50] Xavier Andrade et al.  “Time-dependent density-functional theory in massively parallel computer architectures: the octopus project”  In _Journal of Physics: Condensed Matter_ 24.23 IOP Publishing, 2012, pp. 233202  DOI: 10.1088/0953-8984/24/23/233202
  * [51] Umberto De Giovannini, Ask Hjorth Larsen and Angel Rubio  “Modeling electron dynamics coupled to continuum states in finite volumes with absorbing boundaries”  In _The European Physical Journal B_ 88.3 Springer ScienceBusiness Media LLC, 2015  DOI: 10.1140/epjb/e2015-50808-0
  * [52] Liangbo Liang et al.  “Electronic Bandgap and Edge Reconstruction in Phosphorene Materials”  In _Nano Letters_ 14.11 American Chemical Society (ACS), 2014, pp. 6400–6406  DOI: 10.1021/nl502892t
  * [53] Can Wang, Qinglin Xia, Yaozhuang Nie and Guanghua Guo  “Strain-induced gap transition and anisotropic Dirac-like cones in monolayer and bilayer phosphorene”  In _Journal of Applied Physics_ 117.12 AIP Publishing, 2015  DOI: 10.1063/1.4916254
  * [54] M. Kahaly, S. Nazir and U. Schwingenschlögl  “Band structure engineering and vacancy induced metallicity at the GaAs-AlAs interface”  In _Applied Physics Letters_ 99.12 AIP Publishing, 2011  DOI: 10.1063/1.3643049
  * [55] Sergei Kühn et al.  “The ELI-ALPS facility: the next generation of attosecond sources”  In _Journal of Physics B: Atomic, Molecular and Optical Physics_ 50.13 IOP Publishing, 2017, pp. 132002  DOI: 10.1088/1361-6455/aa6ee8
  * [56] Dimitris Charalambidis et al.  “The Extreme Light Infrastructure—Attosecond Light Pulse Source (ELI-ALPS) Project”  In _Springer Series in Chemical Physics_ Springer International Publishing, 2017, pp. 181–218  DOI: 10.1007/978-3-319-64840-8˙10
  * [57] Nicolas Tancogne-Dejean, Oliver D. Mücke, Franz X. Kärtner and Angel Rubio  “Impact of the Electronic Band Structure in High-Harmonic Generation Spectra of Solids”  In _Physical Review Letters_ 118.8 American Physical Society (APS), 2017  DOI: 10.1103/physrevlett.118.087403
  * [58] Hanzhe Liu et al.  “High-harmonic generation from an atomically thin semiconductor”  In _Nature Physics_ 13.3 Springer ScienceBusiness Media LLC, 2016, pp. 262–265  DOI: 10.1038/nphys3946
  * [59] Lei Jia et al.  “Optical high-order harmonic generation as a structural characterization tool”  In _Physical Review B_ 101.14 American Physical Society (APS), 2020  DOI: 10.1103/physrevb.101.144304
  * [60] Rui Qin and Zi-Yu Chen  “Strain-controlled high harmonic generation with Dirac fermions in silicene”  In _Nanoscale_ 10.47 Royal Society of Chemistry (RSC), 2018, pp. 22593–22600  DOI: 10.1039/c8nr07572g
  * [61] Nicolas Tancogne-Dejean, Oliver D. Mücke, Franz X. Kärtner and Angel Rubio  “Ellipticity dependence of high-harmonic generation in solids originating from coupled intraband and interband dynamics”  In _Nature Communications_ 8.1 Springer ScienceBusiness Media LLC, 2017  DOI: 10.1038/s41467-017-00764-5

##  Supplementary Information Supplementary Information 1 Projected bands and
partial density of states

Figure SI 1: (a-c) Projected band structure and (d-f) partial density of
states of 0%, -10%, and 2% strained phosphorene systems.

##  Supplementary Information Supplementary Information 2 Induced electronic
currents in -10% strained phosphorene.

Figure SI 2: Induced electronic currents in $-10\%$ strained phosphorene that
are perpendicular (x) and parallel (y) to the incident laser polarization
along ZZ direction.