Also at ]Quantum Optics, Quantum Nanophysics and Quantum Information Faculty
of Physics, University of Vienna, Boltzmanngasse 5 1090, Vienna, Austria

# Generation of Half-Integer Harmonics and Efficient THz-to-Visible Frequency
Conversion in Strained Graphene

Marco Ornigotti<EMAIL_ADDRESS>Faculty of Engineering and Natural
Sciences, Photonics, Tampere University, Tampere, FI-33720 Finland  Luca
Ornigotti  [  Department of Optics, Palacký University, 17. listopadu 1192/12,
711 46 Olomouc, Czech Republic  Fabio Biancalana  School of Engineering and
Physical Sciences, Heriot-Watt University, Edinburgh, UK

###### Abstract

We study the generation of harmonics from graphene under the influence of an
artificial magnetic field, generated via bending of a graphene flake. We show
how the Landau level structure induced by the pseudomagnetic field breaks the
centrosymmetry of graphene, thus allowing the generation of even harmonics. We
also show, that depending on the impinging pulse duration, the nonlinear
signal does not only contain the integer harmonics of the impinging pulse, but
also its half-integer ones, due to the peculiar square-root-like nature of
Landau levels in graphene.

††preprint: AIP/123-QED

##  I introduction

Gauge fields are ubiquitous in Nature, and regulate the dynamics of several
different fields of physics. Perhaps the most common example of a gauge field
is the electromagnetic field, whose quantum, the photon, mediates the
interaction between charged matter ref1 . Besides electrodynamics, gauge
fields make their appearance in the Standard Model of particle physics, for
example, as a unified way to describe the interaction of matter with the
fundamental forces of Nature, excluding gravityref2 . In quantum mechanics,
the gauge invariance of the electromagnetic field is at the core of the
celebrated Aharonov-Bohm effect ref3 , which paved the way for a deeper
understanding of gauge fields in terms of fibre bundles ref4 ; ref5 ; ref6 .
Moreover, non-Abelian gauge theories are also an essential ingredient to
understand geometrical and Berry phases ref7 ; ref8 . In condensed matter
physics, gauge fields play a crucial role in understanding long-range
interaction, and the emergence of collective phenomena, such as the appearance
of Abrikosov vortices ref9  or topological states of matter ref5 ; ref10 ;
ref11 .

In recent years, artificial (or synthetic) gauge fields (AGF) have started to
attract considerable attention. Contrary to gauge fields, which arise from
_real_ fields or geometric connections, AGFs can be tuned ad-hoc by acting on
a physical system in a certain, predetermined and tunable way. Typically, the
physical mechanism generating AGFs is also different, in nature, than the one
generating the actual gauge field it mimics. For example, an artificial
magnetic field is experiences by cold atoms in a rotating frame, due to the
mathematical equivalence between the Coriolis force and the Lorentz force
ref13 , or for light propagating in a waveguide with a twisted propagation
direction, due to the topological equivalence between the twisting and an
actual magnetic field applied to the waveguide ref14 . The advent of AGFs
paved the way for completely new research fields, which merge together aspects
from topology and set theory with different concepts and methods from various
physical disciplines, such as topological mechanics ref15 , topological
condensed matter physics ref16 , topological atomic physics ref17 ; ref18 ;
ref19 ; ref20 , and topological photonics ref21 ; ref22 , where both concepts
from lattice gauge field theories ref23  and condensed matter physics ref16
were used to create novel ways to control the propagation and coupling of
light in wave-guiding structures, such as photonic topological insulators
ref24 ; ref24bis ; ref25 , non-Abelian-like dynamics in engineered waveguide
lattices ref26 , AGF switching using the angular momentum of light ref27 , and
topological protection ref28 ; ref29 ; ref30 ; ref31 .

Aside from photonics, another platform that in the last decade represented a
rich playground for testing the effects of different classes of AGFs has been
graphene. Since its experimental discovery in 2004 by Novoselov and Geim ref32
, graphene has in fact attracted considerable attention not only for its
peculiar band structure ref33 , its anomalous quantum Hall effect ref34  or
its minimal conductivity ref35  and universal absorbance ref35a , but it also
represented an unexpected connection between condensed matter physics and
gauge field theory. The application of strain, stress, or bending on graphene,
in fact, gives rise to AGFs in the form of effective electric and magnetic
fields ref36 ; ref37 ; ref38 ; ref39 . Moreover, out-of-plane bending of
graphene flakes is analog, for the electrons in graphene, to consider their
evolution on a curved background, under the action of gravity ref40 ; ref41 .

Graphene also represents a very interesting platform for photonics, mainly for
its very large nonlinear response, compared to bulk materials. Recent
theoretical extra1 ; extra7  and experimental extra8 ; extra9 ; extra10
results, in fact, estimate the bulk-equivalent third-order nonlinear
susceptibility of graphene to go from $\chi_{3D}^{(3)}\simeq 10^{-15}$
$m^{2}/V^{2}$ in the visible region extra8  up to $\chi_{3D}^{(3)}\simeq
10^{-9}$ $m^{2}/V^{2}$ in the terahertz region extra9 . For comparison, the
bulk third-order nonlinearity of a typical silica glass is of the order of
$\chi_{crystal}^{(3)}\simeq 10^{-23}$ $m^{2}/V^{2}$ (see Ref. extra11 ).
Graphene, therefore, possesses an extremely high nonlinear response, several
orders of magnitude higher than that of a normal nonlinear material . A
typical way to obtain such results experimentally is by means of the Z-scan
technique, which is a standard characterisation method for the nonlinear
response of bulk materials, and can be applied for graphene as well under
certain assumptions (see Ref. extra10  for details). The results of such
measurement leads naturally to a 3D susceptibility, but in case of graphene,
these numbers should be taken with a grain of salt, since the notion of bulk
susceptibility for monolayer graphene does not make much sense (graphene is ,
indeed, a true 2D material), and they should be used only for comparison with
other materials extra17 .

The nonlinear response of graphene in the limit of strong magnetic field
(i.e., when the magnetic length is much smaller than the wavelength of the
impinging electromagnetic pulse) has also been estimated theoretically extra12
to be inversely proportional to the external magnetic field, i.e.,
$\chi_{3D}^{(3)}\simeq 5\times 10^{-9}/B(T)$ $m^{2}/V^{2}$. In the same work
extra12 , moreover, the maximum intensity of the nonlinear signal has also
been estimated to be linearly proportional to the applied magnetic field,
i.e., $I_{3}^{(max)}\simeq B$. From these results we see how controlling the
magnetic field results in an overall increase in the intensity of the
nonlinear signal. This feature, combined with the possibility of creating
pseudomagnetic fields through AGFs, could be exploited to reduce the necessary
pump intensity to trigger nonlinear phenomena in graphene, and might lead, in
the future, to a new generation of integrated nonlinear devices.

Amidst the vast literature on AGFs in graphene, some works have been focusing
the attention on their effects on the interaction of graphene with external
electromagnetic fields, pointing out how a constant external electric field
can drastically modify the arrangement of Landau levels and edge states of
bent graphene flakes, introducing a squeezing of the Landau states in graphene
ref42 , while a time-dependent electric field can induce modulations of the
angular momentum transfer between light and graphene ref43 , or induce
coherent population transfer between different Landau levels of the valence
and conduction band of graphene ref44 . To the best of our knowledge, however,
the influence of AGFs on the nonlinear response of graphene interacting with
an ultrashort laser pulse has not been investigated yet.

In this work, therefore, we present a comprehensive analysis of the
interaction of ultrashort light pulses with graphene, in the presence of an
AGF. In particular, we consider the case of strained graphene, as presented in
Ref. ref39 , which implements an effective constant, uniform, magnetic field
orthogonal to the graphene plane. The presence of an AGF induces Landau levels
in both the valence and conduction band of graphene, thus introducing
selection rules on the impinging pulse polarisation. Under these assumption,
we study different interaction configurations, for impinging pulses of
different resonant frequencies, and we show, that the nonlinear signal
produced by electrons near the Dirac point of graphene differ substantially
from the case of unbent graphene, since the presence of the AGF breaks the
symmetry of graphene, thus allowing the appearance of even harmonics in the
nonlinear optical response.

This work is organised as follows: in Sect. II, we briefly review how bending
graphene introduces an artificial gauge field, and what are the consequences
of that for an electron in graphene. Then, we present the theoretical
framework necessary to describe the interaction of ultrashort pulses with
graphene in the presence of AGF in Sect. III. Section IV is dedicated to the
main results of our work, i.e., the calculation and discussion of the
nonlinear response of bent graphene to an external ultrashort pulse. Finally,
conclusions and future prospects are discussed in Sect. VI.

##  II Artificial Gauge Field from Bending

When a mechanical strain is applied to graphene, the immediate result is the
appearance of an artificial gauge field, whose spatial distribution and
orientation depends on the nature of the strain applied ref46 . The general
expression for the induced AGF in the absence of out-of-plane modulations can
be then written as ref39

$\textbf{A}^{(s)}=\pm\frac{s\beta}{a}\left(u_{xx}-u_{yy}\right)\hat{\textbf{x}}\mp\frac{2s\beta}{a}u_{xy}\hat{\textbf{y}},$
(1)

Figure 1: (a) Pictorial representation of a rectangular graphene flake
deformed into an arc. The radii of the lower and upper edge are, respectively,
$R$, and $R^{\prime}$. With a graphene flake of width $W=200$ nm, and length
$L=192$ nm, and an inner and outer radius of, respectively, $R=5L=960$ nm, and
$R^{\prime}=R+W=1.16$ $\mu$m, the maximum achievable magnitude of the
pseudomagnetic field in the central region of the flake is $B=10$ T ref39 .
(b) Flattened equivalent geometry for the bent graphene flake in panel (a).
The curvature induced by the bending is replaced with an artificial gauge
field, which gives rise to a uniform pseudomagnetic field
$\textbf{B}=B\,\hat{\textbf{z}}$ parallel to the $z$-axis. The width $w$ and
length $\ell$ of the flattened flake might be extended to infinity, without
changing the essential role of the pseudomagnetic field. in the interaction
dynamics of the flake with an external electromagnetic pulse.

where $a\simeq 1.42$ $\mathrm{\SIUnitSymbolAngstrom}$ is the carbon-carbon
interatomic distance, $\beta=\simeq 2$ is the electron Grüneisen parameter
ref47 , $s$ accounts for the strength of the strain, and
$u_{\mu\nu}=\left(\partial_{\mu}u_{\nu}+\partial_{\nu}u_{\mu}\right)$ is the
strain tensor, with $u_{\mu}$ being the displacement vector. Notice, that
strain also induces a scalar potential $V(x,y)\propto u_{xx}+u_{yy}$, which
can be anyway neglected, by choosing an appropriate gauge, where the scalar
potential is set to zero. A careful choice of the strain tensor can lead to
different bending and deforming geometries, corresponding to different AGFs.
Amongst the various choices available ref46 , we choose the bending profile
discussed in Ref. ref39 , which allows the creation of a uniform
pseudomagnetic field in a rectangular graphene flake by introducing only one
deformation parameter, i.e., the bending radius $R$. The explicit expression
of the displacement vector that implements this geometry is then given as

$\displaystyle u_{x}(x,y)$ $\displaystyle=\frac{xy}{R},$ (2a) $\displaystyle
u_{y}(x,y)$ $\displaystyle=-\frac{x^{2}}{2R},$ (2b)

which corresponds to an AGF $\textbf{A}^{(s)}=\pm s\beta
y/aR\hat{\textbf{x}}\equiv-By\hat{\textbf{x}}$, and, therefore, to a
pseudomagnetic field $\textbf{B}=B\hat{\textbf{z}}$. A schematic
representation of the bent graphene flake, together with a set of
experimentally realisable parameters, is given in Fig. 1. Notice, that for the
choice of parameters as in Fig. 1, a uniform magnetic field of magnitude
$B=10$ T can be generated within the flake. In general, however, since the
magnitude of the pseudomagnetic field is inversely proportional to the bending
radius of the graphene flake, i.e., $B=s\beta/R$, a smaller bending radius,
i.e., a bigger bending angle, will result in a higher pseudomagnetic field. It
is also worth noticing, that introducing a strain also has two other effects,
namely it shifts the position of the Dirac points in $k$ space by a quantity
proportional to $\textbf{A}^{(s)}$ ref47 , and it also renders Fermi velocity
anisotropic ref48 , according to the relation
$\textbf{v}_{F}=v_{F0}\left(\textbf{I}-\beta\bar{\textbf{u}}+\bar{\textbf{u}}\right)$,
where $\bar{\textbf{u}}$ is the strain tensor. However, for the purpose of
this work, we assume that, at the leading order in the bending radius $R$, the
Fermi velocity is not affected by the deformation and it remains approximately
constant.

The dynamics of electrons in graphene in the presence of the AGF
$\textbf{A}^{(s)}$ can be studied by replacing the kinetic momentum p in the
low-energy graphene Hamiltonian ref47
$\mathcal{H}=v_{F}\bm{\sigma}\cdot\textbf{p}$ (with
$\bm{\sigma}=\sigma_{x}\hat{\textbf{x}}+\sigma_{y}\hat{\textbf{y}}$, and
$\sigma_{\mu}$ are Pauli matrices), with the canonical momentum
$\bm{\Pi}=\textbf{p}+e\textbf{A}^{(s)}$ deriving from minimal coupling of the
electron field with the AGF ref2 , which results in the following Dirac
equation for electrons in the vicinity of the Dirac point

$i\hbar\partial_{t}\phi=-i\hbar
v_{F}\sigma^{j}\left(\partial_{j}+\frac{ie}{\hbar}A_{j}^{(s)}\right)\phi$ (3)

where $j=\\{x,y\\}$. The above equation can be solved analytically, as it
represents the well-known problem of a relativistic electron in a magnetic
field, whose solution can be cast in terms of the Landau eigenfunctions

$\phi_{n}^{\pm}(y,t,k)=\mathcal{N}\,e^{i\left(kx-\frac{\mathcal{E}_{n}}{\hbar}t\right)}\left(\begin{array}[]{c}\text{sign}\left(n\right)\phi_{|n|-1}(\xi),\\\
\\\ \phi_{|n|}(\xi)\end{array}\right),$ (4)

Figure 2: (a) Pictorial representation of the band structure of graphene in
the vicinity of a Dirac point (solid cones). The presence of a pseudomagnetic
field generates Landau levels in both the valence band (red lines) and
conduction band (blue lines), as well as at the Dirac point itself (green
point). The energies of the Landau levels created in this manner are the same
in modulus for the two bands, and differ only in sign, i.e., positive energies
are associated with the conduction band, while negative energies with the
valence band. (b) Schematic representation of the level structure induced by
the pseudomagnetic field, and the correspondent selection rules for the first
3 Landau levels in valence and conduction band. The green, solid arrows
correspond to the transition allowed for incoming left-handed circular
polarisation (LHC) photons, while the red, dashed arrow correspond to
transitions allowed for incoming right-handed circular polarisation (RHC)
photons.

where $\mathcal{N}$ is a normalisation constant, which equals one for $n=0$,
and $1/\sqrt{2}$ otherwise, $\phi_{n}(\xi)$ are harmonic oscillator
eigenstates, $\xi=\left(y+L_{c}^{2}k\right)/L_{c}$, with
$L_{c}=\sqrt{\hbar/eB}$ being the magnetic length (for a pseudomagnetic field
of $B=10$ T we get $L_{c}\simeq 8.11$ nm), and
$\mathcal{E}_{n}=\text{sign}\left(n\right)\hbar\omega_{C}\sqrt{|n|}$, with
$\omega_{C}=v_{F}/L_{c}$ being the cyclotron frequency. The above solution is
moreover equipped with the constraint $\phi_{-1}(\xi)=0$ ref47 . A sketch of
the structure of Landau levels in graphene is depicted in Fig. 2. Notice, that
contrary to the case of a non-relativistic electron in magnetic field, where
the spacing of Landau levels is constant, i.e., Landau levels have the full
structure of a harmonic oscillator, in graphene they are not equally spaced,
as the energy eigenvalue $\mathcal{E}_{n}$ scales as $\sqrt{|n|}$. Moreover,
that due to the peculiar band structure of graphene in the vicinity of the
Dirac point, two set of Landau levels are created by the magnetic field,
corresponding to Landau levels in the valence band [associated to negative
values of the index $n$, and appearing as red lines in Fig. 2 (a) ], and in
the conduction band [associated to positive values of the index $n$, and
appearing as blue lines in Fig. 2 (a)]. The two set of oscillator states are
almost disjoint from each other, with the exception of the ground state $n=0$
[green dot in Fig. 2 (a)], which sits at exactly the Dirac point and it is
common to both sets.

The creation of Landau levels in graphene introduces selection rules for the
dipole-allowed transitions, i.e., an incoming photon can only excite the
transition $|n_{i}\rangle\rightarrow|n_{f}\rangle$ if and only if
$|n_{f}|=|n_{i}|\pm 1$, where the plus sign holds for right-handed circularly
polarised photons, while the minus sign holds for left-handed circularly
polarised photons only. These selection rules are a consequence of the natural
spin-orbit coupling arising from the interaction of graphene’s band structure
with the pseudomagnetic field ref49 . The structure of these selection rules
is shown in Fig. 2 (b) for the first three Landau levels in both conduction
and valence band.

##  III Interaction with an Ultrashort Electromagnetic Pulse

To reach our goal, i.e., investigate the nonlinear response of graphene to an
external ultrashort pulse, we first need to understand how to describe its
interaction with an external, time-dependent field. To this aim, let us assume
that a flake of bent graphene is interacting with an external field within the
electric dipole interaction approximation. In this case, therefore, the
equation of motion for an electron in the vicinity of a Dirac point reads as
follows

$i\hbar\partial_{t}\psi=\left[-i\hbar
v_{F}\sigma^{j}\left(\partial_{j}+\frac{ie}{\hbar}A^{(s)}_{j}\right)+e\,E_{j}(t)x^{j}\right].$
(5)

Without any loss of generality, we can assume, that the impinging electric
field is linearly polarised along the $x$-direction, and that is characterised
by a central frequency $\omega_{L}$. We also assume, that the electric field
is normally impinging on the graphene flake, i.e., that the field is
propagating along the $z$ direction, as defined in Fig. 1. We then work in the
Landau gauge, where the scalar potential is zero. In this case, we can
entirely describe the electric field by means of its vector potential
$A(t)=-\int\,dt\,E(t)\equiv\mathcal{A}(t)\exp{[-i\omega_{L}t]}$, where
$\mathcal{A}(t)$ accounts for the temporal shape of the pulse, so that the
total vector potential, given by the combination of the true, $A(t)$, and the
artificial, $A^{(s)}(y)$ gauge field, can be written as
$\textbf{A}(y,t)=\left[A(t)-By\right]\hat{\textbf{x}}$. In this gauge, Dirac
equation reduces to the following, manifestly covariant, compact form

$i\hbar\gamma^{\mu}\left(\partial_{\mu}+\frac{ie}{\hbar}A_{\mu}\right)\psi=0,$
(6)

where $\mu=\\{0,1,2\\}\rightarrow\\{v_{F}t,x,y\\}$, $\gamma^{0}=\sigma^{z}$,
$\gamma^{1}=i\sigma^{y}$, and $\gamma^{2}=-i\sigma^{x}$. In general, this
equation does not admit an explicit analytical solution with an arbitrary
time-dependent vector potential. However, following the procedure described in
Ref. ref35bis , if we know the form of the instantaneous eigenstates for the
above equation, we can then make an educated guess at the true form of the
solution for any time $t$. To this aim, we first take the Fourier transform
with respect to the $x$-variable and operate the following phase
transformation, i.e.,

Figure 3: Temporal evolution of the population of the Landau levels
$|0\rangle$ (black, dot-dahsed line), $|1\rangle$ (blue, dashed line), and
$|-1\rangle$ (red, solid line), for different values of the time duration
$\tau$ of the impinging ultrashort pulse, i.e., (a) $\tau=20$ fs, and (b)
$\tau=50$ fs. The central frequency of the impinging electric field is
$\omega_{L}=\omega_{1}=174$ THz, and the pulse amplitude is $E_{0}=10^{7}$
$V/m$. As it can be seen, as $\tau$ increases, the rate of population transfer
also changes accordingly (longer pulses correspond to a higher Rabi frequency
ref51 ). However, for the selected configuration, the interaction between
levels $|1\rangle$ and $|-1\rangle$ is always mediated by the ground state
$|0\rangle$ at the Dirac point. To obtain these plots, the value $v_{F}=c/300$
m/s for the Fermi velocity has been assumed.

$\psi(x,y,t)=\int\,dk\,e^{ikx}e^{iG(t)\sigma_{x}}\phi(y,t,k),$ (7)

with

$G(t)=\frac{ev_{F}}{\hbar}\int_{0}^{t}\,d\tau\,A(\tau),$ (8)

and notice that Eq. (6) reduces to

$\left[\frac{1}{v_{F}}\gamma^{0}\partial_{t}+\gamma^{1}\left(-ik-\frac{ieB}{\hbar}y\right)+\partial_{y}\right]\phi=0,$
(9)

which is equivalent to Eq. (3). This indicates, that the Landau eigenstates in
Eq. (4) can be taken as the instantaneous eigenstates for the problem at-hand,
and that the general solution of Eq. (6) can be written as

$\psi(y,t,k)=\sum_{n}e^{iG(t)\sigma_{x}}\left[c_{n}^{+}(t)\phi_{n}^{+}(y,k)+c_{n}^{-}(t)\phi_{n}^{-}(y,k)\right].$
(10)

Substituting the above Ansatz into Eq. (6) leads to a set of coupled mode
equations for the expansion coefficients $c_{n}^{\pm}(t)$, which, for a
linearly polarised impinging electromagnetic field read

$\displaystyle\dot{c}_{m}^{+}$
$\displaystyle=i\left|\mathcal{N}\right|^{2}\Omega(t)e^{-i\omega_{L}t}\Big{[}f_{m,m-1}^{-}(t)c_{m-1}^{+}+f_{m,m+1}^{-}(t)c_{m+1}^{+}$
$\displaystyle+\sum_{n<0}f^{+}_{m,n}(t)\langle\phi_{m}^{-}|\sigma^{x}|\phi_{n}^{+}\rangle
c_{n}^{-}\Big{]},$ (11a) $\displaystyle\dot{c}_{m}^{-}$
$\displaystyle=-i\left|\mathcal{N}\right|^{2}\Omega(t)e^{-i\omega_{L}t}\Big{[}f_{m,m-1}^{-,*}(t)c_{m-1}^{+}+f_{m,m+1}^{-,*}(t)c_{m+1}^{+}$
$\displaystyle+\sum_{n<0}f^{+,*}_{m,n}(t)\langle\phi_{m}^{+}|\sigma^{x}|\phi_{n}^{-}\rangle
c_{n}^{-}\Big{]},$ (11b)

where $\Omega(t)=(ev_{F}/\hbar)\mathcal{A}(t)$ is the Rabi frequency,
$f^{\pm}_{m,n}(t)=e^{i(\omega_{m}\pm\omega_{n})t}$ accounts for the eigenvalue
mismatch between the states participating in the temporal evolution of
$c_{m}^{\pm}(t)$, $\omega_{m}=\omega_{C}\sqrt{|m|}$ is the eigenfrequency
associated to the state $|\phi_{m}^{\pm}\rangle$ (notice, that the sign of
$\omega_{m}$ has been explicitly taken care of during the calculations that
lead to the above equations already), and
$\langle\phi_{m}^{-}|\sigma^{x}|\phi_{n}^{+}\rangle$
($\langle\phi_{m}^{+}|\sigma^{x}|\phi_{n}^{-}\rangle$) is the dipole matrix
element of graphene ref49 , whose explicit expression in this case is given by

$\langle\phi_{m}^{-}|\sigma^{x}|\phi_{n}^{+}\rangle=\delta_{|n|,|m|-1}+\delta_{|n|,|m|+1}.$
(12)

These equations are the first result of our work. They, in fact, describe the
interaction of an arbitrarily shaped linearly polarised impinging field.
Notice, how the polarisation of the field enters only in the dipole matrix
element, and therefore the above equation can be easily generalised for an
arbitrary polarisation, by replacing $\sigma^{x}$ in the expression of the
dipole matrix element, with the Pauli matrix (or combination thereof)
correspondent to the impinging polarisation.

We can simplify the coupled mode equations above by assuming that the
impinging electromagnetic field is nearly resonant with one specific
transition. Since the polarisation of the impinging field is linear, rather
than circular, both green and red transitions in Fig. 2(b) will be allowed,
once the frequency of the incoming field has been chosen to match one of the
transitions between the Landau levels. A closer inspection to the structure of
the selection rules depicted in Fig. 2(b), moreover, reveals that the only
nontrivial dynamics that can be induced in bent graphene corresponds to an
impinging field resonant with the transition $|0\rangle\leftrightarrow|\pm
1\rangle$, and a consequent three-level intra-band dynamics, rather than a
simple two-level one, either inter-band or intra-band.

If we then now assume, that the impinging field is resonant with the
$|0\rangle\leftrightarrow|\pm 1\rangle$ transitions, i.e.,
$\omega_{L}=\omega_{1}=\omega_{C}$, Eqs. (11) become

$i\frac{d}{dt}\left(\begin{array}[]{c}c_{1}^{-}\\\ c_{0}\\\
c_{1}^{+}\end{array}\right)=\left(\begin{array}[]{ccc}0&\Gamma(t)&0\\\
\Gamma^{*}(t)&0&-\Gamma(t)\\\
0&-\Gamma^{*}(t)&0\end{array}\right)\left(\begin{array}[]{c}c_{1}^{-}\\\
c_{0}\\\ c_{1}^{+}\end{array}\right).$ (13)

where $\Gamma(t)=\Omega(t)\exp{\left(-i\,\omega_{1}\,t\right)}$. Notice, that
the structure of the coefficient matrix in the equations above admits the
existence of a dark state, corresponding in this case to no population being
at the Dirac point at any given time, i.e., $c_{0}(t)=0$. The consequences of
this, and its benefits for coherent population control of graphene in the
presence of an external magnetic field have been recently investigated in Ref.
ref44 .

Figure 4: Temporal evolution (a) and Fourier transform (b) of the components
of the Dirac current $J_{x}(t)$ (red, solid line), and $J_{y}(t)$ (blue,
dashed line), for an impinging pulse of duration $\tau=20$ fs. As it can be
seen from panel (a), despite the fact the impinging electric field is
polarised along the $x$-direction, a nonzero current is also generated in the
$y$-direction. In particular, panel (b) reveals how while
$\tilde{J}_{y}(\omega)$ contains essentially only a peak at
$\omega=\omega_{1}$, $\tilde{J}_{x}(\omega)$ has a richer structure. The extra
peaks in $\tilde{J}_{y}(\omega)$ are localised around the fundamental
frequency $\omega_{1}$, and are the source of half-integer harmonics seen in
Fig. 5 and 6. For these plots, $\omega_{L}=\omega_{1}=174$ THz, $E_{0}=10^{7}$
V/m, and $v_{F}=c/300$ m/s have been used, to reproduce results compatible
with Fig. 3.

We solve the above equations with the initial condition $c_{1}^{-}(0)=1$,
i.e., with the electron initially in the valence band, and for a vector
potential described by a Gaussian pulse of duration $\tau$, and central
frequency $\omega_{L}=\omega_{1}$, so that $\Omega(t)$ becomes

$\Omega(t)=\frac{ev_{F}E_{0}\tau}{\hbar}e^{-\frac{(t-t_{0})^{2}}{\tau^{2}}}\cos\left(\omega_{L}t\right),$
(14)

where $E_{0}$ is measured in $V/m$ and $t_{0}$ is an arbitrary temporal delay.
The evolution of the expansion coefficients $c_{n}^{\pm}(t)$ for an impinging
pulse of amplitude $E_{0}=10^{7}$ $V/m$, and time duration $\tau=10$ fs,
corresponding to a nearly single-cycle pulse (i.e., $\omega_{L}\tau=1.74$),
and $\tau=50$ fs, corresponding to a pulse with several optical cycles, i.e.,
$\omega_{L}\tau=8.7$, is depicted in Fig. 3 (a) and (b), respectively. Notice,
how a very short pulse ($\tau=10$ fs induces an almost instantaneous change in
the population of the three interested levels, accompained by very faint
Bloch-Siegert oscillations extra15 , as the interaction is too fast, and the
system does not have the time to adapt to it. For longer pulses, as in the
case of Fig. 3 (b), instead, the population dynamics appears more complicated,
but leads, at equilibrium, to a situation in which all population returns back
to the initial state $|-1\rangle$. In both cases, however, the dynamics always
involves the ground state $|0\rangle$ at the Dirac point.

Our simulations, moreover, have been conducted at $T=0$ K, i.e., the
distribution of carriers in the valence and conduction band has not been
explicitly taken into account. However, as it has been discussed by one of the
authors in previous publications extra13 ; extra14 , the impact of temperature
on the nonlinear response of graphene is negligible.

##  IV Nonlinear Signal and Dirac Current

The nonlinear response of graphene can be estimated by evaluating, as a
function of frequency, the intensity of the nonlinear radiation emitted by it,
as a consequence of the interaction with an impinging, time-dependent electric
field, i.e. ref35bis

$I(\omega)\propto\left|\omega\,\tilde{\textbf{J}}(\omega)\right|^{2},$ (15)

where $I(\omega)$ is the spectrum of the emitted radiation, and
$\tilde{\textbf{J}}(\omega)$ is the Fourier transform of the Dirac current
$j^{\mu}(x,y,t)=\bar{\psi}(x,y,t)\gamma^{\mu}\psi(x,y,t)=\psi^{\dagger}(x,y,t)\sigma^{\mu}\psi(x,y,t)$
(where the last equality holds because of our definition of gamma matrices
given in the previous section).

Since we are only interested in the temporal features of the current, we can
integrate it with respect to the transverse space, to obtain

$\textbf{J}(t)=\int\,d^{2}R\,\textbf{j}(x,y,t)=\int\,d^{2}R\,\left(\bm{\sigma}\cdot\textbf{R}\right)\left|\psi(x,y,t)\right|^{2},$
(16)

Figure 5: Nonlinear signal, as defined in Eq. (15), as a function of the
normalised frequency $\omega/\omega_{1}$, for different values of the pulse
length, i.e., (a) $\tau=20$ fs (correspondent to a fluence of $F=2.65$
$m\,J/m^{2}$), (b) $\tau=50$ fs (correspondent to a fluence of $F=6.64$
$m\,J/m^{2}$), and (c) $\tau=100$ fs (correspondent to a fluence of $F=13.27$
$m\,J/m^{2}$). The red, dashed line in each panel represents the position of
the fundamental frequency $\omega_{L}=\omega_{1}$. Notice, that the peak that
should correspond to $\omega=\omega_{1}$ is indeed slightly blue-shifted, due
to the different weight of the $x$\- and $y$-components of the Dirac current,
as it can be seen in Fig. 4(b). This effect only manifests significantly for
short pulses, and tends to disappear for longer ones, as it can be seen by
comparing panel (a) and (c), above. For these plots,
$\omega_{L}=\omega_{1}=174$ THz, $E_{0}=10^{7}$ V/m, and $v_{F}=c/300$ m/s
have been used, to reproduce results compatible with Fig. 3.

Notice, that since the expansion coefficients $c_{n}^{\pm}(t)$ appearing in
the definition of $\psi(x,y,t)$ given by Eq. (10) only depend on time, and
that the solution along the $x$-direction can be expressed in terms of plane
waves, we can safely perform the integration only over the Landau eigenstates,
and exploit their orthogonality relation to compute the spatial integral in
Eq. (16). In doing this, moreover, the resulting time-dependent current
$\textbf{J}(t)$ will have the same functional form for both the case of a
finite graphene flake undergoing bending, and an infinite sheet of graphene
subjected to a uniform magnetic field orthogonal to the graphene plane. The
only difference between the two cases will then be an overall multiplicative
constant, which will account for the actual arrangement of the system. Since
this constant does not change the overall functional for of $\textbf{J}(t)$,
and neither impacts the overall form of its Fourier transform, we then treat
it as a simple normalisation factor, and scale everything to it.

If we then substitute Eq. (10) into the above expression, and limit ourselves
to the case of the impinging pulse being resonant with the
$|0\rangle\leftrightarrow|\pm 1\rangle$, the $x$\- and $y$-components of the
integrated current have the following explicit form

$\displaystyle J_{x}(t)$
$\displaystyle=c_{0}^{*}(t)\left[c_{1}(t)^{+}e^{-i\omega_{1}t}-c_{1}(t)^{-}e^{i\omega_{1}t}\right]+\text{c.c.},$
(17a) $\displaystyle J_{y}(t)$
$\displaystyle=i\,c_{0}^{*}(t)\left[c_{1}(t)^{+}e^{-i\omega_{1}t}-c_{1}(t)^{-}e^{i\omega_{1}t}\right]+\text{c.c.}.$
(17b)

The temporal evolution of the Dirac current, as well as its Fourier transform,
is shown in Fig. 4 for the case of an impinging pulse with $\tau=10$ fs.
Notice how, despite the fact the impinging electric field is polarised along
the $x$-direction, a nonzero component of the current along the $y$-direction
arises. This is a consequence of the broken centrosymmetry, induced by the
artificial gauge field. Since the current enters with its whole vectorial
character in the definition of the nonlinear signal as given by Eq. (15), this
has significant consequences on its spectrum. In fact, for the case of a short
pulse as the one used in Fig. 4, the shape of the nonlinear signal is almost
entirely determined by $\tilde{J}_{x}(\omega)$, while for longer pulse widths,
the interplay between the current components $J_{\mu}(t)$ becomes more
prominent, giving rise to a richer structure of $I(\omega)$.

To prove this, in Fig. 5 we plot the nonlinear signal defined in Eq. (15) for
different values of the impinging pulse width, i.e., $\tau=10$ fs [panel (a)],
$\tau=50$ fs [panel (b)], and $\tau=100$ fs [panel (c)]. As it can be seen,
for very short pulses [panel (a)], a set of equally spaced harmonics appear,
with spacing $n\omega_{1}/2$.This can be explained by noticing, that the
impinging pulse effectively sees a 3-level system with equally spaced levels,
as $|\pm 1\rangle$ have the same distance in energy from $|0\rangle$. This
allows us to approximate the square-root behaviour of the Landau ladder to the
more traditional harmonic oscillator behaviour which, in turn, gives the
equally spaced peaks in Fig. 5(a), with spacing $\omega_{1}/2$.

This is the second result of our work. For short enough pulses, the local
structure of Landau levels around the Dirac point can be approximated with
that of a traditional harmonic oscillator, and thus the corresponding
nonlinear signal contains all the integer and half-integer harmonics of the
pulse carrier frequency $\omega_{L}$.

For longer pulse durations, on the other hand, the nonlinear signal shows a
more rich and complicated spectrum. From Fig. 5(b) and (c), in fact, it is
possible to see how the spectrum broadens, with respect to the situation
depicted in Fig. 5(a), and higher harmonics appear, containing both even and
odd contributions. Notice, moreover, how in both cases of Fig. 5(b) and (c),
although the highest intensity is reached for the 7th [panel (b)] and the 11th
[panel (c)] harmonic, higher ones still have considerable intensity. For
example, for long pulses [panel (c)], the 20th can be generated with
significant intensity.

If we recall that $\omega_{1}=174$ THz, $\omega=20\,\omega_{1}$ would
correspond to radiation in the visible region, with a wavelength of
approximately $\lambda_{20th}=540$ nm. Moreover, harmonics as high as the
25th-26th can be also generated (although with small intensity). This
corresponds to the blue side of the visible spectrum, as $\lambda_{25th}=433$.
nm. The nonlinear signal in Fig. 5(c), therefore, spans the whole spectrum
between THz and the visible region, and could also be pushed into the near UV.
This suggests that artificial gauge fields in graphene could be used as a mean
to efficiently convert signals within these spectral region.

##  V Effect of the Pseudomagnetic Field on the Nonlinear Signal

We now briefly discuss what is the effect of the pseudomagnetic field on the
nonlinear signal. To do that, we focus our attention on the case of a
$\tau=50$ fs pulse, and repeat our simulations using different values of the
pseudomagnetic field. This, in practice, would correspond to a larger or
smaller bending radius (and, consequently, bending angle) for smaller or
greater values of the pseudomagnetic field, respectively (See Fig. 1).

We show the results of these simulations in Fig. 6, where the nonlinear signal
il plotted against three different values of the pseudomagnetic field, namely
$B=2$ T [Fig. 6 (a)], $B=5$ T [Fig. 6 (b)], and $B=15$ T [Fig. 6 (c)].

As it can be seen, for large values of $B$ [Fig. 6 (c)] with respect to those
employed in Fig. 5, half-integer harmonics centred around
$\omega/\omega_{1}\simeq 2-4$ and $\omega/\omega_{1}\simeq 6-8$ distinctively
appear. For $B=2$ T, on the other hand, a very efficient transfer of energy
between the impinging field at $\omega_{L}=\omega_{1}$ and its 20th harmonics
takes place, which means that in the presence of small pseudomagnetic fields
(i.e., large bending angles), graphene behaves as a very efficient frequency
converter between the impinging field oscillating at $\omega_{L}=174$ THz and
its 20th harmonics, that sits well-within the visible region, at about
$\lambda=540$ nm.

Figure 6: Nonlinear signal corresponding to an impinging electromagnetic pulse
of duration $\tau=50$ fs, for different values of the pseudomagnetic field,
namely (a) $B=2$ T, (b) $B=5$ T, and (c) $B=15$ T. For low pseudomagnetic
fields, we observe an almost resonant transfer of energy between the impinging
field as its 20th harmonics [panel (a)], or an almost equal distribution of
energy of the impinging pulse across its 10th-15th harmonics [panel(b)]. For
high values of the pseudomagnetic field [panel (c)], instead, we observe the
appearance of half-integer harmonics centreda round $\omega/\omega_{1}\simeq
2-3$ and $\omega/\omega_{1}\simeq 6-8$. For these plots,
$\omega_{L}=\omega_{1}=174$ THz, $E_{0}=10^{7}$ V/M (correspondent to a
fluence $F=6.64$ $m\,J/m^{2}$), and $v_{F}=c/300$ m/s have been used, to
reproduce results compatible with Fig. 5.

##  VI Conclusion and Outlook

In conclusion, we have investigated the nonlinear response of a bent graphene
flake to an impinging, linearly polarised, time-dependent electromagnetic
field. We have shown, that the AGF induced by bending, and the consequent
appearance of Landau levels significantly modify the harmonic signal generated
by graphene. in particular, we have shown, that for the case of an impinging
field resonant with the transition $|0\rangle\leftrightarrow|\pm 1\rangle$,
i.e., with the Landau ladder in the vicinity of the Dirac point, the nonlinear
signal shows a spectrum containing integer, as well as non-integer harmonics
of the pulse carrier frequency. In particular, for ultrashort pulses, the
nonlinear signal contains integer harmonics $\omega=n\omega_{L}$ as well as
half-integer harmonics $\omega=(n+1/2)\omega_{L}$, while as the pulse gets
longer, a more complicated scenario arises, with a nonlinear signal containing
up to the 26th harmonics in its spectrum. Lastly, we have shown that the
pseudomagnetic field induced by the AGF breaks the central symmetry typical of
graphene, thus allowing the emission of even harmonics of the impinging pulse
carrier frequency. The magnitude of the pseudomagnetic field can be controlled
by changing the radius of curvature (or, equivalently, the bending angle) of
the graphene flake. By doing so, we observe, for “small" magnetic fields, the
almost resonant transfer of energy between the impinging electromagnetic pulse
at $\omega=\omega_{L}$ and its 20th harmonics. For an impinging pulse centred
around $174$ THz, as in our case, the 20th harmonics lies well-within the
visible spectrum, around $\lambda=540$ nm. This effect, controllable by
controlling the magnitude of the pseudomagnetic field, could potentially pave
the way for novel graphene-based photonic devices such as THz-to-visible
frequency converters, and frequency generators.

In this manuscript, we have focused our attention on the role AGFs have on the
nonlinear optical response of graphene to an ultrashort pulse impinging upon
it at zero temperature. However, our model can be readily extended to account
for a finite temperature, using the more standard Dirac-Bloch model, which
will give us the possibility to include other competing effects, such as
dephasing, relaxation dynamics, and Coulomb interactions.

In this work, we have not explored the role of the valley degree of freedom on
the nonlinear response of strained graphene. It is known, however, that
introducing such a strain can break the valley symmetry extra1 ; extra2 . This
could then possibly lead to a valley-contrasting nonlinear optical response,
which could have significant impact on valleytronic applications. To account
for this, one would need to explicitly take spin-orbit coupling into account,
such that the combination of AGF and spin-orbit coupling will lift the
pseudospin degeneracy in both bulk and edge states, leading to time reversal
symmetry breaking extra3  and, ultimately, to a valley-dependent nonlinear
response.

The interplay between AGFs and spin-orbit coupling and their effect on the
nonlinear response of strained graphene, together with the extension of our
analysis to other 2D materials, such as transition metal dichalcogenides
(where we can directly probe the effect of Coulomb interactions on their
nonlinear response in the presence of AGFs), and the extension of our
formalism to 3D Dirac materials, such as topological semimetals extra4 ;
extra5 ; extra6 , will be the subject of our next research.

###### Acknowledgements.

M.O. acknowledges the support of the Academy of Finland Flagship Programme,
Photonics research and Innovation (PREIN), decision 320165.

## Data Availability

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

## References

## References

  * (1) J.D. Jackson, “Classical Electrodynamics", 3rd Edition (Wiley, New York, 1998).
  * (2) M. D. Schwartz, “Quantum Field Theory and the Standard Model" (Cambridge university Press, Cambridge, 2013).
  * (3) Y. Aharonov, and D. Bohm, “Significance of Electromagnetic Potentials in Quantum Theory", Phys. Rev. 115, 485 (1959).
  * (4) T. T. Wu, and C. N. Yang, “Concept of nonintegrable phase factors and global formulation of gauge fields", Phys. Rev. D 12 ,3845 (1975).
  * (5) B. Simon, “Holonomy, the quantum adiabatic theorem, and berry’s phase", Phys. Rev. Lett. 51, 2167 (1983)
  * (6) M. Nakahara, “Geometry, topology and physics", (CRC Press, Boca Raton, 2003).
  * (7) M. V. Berry, "Quantal Phase Factors Accompanying Adiabatic Changes". Proc. R. Soc. Lond. A. 392, 45 (1984).
  * (8) S. Pancharatnam, “Generalized theory of interference and its applications", in “Proceedings of the Indian Academy of Sciences-Section A", Vol. 44, pp.398-417 (Springer, Berlin, 1956).
  * (9) A. A. Abrikosov, “The magnetic properties of superconducting alloys", J. Phys. Chem. Solids. 2, 199 (1957).
  * (10) D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, “Quantized hall conductance in a two-dimensional periodic potential", Phys. Rev. Lett. 49, 405 (1982).
  * (11) M. Kohmoto, “Zero modes and the quantized hall conductance of the two-dimensional lattice in a magnetic field",Phys. Rev. B 39, 11943 (1989).
  * (12) E. J. Yarmchuk, M. J. V. Gordon, and R. E. Packard, “Observation of Stationary Vortex Arrays in Rotating Superfluid Helium", Phys. Rev. Lett. 43, 214 (1979).
  * (13) M. Ornigotti, G. Della Valle, D. Gatti, and S. Longhi, “Topological suppression of optical tunneling in a twisted annular fiber", Phys. Rev. A 76, 023833 (2007).
  * (14) S. D. Huber, “Topological Mechanics", Nat. Phys. 12, 621 (2016).
  * (15) B. A. Bernevig, and T. L. Hughes, “Topological insulators and topological superconductors", (Princeton University Press,Princeton,, 2013).
  * (16) N. Goldman, J. C. Budich, and P. Zoller, “Topological quantum matter with ultracold gases in optical lattices", Nat. Phys. 12, 639 (2016).
  * (17) J. Dalibard, F. Gerbier, G. Juzeliun̄as, and P. Öhberg, “Colloquium: Artificial gauge potentials for neutral atoms", Rev. Mod. Phys. 83, 1523 (2011).
  * (18) N. Goldman, G. Juzeliun̄as, P.Öhberg, and I. B. Spielman, “Light-induced gauge fields for ultracold atoms", Rep. Prog. Phys. 77, 126401 (2014).
  * (19) J. Dalibard, “Introduction to the physics of artificial gauge fields", arXiv:1504.05520 (2015).
  * (20) L. Lu, J. D. Joannopoulos, and M. Soljačić,“Topological Photonics", Nat. Photon. 8, 821 (2014).
  * (21) T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, “Topological Photonics", Rev. Mod. Phys. 91, 015006 (2019).
  * (22) H. J. Rothe, “Lattice gauge theories", 4th Edition (World Scientific, Singapore, 2012).
  * (23) M. C. Rechtsmann, J. Zeuner, Y. Plotnik, Y. Lumer, D. Podolsky, F. Dreisow, S. Nolte, M. Segev, and A. Szameit, “Photonic Floquet topological insulator", Nature 496, 196 (2013).
  * (24) L. J. Maczewsky, M. Heinrich, M. kremer, S. K. Ivanov, M. Ehrhardt, F. Martinez, Y. V. Kartashov, V. V. Konotp, L. Torner, D. bauer, and A. Szmaiet, “Nonlinearity-induced photonic topological insulator", Science 370, 701 (2020).
  * (25) A. B. Khanikaev, S. H. Mousavi, W.-K. Tse, M. Kargarian, A. H. MacDonald, and G. Shvets, “Photonic topological insulators", Nat. Mater. 12, 233 (2013).
  * (26) M. Kremer, L. Teuber, A. Szameit, and S. Scheel, “Optimal design strategy for non-Abelian geometric phases using Abelian gauge fields based on quantum metric", Phys. Rev. Research 1, 033117 (2019).
  * (27) C. Jörg, G. Queraltó, M. Kremer, G. Pelegrí, J. Schulz, A. Szameit, G. von Freymann, J. Mompart, and V. Ahufinger, “Artificial gauge field switching using orbital angular momentum modes in optical waveguides", Light: Science & Applications 9, 150 (2020).
  * (28) M. C. Rechtsman, J. M. Zeuner, A. Tünnermann, S. Nolte, M. Segev, and A. Szameit, “Strain-induced pseudomagnetic field and photonic Landau levels in dielectric structures", Nat. Phot. 7, 153 (2013).
  * (29) M. C. Rechtsman, Y. Lumer, Y. Plotnik, A. Perez-Leija, A. Szameit, and M. Segev, “Topological protection of photonic path entanglement", Optica 3, 925 (2016).
  * (30) J. Noh, W. A. Benalcazar, S. Huang, M. J. Collins, K. P. Chen, T. L. Hughes, and M. C. Rechtsman, “Topological protection of photonic mid-gap defect modes", Nat. Phot. 12, 408 (2018).
  * (31) M. Wang, C. Doyle, B. Bell, M. J. Collins, E. Magi, B. J. Eggleton, M. Segev, and A. Blanco-Redondo, “Topologically protected entangled photonic states", Nanophot. 8, 1327 (2019).
  * (32) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, “Electric Field Effect in Atomically Thin Carbon Films", Science, 306, 666 (2004).
  * (33) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, “Two-dimensional gas of massless Dirac fermions in graphene", Nature 438, 197 (2005).
  * (34) Yuanbo Zhang, Yan-Wen Tan, Horst L. Stormer, and P. Kim, “Experimental observation of the quantum Hall effect and Berry’s phase in graphene", Nature 438, 201 (2005).
  * (35) M. I. Katsnelson, “Zitterbewegung, chirality, and minimal conductivity in graphene,” Eur. Phys. J. B.51, 157 (2006).
  * (36) R.R.Nair, P.Blake, A.N.Grigorenko, K.S.Novoselov, T.J.Booth, T.Stauber, N.M.R.Peres, and A.K.Geim, “Fine Structure Constant Defines Visual Transparency of Graphene,” Science 320, 1308 (2008).
  * (37) M.A.H. Vozmediano, M.I. Katsnelson, and F. Guinea, “Gauge fields in graphene", Phys. Rep. 496 109 (2010).
  * (38) H. Suzuura, and T. Ando, “Phonons and electron-phonon scattering in carbon nanotubes", Phys. Rev. B 65, 235412 (2002).
  * (39) J. L. Ma$\tilde{\text{n}}$es, “Symmetry-based approach to electron-phonon interactions in graphene", Phys. rev. B 76, 045430 (2007).
  * (40) F. Guinea, A. K. Geim, M. I. Katsnelson, and K. S. Novoselov, “Generating quantizing pseudomagnetic fields by bending graphene ribbons", Phys. Rev. B 81, 035408 (2010).
  * (41) F. Guinea, M. I. Katsnelson, and M. A. H. Vozmediano, “Midgap states and charge inhomogeneities in corrugated graphene", Phys. Rev. B 77, 075422 (2008).
  * (42) A. Cortijo, and M. A. H. Vozmediano, “Electronic properties of curved graphene sheets", EPL 77, 47002 (2007).
  * (43) A. R. Wright, X. G. Xu, J. C. Cao, and Z. Chang, “Strong nonlinear response in the terahertz regime", Appl. Phys. Lett. 95, 072101 (2009).
  * (44) E. Hendry, P. J. Hale, J. Moger, and A. K. Savchenko, “Coherent nonlinear optical response of graphene", Phys. Rev. Lett. 105, 097401 (2010).
  * (45) H. A. Hafez, S. Kovalev, K.-J. Tielrooij, M. Bonn, M. Gensch, and D. Turchinovich, “The nonlinear optics of graphene: from saturable absorption to high-harmonic generation", Adv. opt. Mater. 8, 1900771 (2020).
  * (46) G. Demetriou, H. T. Bookey, F. Biancalana, E. Abraham, Y. Wang, W. Ji, and A. K. Kor, “Nonlinear optical properties of multilayer graphene in the infrared", Opt. Express 24, 13033 (2016).
  * (47) J. B. Khurgin, “Graphene - A rather ordinary nonlinear optical material", Appl. Phys. Lett. 104, 161116 (2014).
  * (48) R. W. Boyd, “Nonlinear Optics", 3rd Edition (Academic Press, Cambridge, MA, 2008).
  * (49) X. Yao, and A. Belyonin, “Nonlinear optics of graphene in a strong magnetic field", J. Phys.: Condens. Matter 25, 054203 (2013).
  * (50) Vinu Lukose, R. Shankar, and G. Baskaran, “Novel Electric Field Effects on Landau Levels in Graphene", Phys. Rev. Lett. 98, 116802 (2007).
  * (51) A. Lòpez, A. Di Teodoro, J. Schliemann, B. Berche, and B. Santos, “Laser-induced modulation of the Landau level structure in single-layer graphene", Phys. Rev. B 92, 235411 (2015).
  * (52) C. Ding, R. Yu, X. Hao, and D. Zhang, “Controllable population dynamics in Landau-quantized graphene", Sci. rep. 8, 1530 (2018).
  * (53) F. Guinea, M. I. Katsnelson, and A. K. Geim, “Energy gaps and a zero-field quantum Hall effect in graphene by strain engineering", nat. Phys. 6, 30 (2010).
  * (54) M. I. Katsnelson, “Graphene: Carbon in Two Dimensions", (Cambridge University Press, Cambridge, 2012).
  * (55) M. Oliva-Leyva, and G. G. Naumis, “Anisotropic AC conductivity of strained graphene", J. Phys.: Condes. Matter 26, 125302 (2014).
  * (56) K. L. Ishikawa, "Nonlinear optical response of graphene in time domain", Phys. Rev. B 82, 201402(R) (2010).
  * (57) M. Wegener, “Extreme Nonlinear optics" (Springer, Berlin, 2005).
  * (58) D. N. Carvalho, A. Marini, and F. Biancalana, “Dynamical centrosymmetry breaking - A novel mechanism for second-harmonic generation in graphene", Ann. Phys. 378, 24 (2017).
  * (59) D. N. Carvalho, F. Biancalana, and A. Marini, “Monolayer graphene can emit SHG waves", Opt. Data Process. Storage, 3, 47 (2017).
  * (60) D. S. L. Abergel, and V. I. Fal’ko, “Optical and magneto-optical far-infrared properties of bilayer graphene", Phys. Rev. B 75, 155430 (2007).
  * (61) R. Loudon, “The Quantum theory of Light", 3rd Edition (Oxford Science Publications, Oxford, 2000).
  * (62) Y. S. Ang, S. Sultan, and C. Zhang, “Nonlinear optical spectrum of bilayer graphene in the terahertz regime", Appl. Phys. Lett. 97, 243110 (2010).
  * (63) Y. Jiang, T. Low, K. Chang, M. I. Katsnelson, and F. Guinea, “Generation of pure valley current in graphene", Phys. Rev. Lett. 110, 046601 (2015).
  * (64) W.-Y. He, and L. He, “Coupled spin and pseudomagnetic field in graphene nanoribbons", Phys. Rev. B 28, 085411 (2013).
  * (65) K. J. A. Ooi, Y. S. Ang, Q. Zhan, D. T. H. Han, L. K. Ang, and C. K. Ong, “Nonlinear plasmonics of three-dimensional Dirac semimetals", APL Photonics 4, 034402 (2019).
  * (66) J. Lim, Y. S. Ang, F. J. Garcia de Abajo, I. Kaminer, L. K. Ang, and L. J. Wong, “Efficient generation of terahertz harmonics in three-dimensional Dirac semimetals", Phys. Rev. Research 2, 043252 (2020).
  * (67) J. Lim, K. J. A. Ooi, C. Zhang, L. K. Ang, and Y. S. Ang, “Broadband strong optical dichroism in topological Dirac semimetals with Fermi velocity anisotropy", Chinese Phys. B 29, 077802 (2020).