# Modulated harmonic wave in series connected discrete Josephson transmission
line: the discrete calculus approach

Eugene Kogan<EMAIL_ADDRESS>Department of Physics, Bar-Ilan
University, Ramat-Gan 52900, Israel  Max-Planck-Institut fur Physik komplexer
Systeme Dresden 01187, Germany

###### Abstract

We consider the modulated harmonic wave in the discrete series connected
Josephson transmission line (JTL). We formulate the approach to the modulation
problems for discrete wave equations based on discrete calculus. We check up
the approach by applying it to the Fermi-Pasta-Ulam-Tsingou type problem.
Applying the approach to the discrete JTL, we obtain the equation describing
the modulation amplitude, which turns out to be the defocusing nonlinear
Schrödinger (NLS) equation. We compare the profile of the single soliton
solution of the NLS with that of the soliton obtained in our previous
publication.

##  I Introduction

In our previous publication kogan  we considered propagation of kinks,
solitons and shocks along the Josephson transmission line (JTL). In the
present short note we want to consider a different kind of excitations –
modulated harmonic waves solitons . In literature, JTLs have been extensively
discussed in connection with travelling wave amplifies. The main interest was
in interaction of small number of harmonic waves (pump, signal and idler; see
e.g. brien  and references there). An exception to this statement is Ref.
grimso , where the wave packets are considered (implicitly) in the analysis,
but single-frequency excitations are used for visualizing results. Hence the
topic of modulated harmonic waves in JTLs is not very well explored. In the
present paper we will show that the modulation amplitude is described by the
defocusing nonlinear Schrödinger (NLS) equation kosevich ; kivshar ;
kevrekidis .

The rest of the paper is constructed as follows. In Section II we write down
equations describing JTL and present the definition of the modulated harmonic
wave. In Section III we present our approach based on discrete calculus (DC)
and show that in linear approximation, equation describing the modulation
amplitude turns out to be linear Schrödinger equation. Nonlinear problem is
considered in Section IV in the framework of the DC approach, and the equation
describing the modulation amplitude turns out to be the defocusing nonlinear
Schrödinger equation. In Section V we show that the DC approach, being applied
to Fermi-Pasta-Ulam-Tsingou type problem gives the results identical to the
known ones. In Section VI we compare dark solitons, known to exist for the
defocusing nonlinear Schrödinger equation, with the solitons in the JTL,
obtained by us previously. We conclude in Section VII. In the Appendix we
compare the results of the DC approach with those obtained in the framework of
the fourier integral representation of the solution.

##  II Discrete Josephson transmission line

Consider the model of JTL constructed from identical Josephson junctions (JJ)
and capacitors, which is shown in Figure 1. We take as dynamical variables the
phase differences (which we for brevity will call just phases) $\varphi_{n}$
across the JJ and the charges $q_{n}$ which have passed through the JJ. The
circuit equations are

$\displaystyle\frac{\hbar}{2e}\frac{d\varphi_{n}}{dt}$
$\displaystyle=\frac{1}{C}\left(q_{n+1}-2q_{n}+q_{n-1}\right)\,,$ (1a)
$\displaystyle\frac{dq_{n}}{dt}$ $\displaystyle=I_{c}\sin\varphi_{n}\,,$ (1b)

where $C$ is the capacitance, and $I_{c}$ is the critical current of the JJ.
Differentiating Eq. (1a) with respect to $t$ and substituting $dq_{n}/dt$ from
Eq. (1b), we obtain closed equation for $\varphi_{n}$ kogan

$\displaystyle\frac{d^{2}\varphi_{n}}{d\tau^{2}}=\sin\varphi_{n+1}-2\sin\varphi_{n}+\sin\varphi_{n-1}\,,$
(2)

where we have introduced the dimensionless time $\tau=t/\sqrt{L_{J}C}$, and
$L_{J}=\hbar/(2eI_{c})$.

Figure 1: Discrete JTL.

In the linear approximation, Eq. (2) takes the form

$\displaystyle\frac{d^{2}\varphi_{n}}{d\tau^{2}}=\varphi_{n+1}-2\varphi_{n}+\varphi_{n-1}\,.$
(3)

The dispersion law for the eigenmodes of Eq. (3) can be easily found kittel

$\displaystyle\omega(k)=2\left|\sin\left(\frac{k}{2}\right)\right|\,,\;\;k\in(-\pi,\pi]\,.$
(4)

We will consider a modulated harmonic wave, that is a carrier wave of wave
number $k_{0}$ (for the sake of definiteness we take $k_{0}>0$) and frequency
$\omega_{0}=2\sin(k_{0}/2)$, modulated by a waveform $a$, which varies slowly
in time and space compared to the variations of the carrier wave

$\displaystyle\varphi(n,\tau)=\psi(n,\tau)+\psi^{*}(n,\tau)=a_{n}b_{n}+a_{n}^{*}b_{n}^{*}\,,$
(5)

where

$\displaystyle
a_{n}=a(n,\tau)\,,\;\;\;\;b_{n}=e^{i(k_{0}n-\omega_{0}\tau)}\,,$ (6)

##  III Linear approximation

Let us start from analysing Eq. (3). Because the equation is linear, we obtain
decoupled equations for $\psi$ and $\psi^{*}$, and equation for $\psi$ becomes

$\displaystyle\frac{d^{2}\psi_{n}}{d\tau^{2}}=\psi_{n+1}-2\psi_{n}+\psi_{n-1}\,.$
(7)

Substituting $\psi_{n}=a_{n}b_{n}$ into (7) we obtain

$\displaystyle
a_{n}\frac{d^{2}b_{n}}{d\tau^{2}}+2\frac{da_{n}}{d\tau}\frac{db_{n}}{d\tau}+\frac{d^{2}a_{n}}{d\tau^{2}}b_{n}$
$\displaystyle=a_{n+1}b_{n+1}-2a_{n}b_{n}+a_{n-1}b_{n-1}\,.$ (8)

We first focus on the r.h.s. of (III) and use the elementary formula for the
discrete second derivative of the product of two functions

$\displaystyle a_{n+1}b_{n+1}-2a_{n}b_{n}+a_{n-1}b_{n-1}$
$\displaystyle=a_{n}\left(b_{n+1}-2b_{n}+b_{n-1}\right)+\left(a_{n+1}-a_{n}\right)\left(b_{n+1}-b_{n}\right)$
$\displaystyle+\left(a_{n}-a_{n-1}\right)\left(b_{n}-b_{n-1}\right)+\left(a_{n+1}-2a_{n}+a_{n-1}\right)b_{n}\,.$
(9)

Then the first term of the r.h.s. of (III) cancels against the
$a_{n}d^{2}b_{n}/d\tau^{2}$-term on the l.h.s. of (III), because $b_{n}$ is by
itself the solution of (3), and we arrive to the equation

$\displaystyle
2\frac{da_{n}}{d\tau}\frac{db_{n}}{d\tau}+\frac{d^{2}a_{n}}{d\tau^{2}}b_{n}=\left(a_{n+1}-a_{n}\right)\left(b_{n+1}-b_{n}\right)$
$\displaystyle+\left(a_{n}-a_{n-1}\right)\left(b_{n}-b_{n-1}\right)+\left(a_{n+1}-2a_{n}+a_{n-1}\right)b_{n},$

or, more explicitly (taking into account (4)),

$\displaystyle-4i\sin\left(\frac{k_{0}}{2}\right)\frac{da_{n}}{d\tau}$
$\displaystyle+$ $\displaystyle\frac{d^{2}a_{n}}{d\tau^{2}}=i\sin
k_{0}\left(a_{n+1}-a_{n-1}\right)$ (11) $\displaystyle+$ $\displaystyle\cos
k_{0}\left(a_{n+1}-2a_{n}+a_{n-1}\right).$

In the lowest order (with respect to the ratio of the effective modulation
wave vector and $k_{0}$) approximation we should equate the first terms in the
l.h.s. and r.h.s. of (11) thus obtaining

$\displaystyle\frac{da_{n}}{d\tau}=-\frac{1}{2}\cos\left(\frac{k_{0}}{2}\right)\left(a_{n+1}-a_{n-1}\right)$
(12)

and, hence,

$\displaystyle\frac{d^{2}a_{n}}{d\tau^{2}}=\frac{1}{4}\cos^{2}\left(\frac{k_{0}}{2}\right)\left(a_{n+2}-2a_{n}+a_{n-2}\right)\,.$
(13)

Substituting (13) into (11) we obtain a discrete linear Schrödinger equation
(with an additional convection term that can be removed by a change of
reference frame)

$\displaystyle-4i\sin\left(\frac{k_{0}}{2}\right)\frac{da_{n}}{d\tau}=i\sin
k_{0}\left(a_{n+1}-a_{n-1}\right)$ $\displaystyle+\cos
k_{0}\left(a_{n+1}-2a_{n}+a_{n-1}\right)$
$\displaystyle-\frac{1}{4}\cos^{2}\left(\frac{k_{0}}{2}\right)\left(a_{n+2}-2a_{n}+a_{n-2}\right)\,.$
(14)

Substituting continuous variable $x$ for the discrete variable $n$ as the
argument of $a$ and expanding around $x=n$, we can present the quantities in
the parentheses of (III) as

$\displaystyle a_{n+1}-a_{n-1}$ $\displaystyle=2\frac{\partial a}{\partial
x}\,,$ (15a) $\displaystyle a_{n+1}-2a_{n}+a_{n-1}$
$\displaystyle=\frac{\partial^{2}a}{\partial x^{2}}\,,\,$ (15b) $\displaystyle
a_{n+2}-2a_{n}+a_{n-2}$ $\displaystyle=4\frac{\partial^{2}a}{\partial
x^{2}}\,,$ (15c)

and after simple algebra the equation takes the form of the linear Schrödinger
equation

$\displaystyle i\left(\frac{\partial a}{\partial\tau}+v_{g}\frac{\partial
a}{\partial x}\right)-D\frac{\partial^{2}a}{\partial x^{2}}=0\,,$ (16)

where the group velocity $v_{g}$ and the coefficient $D$ are

$\displaystyle
v_{g}=\cos\left(\frac{k_{0}}{2}\right)\,,\;\;\;\;\;D=\frac{1}{4}\sin\left(\frac{k_{0}}{2}\right)\,.$
(17)

When one looks at (16) a natural question appears: If the approximation of the
underlying chain of coupled junctions by the quasi-continuum limit is
extended, what will be the next-order terms in the resulting equation? We can
answer this question by inspection of (15). In case of the extension mentioned
above, in the derived linear Schrödinger equation there would appear
additional (quasi) drift term with the third derivative with respect to $x$,
and additional term with the forth derivative with respect to $x$.

##  IV The nonlinear Schrödinger equation

Now let us return to Eq. (2). Presenting $\varphi$ as in (5) and expanding the
sine in series we get

$\displaystyle\sin\varphi=\sum_{m=0}^{\infty}\frac{(-1)^{m}}{(2m+1)!}(\psi+\psi^{*})^{2m+1}\,.$
(18)

Keeping only the first harmonics (the rotating wave approximation), we will
present $\sin\varphi$ as

$\displaystyle\sin\varphi$ $\displaystyle=$
$\displaystyle\left(\psi+\psi^{*}\right)\sum_{m=0}^{\infty}\frac{(-1)^{m}}{m!(m+1)!}|\psi|^{2m}$
(19) $\displaystyle=$
$\displaystyle\left(\psi+\psi^{*}\right)\frac{J_{1}(2|\psi|)}{|\psi|}\,,$

where $J_{1}$ is the Bessel function. Thus we again obtain decoupled equations
for $\psi$ and $\psi^{*}$, and equation for $\psi$ becomes (compare with (7))

$\displaystyle\frac{d^{2}\psi_{n}}{d\tau^{2}}=\frac{J_{1}(2|\psi_{n+1}|)}{|\psi_{n+1}|}\psi_{n+1}-2\frac{J_{1}(2|\psi_{n}|)}{|\psi_{n}|}\psi_{n}$
$\displaystyle+\frac{J_{1}(2|\psi_{n-1}|)}{|\psi_{n-1}|}\psi_{n-1}\,.$ (20)

Note that in the lowest nontrivial order (IV) takes the form

$\displaystyle\frac{d^{2}\psi_{n}}{d\tau^{2}}=\left(1-\frac{|\psi_{n+1}|^{2}}{2}\right)\psi_{n+1}-2\left(1-\frac{|\psi_{n}|^{2}}{2}\right)\psi_{n}$
$\displaystyle+\left(1-\frac{|\psi_{n-1}|^{2}}{2}\right)\psi_{n-1}\,.$ (21)

Using (5) we obtain

$\displaystyle
2\frac{da_{n}}{d\tau}\frac{db_{n}}{d\tau}+\frac{d^{2}a_{n}}{d\tau^{2}}b_{n}=\left(a_{n+1}-a_{n}\right)\left(b_{n+1}-b_{n}\right)$
$\displaystyle+\left(a_{n}-a_{n-1}\right)\left(b_{n}-b_{n-1}\right)+\left(a_{n+1}-2a_{n}+a_{n-1}\right)b_{n}$
$\displaystyle-d.s.d.\,,$ (22)

where the last term in the r.h.s. of (IV) is the discrete second derivative
(d.s.d.) of the product of three quantities

$\displaystyle
d.s.d.=a_{n+1}g_{n+1}b_{n+1}-2a_{n}g_{n}b_{n}+a_{n-1}g_{n-1}b_{n-1},$ (23)

and

$\displaystyle g_{n}\equiv
g\left(|a_{n}|^{2}\right)=1-\frac{J_{1}(2|a_{n}|)}{|a_{n}|}\,.$ (24)

The first three terms in the r.h.s. of (IV) are identical to the r.h.s. of
(III). The forth term can be treated in a simpler way. In fact, the difficulty
of treating the r.h.s. of (III) was connected with the fact that in the
expansion with respect to the ratio of the effective modulation wave vector
and $k_{0}$, the lowest order term was canceled with the appropriate term in
the l.h.s. Thus we have to take into account the next order terms. There is no
such cancellation for the d.s.d. term in the r.h.s. of (IV), so the term can
be considered in the lowest order approximation. Among the three quantities
$a_{n}$, $g_{n}$ and $b_{n}$, two ($a_{n}$ and $g_{n}$) change slowly with
$n$, and the third quantity ($b_{n}$) changes fast. This is why in the r.h.s.
of (23) we can ignore the difference between $a_{n}$, $a_{n-1}$ and $a_{n+1}$
(and between $g_{n}$, $g_{n-1}$ and $g_{n+1}$) and present the equation as

$\displaystyle d.s.d.=a_{n}g_{n}\left(b_{n+1}-2b_{n}+b_{n-1}\right)\,.$ (25)

Proceeding as in Section III, we obtain, instead of (16), the defocusing NLS
karpman ; kosevich ; kivshar ; kevrekidis

$\displaystyle i\left(\frac{\partial a}{\partial\tau}+v_{g}\frac{\partial
a}{\partial x}\right)-D\frac{\partial^{2}a}{\partial
x^{2}}+4Dg\left(|a|^{2}\right)a=0\,.$ (26)

In the lowest nontrivial order with respect to $|a|^{2}$, Eq. (26) is reduced
to

$\displaystyle i\left(\frac{\partial a}{\partial\tau}+v_{g}\frac{\partial
a}{\partial x}\right)-D\frac{\partial^{2}a}{\partial x^{2}}+2D|a|^{2}a=0\,.$
(27)

As an additional support for the validity of Eq. (27) we will ”rederive” it in
the Appendix.

##  V The Fermi-Pasta-Ulam-Tsingou problem

To check up our DC method, it would be appropriate to apply it to the Fermi-
Pasta-Ulam-Tsingou (FPUT) problem gallavotti . The FPUT analog of (2) would be
james ; james2

$\displaystyle\frac{d^{2}\varphi_{n}}{d\tau^{2}}=\sin\left(\varphi_{n+1}-\varphi_{n}\right)+\sin\left(\varphi_{n-1}-\varphi_{n}\right)\,.$
(28)

In the rotating wave approximation we again have the decoupling of the
equations for $\psi$ and $\psi^{*}$. Expanding sine in Taylor series and
keeping the two lowest order nonzero terms, we obtain equation for $\psi$ in
the form

$\displaystyle\frac{d^{2}\psi_{n}}{d\tau^{2}}=\psi_{n+1}-\psi_{n}+\psi_{n-1}-\psi_{n}$
$\displaystyle-\frac{1}{2}\left[\left|\psi_{n+1}^{2}\right|\psi_{n+1}-\psi_{n+1}^{2}\psi_{n}^{*}-2\left|\psi_{n}^{2}\right|\psi_{n}\right.$
$\displaystyle+\psi_{n}^{2}\psi_{n+1}^{*}-2\left|\psi_{n+1}^{2}\right|\psi_{n}+2\left|\psi_{n}^{2}\right|\psi_{n+1}$
$\displaystyle+\left|\psi_{n-1}^{2}\right|\psi_{n-1}-\psi_{n-1}^{2}\psi_{n}^{*}$
$\displaystyle\left.+\psi_{n}^{2}\psi_{n-1}^{*}-2\left|\psi_{n-1}^{2}\right|\psi_{n}+2\left|\psi_{n}^{2}\right|\psi_{n-1}\right]\,.$
(29)

Following the example of Section IV we again obtain the defocusing NLS

$\displaystyle i\left(\frac{\partial a}{\partial\tau}+v_{g}\frac{\partial
a}{\partial x}\right)-D\frac{\partial^{2}a}{\partial
x^{2}}+D^{\prime}|a|^{2}a=0\,,$ (30)

only this time

$\displaystyle
D^{\prime}=\frac{1}{8b_{n}\sin\left(\frac{k_{0}}{2}\right)}\left(3b_{n+1}+3b_{n-1}-b_{n+1}^{2}b_{n}^{*}\right.$
$\displaystyle\left.-b_{n-1}^{2}b_{n}^{*}-6b_{n}+b_{n}^{2}b_{n+1}^{*}+b_{n}^{2}b_{n-1}^{*}\right)\,.$
(31)

After simple algebra we obtain

$\displaystyle D^{\prime}=2\sin^{3}\left(\frac{k_{0}}{2}\right)\,.$ (32)

Equation (30) exactly coincides with the result of Ref. fly , in the
appropriate particular case. (In that Reference, equation more general than
(28) is considered.)

##  VI Dark solitons

The defocusing NLS equation has an interesting type of solutions, called dark
solitons zakharov ; veks ; nateo ; triki . Using the opportunity, we would
like to present here the derivation of these solutions, borrowed to some
extent from the book kosevich  (which, to the best of our knowledge, was never
translated into English).

Making the transformation

$\displaystyle\xi=\left(x-v_{g}\tau\right),\;\;\;\tau^{\prime}=D\tau,$ (33)

we can present (27) as

$\displaystyle i\frac{\partial
a}{\partial\tau^{\prime}}-\frac{\partial^{2}a}{\partial\xi^{2}}+2|a|^{2}a=0\,.$
(34)

Looking for the solutions of (34) in the form

$\displaystyle a=Ae^{i\Phi}e^{i\Omega\tau^{\prime}},$ (35)

where $\Omega$ is a constant, and the amplitude $A$ and the phase $\Phi$ of
$a$ are some functions of $\xi$, we obtain for those functions the system of
equations

$\displaystyle\frac{d}{d\xi}\left(A^{2}\frac{d\Phi}{d\xi}\right)$
$\displaystyle=0\,,$ (36a)
$\displaystyle\frac{d^{2}A}{d\xi^{2}}-A\left(\frac{d\Phi}{d\xi}\right)^{2}-2A^{3}+\Omega
A$ $\displaystyle=0\,.$ (36b)

From (36a) we obtain

$\displaystyle\frac{d\Phi}{d\xi}=\frac{c}{A^{2}}\,,$ (37)

where $c$ is the integration constant. Substituting this expression into (36b)
we obtain

$\displaystyle\frac{d^{2}A}{d\xi^{2}}-\frac{c^{2}}{A^{3}}-2A^{3}+\Omega A$
$\displaystyle=0\,.$ (38)

Integrating, we obtain the differential equation of the first order

$\displaystyle\left(\frac{dA}{d\xi}\right)^{2}+\frac{c^{2}}{A^{2}}-A^{4}+\Omega
A^{2}=E\,,$ (39)

where $E$ is another integration constant. We are looking for the solutions
satisfying the conditions

$\displaystyle\lim_{\xi\to\pm\infty}A=A_{1}.$ (40)

Expressing the constants of integration through $\Omega$ and $A_{1}$

$\displaystyle c^{2}=A_{1}^{4}\left(\Omega-2A_{1}^{2}\right)\,,\hskip
28.45274ptE=A_{1}^{2}\left(2\Omega-3A_{1}^{2}\right)\,,$ (41)

we may present (39) as

$\displaystyle\left(\frac{d\rho}{d\xi}\right)^{2}=4\left(\rho-\rho_{1}\right)^{2}\left(\rho+2\rho_{1}-\Omega\right),$
(42)

where we have introduced $\rho=A^{2}$. Equation (42) can be easily integrated,
and we obtain

$\displaystyle\rho=\rho_{1}-\frac{\kappa^{2}}{\cosh^{2}\kappa\xi},$ (43)

where $\kappa^{2}=3\rho_{1}-\Omega$. It is convenient to present Eq. $\kappa$
as

$\displaystyle\kappa^{2}=A_{1}^{2}-A_{0}^{2},$ (44)

where $A_{0}$ is the minimal value of $A$ (achieved at $\xi=0$). The function
$A(\xi)$ defined by Eq. (43) is presented in Figure 2.

Figure 2: Dark soliton as given by Eq. (43. We have chosen $A_{0}=.5A_{1}$.

Using immobile soliton solution and the property of Galilean invariance of
(34), it is easy to obtain a solution describing moving soliton. In general,
if $a(\xi,\tau^{\prime})$ is a solution of (34), so is
$e^{i(2V\xi+V^{2}\tau^{\prime})}a(\xi-V\tau^{\prime},\tau^{\prime})$ for
arbitrary $V$.

It would be appropriate to compare the soliton given by Eq. (43) with the
solitons in the JTL presented in our previous publication kogan . Those
solitons were characterized by the Josephson phase, asymptotically constant at
both ends of the line

$\displaystyle\lim_{x\to\pm\infty}\varphi=\varphi_{1}.$ (45)

For the solitons described by (43), the phase at both sides of the
transmission line asymptotically coincides with that of a (high frequency)
harmonic wave. In general, while carrier wave is all important in the present
paper, there was no such wave whatsoever in our previous publication kogan .

However, it is worth to compare Figure 2 with Figure 3 borrowed from our
previous publication, and showing the soliton profile calculated there. Though
we didn’t use there the term dark soliton previously kogan , we see that the
curves are very much similar: in both cases the amplitude goes to a constant
value when $x$ ($\xi$) goes to infinities and decreases in between.

Figure 3: The soliton profile calculated in Ref. kogan . ($\Lambda$ is the
period of the JTL, which in the present publication was chosen to be equal to
one.)

##  VII Conclusions

We consider the modulated harmonic wave in the discrete series connected
Josephson transmission line (JTL). We formulate the approach to the modulation
problems for discrete wave equations based on discrete calculus. We check up
the approach by applying it to the Fermi-Pasta-Ulam-Tsingou type problem.
Applying the approach to the discrete JTL, we obtain the equation describing
the modulation amplitude, which turns out to be the defocusing nonlinear
Schrödinger (NLS) equation. The NLS, being the normal form for envelope waves,
represents a universal model at the root of an extremely wide range of
physical and other natural phenomena and applications. Furthermore, due to its
rich and complex phenomenology, the NLS is also a paradigm for nonlinear
spatio-temporal dynamics and is at the forefront of intense and challenging
mathematical research kevrekidis . We presented a new derivation of the single
soliton solution of the NLS and compared its profiles with the profile of the
soliton, obtained in our previous publication.

###### Acknowledgements.

We are grateful to M. Goldstein, G. James and B. Malomed for their insightful
comments. We are also very grateful to the anonymous Referees. The work on the
problem was initiated by my participation in the workshop Coherent Structures:
Current Developments and Future Challenges 4-8 July @ Oort. I would like to
express my gratitude to the Lorentz Center for the hospitality and for the
stimulating atmosphere.

##  Appendix A The approach based on Fourier integral

To additionally check up the discrete calculus approach we present here an
alternative derivation of some of the results obtained in the main body of the
paper, based on Fourier integral representation of the solution.

First consider the linear case. Equation (16) can be put in a more general
context solitons . Let us present the solution $\psi(n,\tau)$ of a general
linear equation as the Fourier integral

$\displaystyle\psi(n,\tau)=\int
dk\psi_{k}\exp\left\\{ikn-i\omega(k)\tau\right\\}\,.$ (46)

From (46) follows

$\displaystyle a(x,\tau)=\int
dk\psi_{k}\exp\left\\{i(k-k_{0})x-i\left[\omega(k)-\omega_{0}\right]\tau\right\\}\,.$

Differentiating with respect to $x$ we obtain

$\displaystyle i\frac{\partial a}{\partial x}=-\int dk(k-k_{0})\psi_{k}$
$\displaystyle\cdot\exp\left\\{i(k-k_{0})x-i\left[\omega(k)-\omega_{0}\right]\tau\right\\}.$
(48)

Differentiating with respect to $\tau$ we obtain

$\displaystyle i\frac{\partial a}{\partial\tau}=\int
dk\left[\omega(k)-\omega_{0}\right]\psi_{k}$
$\displaystyle\cdot\exp\left\\{i(k-k_{0})x-i\left[\omega(k)-\omega_{0})\right]\tau\right\\}.$
(49)

Expanding $\omega(k)$ with respect to $k$ up to the second order we get

$\displaystyle\omega(k)-\omega_{0}=\left(\frac{d\omega(k)}{dk}\right)_{0}(k-k_{0})$
$\displaystyle+\frac{1}{2}\left(\frac{d^{2}\omega(k)}{dk^{2}}\right)_{0}(k-k_{0})^{2}$
(50)

(the derivatives are calculated at $k=k_{0}$). Comparing (A) with (A) we
obtain the equation for the modulation amplitude

$\displaystyle i\left[\frac{\partial
a}{\partial\tau}+\left(\frac{d\omega(k)}{dk}\right)_{0}\frac{\partial
a}{\partial
x}\right]+\frac{1}{2}\left(\frac{d^{2}\omega(k)}{dk^{2}}\right)_{0}\frac{\partial^{2}a}{\partial
x^{2}}=0\,.$ (51)

For the dispersion law (4), Eq. (51) exactly coincides with (16).

The mnemonic rule for obtaining the equation for amplitude, modulating the
harmonic wave $\exp\left\\{ik_{0}n-i\omega(k_{0})\tau\right\\}$, can be
formulated as follows. Take the expansion of $\omega(k)$ with respect to $k$
up to the second order (A), replace $k-k_{0}$ by a spacial operator
$-i\partial/\partial x$, and $\omega-\omega_{0}$ by a temporal operator
$i\partial/\partial\tau$, and let (A) operate on the complex amplitude
function $a$.

If we present the complex amplitude $a$ as

$\displaystyle a=|a|e^{i\theta}\,,$ (52)

the complex equation (51) can be presented as two real equations

$\displaystyle\frac{\partial|a|^{2}}{\partial\tau}+\left(\frac{d\omega(k)}{dk}\right)_{0}\frac{\partial|a|^{2}}{\partial
x}$
$\displaystyle+\left(\frac{d^{2}\omega(k)}{dk^{2}}\right)_{0}\frac{\partial}{\partial
x}\left(|a|^{2}\frac{\partial\theta}{\partial x}\right)=0\,,$ (53a)
$\displaystyle\frac{\partial\theta}{\partial\tau}+\left(\frac{d\omega(k)}{dk}\right)_{0}\frac{\partial\theta}{\partial
x}$
$\displaystyle+\frac{1}{2}\left(\frac{d^{2}\omega(k)}{dk^{2}}\right)_{0}\left[\left(\frac{\partial\theta}{\partial
x}\right)^{2}-\frac{\partial^{2}|a|}{\partial x^{2}}\right]=0\,.$ (53b)

It is interesting to compare Eq. (53) with the equations borrowed from
geometric optics whitham . These equations describe a slowly varying wavetrain
by equations determining the propagation of wave number and frequency

$\displaystyle\frac{\partial|a|^{2}}{\partial\tau}+\frac{d\omega(k)}{dk}\frac{\partial|a|^{2}}{\partial
x}+\frac{d^{2}\omega(k)}{dk^{2}}|a|^{2}\frac{\partial k}{\partial x}$
$\displaystyle=0\,,$ (54a) $\displaystyle\frac{\partial
k}{\partial\tau}+\frac{d\omega(k)}{dk}\frac{\partial k}{\partial x}$
$\displaystyle=0\,.$ (54b)

It is convenient to write down Eq. (54) in a more explicit form by
substituting $k=k_{0}+\partial\theta/\partial x$

$\displaystyle\frac{\partial|a|^{2}}{\partial\tau}$
$\displaystyle+\left(\frac{d\omega(k)}{dk}\right)_{0}\frac{\partial|a|^{2}}{\partial
x}+\left(\frac{d^{2}\omega(k)}{dk^{2}}\right)_{0}\frac{\partial\theta}{\partial
x}\frac{\partial|a|^{2}}{\partial x}$
$\displaystyle+\left(\frac{d^{2}\omega(k)}{dk^{2}}\right)|a|^{2}\frac{\partial^{2}\theta}{\partial
x^{2}}=0\,,$ (55a) $\displaystyle\frac{\partial^{2}\theta}{\partial
x\partial\tau}$
$\displaystyle+\left(\frac{d\omega(k)}{dk}\right)_{0}\frac{\partial^{2}\theta}{\partial
x^{2}}+\left(\frac{d^{2}\omega(k)}{dk^{2}}\right)_{0}\frac{\partial\theta}{\partial
x}\frac{\partial^{2}\theta}{\partial x^{2}}=0\,.$ (55b)

Equation (55a) is identical to (53a). To compare (55b) with (53b), let us
differentiate the latter with respect to $x$. We obtain

$\displaystyle\frac{\partial\theta}{\partial
x\partial\tau}+\left(\frac{d\omega(k)}{dk}\right)_{0}\frac{\partial^{2}\theta}{\partial
x^{2}}+\left(\frac{d^{2}\omega(k)}{dk^{2}}\right)_{0}\frac{\partial\theta}{\partial
x}\frac{\partial^{2}\theta}{\partial x^{2}}$
$\displaystyle-\frac{1}{2}\left(\frac{d^{2}\omega(k)}{dk^{2}}\right)_{0}\frac{\partial^{3}|a|}{\partial
x^{3}}=0\,.$ (56)

In the r.h.s. of Eq. (A) there is an additional term in comparison with (55b),
but this term is probably negligible within the framework of the
approximations made.

Now we can turn to the nonlinear case. Equation (27) can be ”rederived”
following the pattern presented above. The solution of (IV) with constant
amplitude $a$ is

$\displaystyle\psi_{n}=ae^{i(kn-\omega t)}\,.$ (57)

Substituting (57) into (IV) we get the nonlinear dispersion law

$\displaystyle\omega(k;|a|^{2})=2\left(1-\frac{|a|^{2}}{4}\right)\left|\sin\left(\frac{k}{2}\right)\right|=\omega(k)+\omega_{2}(k)|a|^{2}\,,$

where $\omega_{2}=-|\sin(k/2)|/2$. Hence (A) should be modified to

$\displaystyle\omega-\omega_{0}=\left(\frac{d\omega(k;0)}{dk}\right)_{0}(k-k_{0})$
$\displaystyle+\frac{1}{2}\left(\frac{d^{2}\omega(k;0)}{dk^{2}}\right)_{0}(k-k_{0})^{2}+\omega_{2}(k)|a|^{2}\,.$
(59)

If we assume that in the nonlinear case, the formula (A) is modified to a
self-consistent equation

$\displaystyle a(x,\tau)=\int dk\psi_{k}$ (60) $\displaystyle\cdot$
$\displaystyle\exp\left\\{i(k-k_{0})x-i\left[\omega(k;|a|^{2})-\omega_{0}\right]\tau\right\\}\,,$

then (A) remains as it was in the linear case. Differentiating (60) with
respect to $\tau$ we obtain

$\displaystyle i\frac{\partial a}{\partial\tau}$ $\displaystyle=$
$\displaystyle\int
dk\psi_{k}\exp\left\\{i(k-k_{0})x-i\left[\omega(k;|a|^{2})-\omega_{0}\right]\tau\right\\}$
(61) $\displaystyle\cdot$
$\displaystyle\left[\left(\frac{d\omega(k)}{dk}\right)_{0}(k-k_{0})+\frac{1}{2}\left(\frac{d^{2}\omega(k)}{dk^{2}}\right)_{0}(k-k_{0})^{2}\right.$
$\displaystyle+$
$\displaystyle\left.\frac{\partial\omega(k;|a|^{2})}{\partial|a|^{2}}|a|^{2}\right]\,.$

Comparing (61) with (A) we get

$\displaystyle i\left[\frac{\partial
a}{\partial\tau}+\left(\frac{\partial\omega(k;a)}{\partial
k}\right)_{0}\frac{\partial a}{\partial x}\right]$ $\displaystyle+$
$\displaystyle\frac{1}{2}\left(\frac{\partial^{2}\omega(k;a)}{\partial
k^{2}}\right)_{0}\frac{\partial^{2}a}{\partial x^{2}}$
$\displaystyle-\frac{\partial\omega(k;|a|^{2})}{\partial|a|^{2}}|a|^{2}$
$\displaystyle=$ $\displaystyle 0\,,$ (62)

which, for the dispersion law (A), exactly coincides with (27). We hope that
the identity of the results obtained by two different approximate methods
gives us additional confidence in their validity.

## References

  * (1) E. Kogan, Phys. Status Solidi B, 2200160 (2022).
  * (2) M. Remoissenet, Waves Called Solitons: Concepts and Experiments, Third Edition, Corrected Second Printing, Springer-Verlag Berlin Heidelberg GmbH (2003).
  * (3) K. O’Brien, C. Macklin, I. Siddiqi, and X. Zhang, Phys. Rev. Lett. 113, 157001 (2014).
  * (4) A. L. Grimsmo and A. Blais, NPJ—Quantum Inf. 3:20 (2017).
  * (5) A. M. Kosevich, A. S. Kovalev, Introduction into Nonlinear Physical Mechanics, Naukova Dumka, Kiev (1989).
  * (6) Yu. S. Kivshar and B. Luther-Davies, Phys. Rep. 298, 81 (1998).
  * (7) P. G. Kevrekidis, D. J. Frantzeskakis, and R. Carretero-Gonzalez, The defocusing nonlinear Schrödinger equation: from dark solitons to vortices and vortex ring, Society for Industrial and Applied Mathematics (2015).
  * (8) C. Kittel, P. McEuen, and P. McEuen, Introduction to solid state physics (Vol. 8), New York: Wiley (1996).
  * (9) V. I. Karpman, Nonlinear Waves in Dispersive Media, Pergamon Press (1975)
  * (10) V. E. Zakharov and A. B. Shabat, Sov. Phys. JETP 37, 823 (1973).
  * (11) D. J. Frantzeskakis, Journal of Physics A: Mathematical and Theoretical 43, 213001 (2010).
  * (12) A. Biswas and D. Milovic, Communications in Nonlinear Science and Numerical Simulation 15, 1473, (2010).
  * (13) H. Triki, T. Hayat, O. M. Aldossary, and A. Biswas, Optics & Laser Technology, 44, 2223 (2012).
  * (14) G. Gallavotti, ed. The Fermi-Pasta-Ulam problem: a status report, Vol. 728. Springer, 2007.
  * (15) G. James, J. Nonlinear Sci. 13, 273 (2003).
  * (16) G. James, Phys. Rev. B 70, 014301 (2004).
  * (17) N. Flytzanis, St. Pnevmatikos, and M. Remoissenet, J. Phys. C: Solid State Phys. 18, 4603 (1985).
  * (18) G. B. Whitham, Linear and Nonlinear Waves, John Wiley & Sons Inc., New York (1999).


# Multiphoton excitation and high harmonic generation in rectangular graphene
quantum dot

A G Ghazaryan∗ and Kh V Sedrakian  Centre of Strong Fields Physics, Yerevan
State University, 1 A. Manukian, Yerevan 0025, Armenia
<EMAIL_ADDRESS>
###### Abstract

The multiphoton excitation and high harmonic generation (HHG) processes are
considered using the microscopic quantum theory of nonlinear interaction of
strong coherent electromagnetic (EM) radiation with rectangular graphene
quantum dot (RGQD). The dynamic Hartree-Fock approximation is developed for
the consideration of the quantum dot-laser field nonlinear interaction at the
nonadiabatic multiphoton excitation regime. The many-body Coulomb interaction
is described in the extended Hubbard approximation. By numerical results, we
show the significance of the RGQD lateral size, shape, and EM wavefield
orientation in RGQD of the zigzag edge compear to the armchair edge in the HHG
process allowing for increasing the cutoff photon energy and the quantum yield
of higher harmonics.

  *

##  1 Introduction

The quantum electrodynamic phenomena caused by the strong coherent EM
radiation, such as multiphoton HHG and related processes, through nonlinear
channels in two-dimensional (2D) atomic systems-nanostructures are of great
interest in lowenergy physics and nanooptoelectronics due to the unique
physical properties of such 2D nanosystems of atomic thickness [1, 2, 3, 4, 5,
6, 7, 8, 9]. It was been studied systematically the HHG in bulk crystals [10,
11, 12, 13, 14, 15, 16] and low-dimensional nanostructures such as graphene
and its derivatives [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32, 33, 34, 35, 36, 37, 38, 39], monolayer transition metal
dichalcogenides [40, 41, 42], hexagonal boron nitride [43], topological
insulator [44], [45], monolayers of black phosphorus [46], curved 2D hexagonal
nanostructures [47], solids [48], [49], as well as in other 2D systems [50,
51, 52]. 2D nanosystems enable the development of important technological
applications [53]. The quantum cascade laser is one such example [54], using
physical phenomena in 2D systems such as the quantum Hall effect [8].

As a nonlinear medium, the graphene quantum dot (GQD) [53], [55, 56, 57, 58],
and such as graphene nanoribbon [1], [59, 60, 61, 62, 63, 64, 65, 66, 67, 68,
69, 70], are of particular interest. These can be closed and convex structures
such as fullerenes of various basic symmetry, as well as plane structures of
various lateral sizes, shapes, edges, and doping levels [59, 60, 61]. As
known, the symmetry of the graphene sublattice in GQD can be controlled by the
lateral size, shape, and type of the GQD edge [57, 58, 59, 60, 61, 62, 63, 64,
65, 66, 67, 68, 69, 70, 71, 72]. The behavior of GQD is quantitatively
different for nanostructures with zigzag and armchair edge [73], [74]. Thus,
GQDs have richer electronic properties than graphene unbounded in space [75,
76, 77]. Moreover, they are of particular interest because their nonlinear
optical properties can be controlled [1].

The 2D graphene unbounded in space formed in the type of a narrow ribbon leads
to the confinement of carriers in quasi-one-dimensional graphene nanoribbons
(with different topologies depending on the shape of the ribbon) [26]. An
important advantage of GQD over graphene nanoribbons [69], [78, 79, 80] is the
complete confinement of quasiparticles in space. For a bounded quantum system,
due to the lack of translational symmetry, the spectrum only consists of a
discrete set of energy levels, not bands. One of the advantages of limited
quantum systems is the ability to change their energy spectra and
wavefunctions by adjusting their size [81], [82] or external parameters such
as gate voltage or magnetic field [83, 84, 85]. Note that the HHG spectra in
GQDs can be affected by the constraint conditions — the lateral size of the
dot and/or the coupling parameters in the bounded GQDs [67], [57], [58]. The
HHG efficiency is expected to increase with increasing confinement since the
latter will limit the propagation of the electron wave packet [86].

Figure 1:  The geometrical structure of ($N_{a}$; $N_{z}$)-sized RGQD with
$N=N_{a}N_{z}$ atoms. (a) shows RGQD with zigzag edge along the $x$ axis with
$N=84$ atoms for $N_{a}=21$; $N_{z}=4$, (b) – RGQD with armchair edge along
the $x$ axis for $N_{a}=12$; $N_{z}=7$. The distance between nearest
neighboring atoms is $b\simeq 1.42\ \mathring{A}$.

In the present paper, we consider the multiphoton HHG in the RGQD–plane quasi-
zero-dimensional GQD of the rectangular shape [87, 88, 89], with zigzag and
armchair edges on the elongated side of the different lateral sizes induced by
intense coherent radiation in the nonlinear regime. The many-body Coulomb
interaction is taken into account in the extended Hubbard approach [60, 61],
[90]. Note that the considered GQDs, where the engineering of optoelectronic
properties is allowed, currently are accessible in practice [59], [68]. A
closed system of differential equations for a single-particle density matrix
at the multiphoton interaction of RGQD with a strong laser field is solved
numerically.

The paper is organized as follows. In Sec. II, we consider multiphoton
excitation and HHG in RGQD of the different lateral sizes, laser intensity,
and polarization of the wave field. The comparison with the case of the
armchair edge on the elonged side of RGQD was done. Conclusions are given in
Sec. III. Finally, in Sec. IV the set of equations for the single-particle
density matrix has been formulated, taking into account the many-body Coulomb
interaction.

##  2 Basic model, numerical results and discussions

We will study in RGQD the HHG of the plane quasi-monochromatic EM wave of
linear polarization, propagating perpendicular to the $xy$ plane, with
homogeneous quasi-periodic electric field strenth:

$\mathbf{E}\left(t\right)=\widehat{\mathbf{e}}E_{0}f\left(t\right)\cos\omega
t,$ (1)

where $\widehat{\mathbf{e}}$ is the unit polarization vector in the $xy$
plane, $E_{0}$ and $\omega$ respectivelly are the amplitude and carrier
frequency, $f\left(t\right)=\sin^{2}\left(\pi t/\mathcal{T}\right)$ is the
slowly changing envelope, $\mathcal{T}=40\pi/\omega$ is the impulse duration.
Fig. 1 shows the geometric structure of the zero-dimensional RGQD in $xy$
plane. There are $N_{a}$ and $N_{z}$ carbon atoms along the $x$ axis and the
$y$ axis, respectively. All carbon atoms are packed into a hexagonal lattice.
The total number of carbon atoms RGQD is $N=N_{a}N_{z}$. Note that, for
nanoribbons with a limited armchair edge width, the electronic structure and
energy spectrum of quasiparticles for $N_{a}>>N_{z}$ and $N_{a}<<N_{z}$, when
the RGQD will be similar to a nanoribbon, critically depend on the width of
the RGQD along the $y$ axis and the width along the $x$ axis, respectively
[70]. In relation to [70], we have the insulator case of RGQD (Fig. 1).

Figure 2:  (Color online) The eigenenergies in RGQD. (a–d) for zigzag edge
along the $x$ axis respectively correspond to
$N=84(N_{a}=21),164(N_{a}=41),244(N_{a}=61),324(N_{a}=81)$, atoms; and
$N_{z}=4$ for all cases. The same for RGQD of armchair edge along the $x$ axis
for (a–d) $N=84(N_{a}=12),168(N_{a}=24),252(N_{a}=36),336(N_{a}=48)$ atoms,
respectively; and $N_{z}=7$ for all cases.

Thus, we will use the microscopic nonlinear quantum theory of the HHG process
in the strong EM wave in RGQD by the tight-binding model (TB) [59], [91] for
the Hamiltonian $\widehat{H}$, taking into account the Coulomb interaction in
the generalized Hubbard approximation. The total Hamiltonian by the empirical
TB model [91] is given in Sec. IV. Using numerical diagonalization, we will
find the eigenstates $\psi_{\mu}\left(i\right)$ and eigenenergies
$\varepsilon_{\mu}$ ($\mu=0,1..N-1$). The results of numerical diagonalization
are shown in Fig. 2 for the different lateral sizes and edges of RGQD. As
shown in Fig. 2, without the tunneling all energy levels have degenerated. It
can also be seen from Fig. 2 that the density of states increases with an
increase in the number of lattice atoms. As will be seen later, this directly
affects the yield of HHG.

The basic theory is given in Seq. Appendixes. The quantum dynamics of RGQD in
a strong EM wave with periodic excitation is determined by a closed system of
differential equations (11) for the density matrix, which must be solved with
the certain initial conditions. We construct mainly the density matrix
$\rho_{0ij}^{\left(\sigma\right)}$ through the filling of electronic states in
valence band in accordance with the Fermi–Dirac distribution at zero
temperature

$\rho_{0ij}^{\left(\sigma\right)}=\sum_{\mu=N/2}^{N-1}\psi_{\mu}^{\ast}\left(j\right)\psi_{\mu}\left(i\right),$
(2)

with an eigenstate $\psi_{\mu}\left(i\right)$ of the Hamiltonian
$\widehat{H}_{0}$ (9). For the density matrix we solve numerically by
integration the equations of motion Eq. 11 over time with the standard fourth-
order Runge-Kutta algorithm.

The emission spectrum of harmonics is defined by the Fourier transform
$\mathbf{a}\left(\Omega\right)$ of the dipole acceleration:

$\mathbf{a}\left(t\right)=d^{2}\mathbf{d}/dt^{2}.$ (3)

The dipole momentum is:

$\mathbf{d}\left(t\right)=\left\langle\sum_{i\sigma}\mathbf{r}_{i}c_{i\sigma}^{\dagger}c_{i\sigma}\right\rangle,$
(4)

normalized by the number $N$ of lattice atoms, and the factor
$a_{0}=\overline{\omega}^{2}{d}_{0},$ where $\overline{\omega}=1\
\mathrm{eV}/\hbar$ and ${d}_{0}=1\ \mathrm{\mathring{A}}$. The power emitted
at a given frequency is proportional to
$\left|\mathbf{a}\left(\Omega\right)\right|^{2}$. We introduce the angle
$\theta$ between the polarization vector
$\widehat{\mathbf{e}}=\left\\{\cos\theta,\sin\theta\right\\}$ and the $x$ axis
directed by the alonged edge (see Fig. 1). The relaxation rate is taken
$\hbar\gamma=50\ \mathrm{meV}$. The electron-electron interaction (EEI) energy
of on-site Coulomb repulsion is $U\simeq 3\ \mathrm{eV}$ (see in Sec. IV). The
inter-site Coulomb repulsion energy is taken $V\simeq 0.9\ \mathrm{eV}$. The
energy between the nearest-neighbor atoms is $t_{ij}=2.7\ \mathrm{eV}$. To
clarify the main aspects of multiphoton HHG in RGQD, we assume that the
excitation frequency is $\omega=0.1\ \mathrm{eV}/\hbar$, which is much smaller
than the typical scales $t_{ij}$, $U$, $V$.

Figure 3: (Color online) The HHG emission spectra in the strong-field regime
via dipole acceleration Fourier transformation
$N^{-1}|a_{x}\left(\Omega\right)|/a_{0}$ (in arbitrary units) versus the
harmonic number for RGQD with zigzag edge on the elongated side of different
number $N$ of atoms. The relaxation rate is $\hbar\gamma=50\ \mathrm{meV}$.
The linearly polarized wave frequency is $\omega=0.1\ \mathrm{eV}/\hbar$ and
wavefield strength is $E_{0}=0.1\ \mathrm{V/\mathring{A}}$ at fixed angle
$\theta=0$. The spectra are shown for moderate typical EEI energies: $U=3\
\mathrm{eV}$, $V=0.9\ \mathrm{eV}$.  Figure 4: (Color online) The same as for
Fig. 3 but for RGQD of armchair edge along the $x$ axis.

Next, we will consider the origin of the HHG in the quantum dot. There are two
contributions to the current: the electron/hole transitions within
unoccupied/occupied states and the electron-hole creation (transitions from
occupied states to unoccupied ones) and subsequent recombination. The first
transitions make a contribution only to low harmonics and is analogous to the
intraband current in a semiconductor, while the latter makes the main
contribution in the high-frequency part corresponding to the interband
current, which represents recombination/creation of electron-hole pairs. This
picture is analogous to HHG in solid-state systems. To separate these
contributions in the dipole acceleration spectrum we made a change of the
basis via formula
$\rho_{ij}=\sum_{\mu^{\prime}}\sum_{\mu}\psi_{\mu^{\prime}}^{\ast}\left(j\right)\rho_{\mu\mu^{\prime}}\psi_{\mu}\left(i\right)$,
where $\rho_{\mu\mu^{\prime}}$ is the density matrix in the energetic
representation. Hence intraband part of dipole acceleration defines as

$\mathbf{d}_{\mathrm{intra}}\left(t\right)=\sum_{\mu,\mu^{\prime}=N/2}^{N-1}\rho_{\mu\mu^{\prime}}\left(t\right)\mathbf{d}_{\mu^{\prime}\mu}+\sum_{\mu,\mu=0}^{N/2-1}\rho_{\mu\mu^{\prime}}\left(t\right)\mathbf{d}_{\mu^{\prime}\mu},$
(5)

and interband part will be

$\mathbf{d}_{\mathrm{inter}}\left(t\right)=2\sum_{\mu^{\prime}=N/2}^{N-1}\sum_{\mu=0}^{N/2-1}Re\rho_{\mu\mu^{\prime}}\left(t\right)\mathbf{d}_{\mu^{\prime}\mu},$
(6)

where dipole transition matrix elements are:

$\mathbf{d}_{\mu^{\prime}\mu}=\sum_{i}\psi_{\mu^{\prime}}^{\ast}\left(i\right)\mathbf{r}_{i}\psi_{\mu}\left(i\right).$
(7)

Figure 5: (Color online) The same as for Fig. 3 but for different angles
$\theta$ composed between in-plane electric field and $x$ axis for $N=164$
atoms of RGQD. Figure 6: (Color online) The same as for Fig. 5 but for
$a_{y}$–component.

In the following we will consider two contributions intraband (5) and
interband (6), to the total current numerically.

In Figs. 3, 4 the HHG spectrum determined by $|a_{x}\left(\Omega\right)|$ is
presented in RGQD of the different lateral sizes of correpondently zigzag and
armchair edge on the elonged side. The linearly polarized EM wavefield is
taken at polarization angle $\theta=0,$ $E_{0}=0.1$ V$\mathrm{/\mathring{A},}$
$\omega=0.1\ \mathrm{eV}/\hbar$. The EEI energies are $U\simeq 3\
\mathrm{eV}$, $V\simeq 0.9\ \mathrm{eV}$. As expected, the appearance of odd
harmonics in the HHG spectrum is associated with the conservation of inversion
symmetry in RGQD, as in ordinary graphene [1], [56]. This is a result of the
interference of the two different contributions, intraband (5) and interband
(6). In Figs. 3, 4 we see the typical nonperturbative behavior of HHG spectra
with a multiple plateau structure. In these cases, the dominant plateau
shifted towards higher frequencies with an increase in the number of atoms. In
particular, as shown in Fig. 3, at $N=164$, $N=244$, and $N=324$ we have an
effective generation of harmonics from the 10th to the 60th, the 20th to the
90th, and from the 20th to the 100th, respectively. When, as shown in Fig. 4,
in RGQD of armchair edge in RGQD, in particular for $N=168$ and $N=336$, we
have an efficient generation of harmonics from 10th to 40th and 70th to 100th,
respectively. This is related to the energy spectrums in Fig. 2, in which,
with an increase in the number of lattice atoms, new higher energy states with
an increasing density of states appear. As shown in Figs. 3 and 4, there are
main differences in the higher harmonics spectra for the close lateral size
RGQDs of different edges on the elonged sides. Thus, the HHG spectra maxima
and cutoff energies are different in RGQDs of zigzag and armchair edges on the
elonged sides. As will be shown by the following (see Fig. 14, red bar), for
RGQD at some angles $\theta>0$ the magnitudes of high-order harmonics in the
elonged zigzag case can become larger than at the armchair edge case,
especially for $a_{y}$.

Figure 7: (Color online) The same as for Fig. 5 but for $N=168$ atoms of RGQD
of armchair edge along the $x$ axis. Figure 8: (Color online) The same as for
Fig. 7 but for $a_{y}$–component.

Furthermore, RGQD also has a strong anisotropic nonlinear response, depending
on the orientation of the EM wavefield relative to the elonged side [57],
[58]. For RGQD at $\theta=0$ and $\theta=\pi/2$, the harmonic polarization
direction coincides with the incident wave polarization direction. At other
angles $0<\theta<\pi/2$, harmonics appear with polarization vectors
perpendicular to the pump wavefield. To reveal the dependence of the HHG
spectra on the orientation of the EM wavefield in Figs. 5 and 6 show the HHG
spectra for different $\theta$ of the pump wavefield relative to the $x$ axis
in the RGQD of zigzag edge for $N=164$, for $|a_{x}|$ and $|a_{y}|$
components, respectively. In contradistinction to the isotropic HHG spectrum
for unbounded in space graphene at the low-frequency pump wave, the spectrum
is anisotropic at a high-frequency wave, when, due to the symmetry of the
carbon hexagonal cell, the optical response with respect to the driven field
polarization is periodic with a period of $\pi/3$ [92]. For the RGQD we have
strong anisotropy. In particular, for HHG process in the case of more than the
first eight harmonics for the $y$ component the angles
$\pi/12\lesssim\theta\lesssim\pi/6$ are preferable, while for the
$x$-component we have a maximum for the angle $\pi/6$. Moreover, different
polarization angles lead to different maxima in the harmonic spectra and
cutoff energies. In Figs. 7, 8 the same is shown for RGQD of armchair edge.
For this case we see a completely different picture. In particular, for the
x-component the angles $0\lesssim\theta\lesssim\pi/12$ are preferable,
meanwhile for the $y$-component we have a maximum near the angle
$\theta\sim\pi/4$. The latter is also a consequent of the rich spectra of
eigenstates in Fig. 2 with different symmetry. Note that during the generation
process, only odd harmonics appear in the RGQD, regardless of its orientation.
This is due to the inversion symmetry of the sublattice of RGQD.

Figure 9: (Color online) The color bars represent the HHG emission rate in
strong field regime in the logarithmic scale via dipole acceleration Fourier
transformation $N^{-1}|a_{x}\left(\Omega\right)|/a_{0}$ (in arbitrary units)
versus harmonic order and the wavefield strength $E_{0}$, at fixed angle
$\theta=0$ for (a) $N=164$ at zigzag edge and (b) $N=168$ at armchair edge of
RGQD on the elonged side. The wave frequency is $\omega=0.1\
\mathrm{eV}/\hbar$ and the EEI energies are $U=3\ \mathrm{eV}$ and $V\simeq
0.9\ \mathrm{eV}$. The relaxation rate is $\hbar\gamma=50\ \mathrm{meV}$.
Figure 10: (Color online) The color bars represent the dipole momenta in the
logarithmic scale $N^{-1}|d_{x}|/d_{0}$ (in arbitrary units) versus eigenstate
index difference and state index. (a–c) respectively correspond to
$N=164,324,500$ atoms of RGQD with a zigzag edge along the $x$ axis. (d–f)
show the same but versus eigenstate energy difference and state index.

Next, we consider the HHG spectra depending on the intensity of the pump wave.
Fig. 9 (a) shows the HHG spectra via the field strength and the harmonic order
for a fixed frequency and EEI energies $U=3$ $\mathrm{eV}$, $V\simeq 0.9\
\mathrm{eV}$ in RGQD of zigzag elonged edge with $N=164$ atoms. We will
mention, that for HHG it is important to increase the HHG emission efficiency
and increase the harmonic cutoff threshold. As shown in Fig. 9 (a), within
each plateau, the cutoff harmonic increases linearly with increasing wavefield
strength. Then, reaching harmonics $n_{\mathrm{cut}}\simeq 110$. Note that the
linear dependence of the cutoff harmonics on the wave intensity is similar to
HHG through discrete levels [93, 94, 95], or in crystals with linear energy
dispersion [45], [96]. To compear in Fig. 9 (b) we show the same for RGQD of
the armchair edge on the elongated side. Fig. 9 has demonstrated the same
multi-plateau behavior. Many differences appear in Fig. 9 (a) and (b), in
particular, the absence of the rates of harmonic orders less $n=30$ in Fig. 9
(a), depending on the dipole momentum transition matrix elements (4) in the
RGQD of the zigzag or armchair edges are defined by the internal static
property of the system of the quantum dot.

To reveal such differences via edge on the elongated side of the RGQD for
different lateral sizes, in Figs. 10 and 11 was shown the $x$–component of the
dipole momentum $|\mathbf{d}_{\mu^{\prime}\mu}|$ (7) via the eigenstate index
difference and energy difference (transition energy), correspondently
demonstrated for a zigzag and armchair edges.

Figure 11: (Color online) The same as for Fig. 10 but for RGQD with an
armchair edge for (a–c), as for as (d–f), corresponding to $N=168,336,504$
atoms, respectively.

In Figs. 10, 11 appear, that with a change in the state index, new higher
energy states with an increasing density of states appear, and the states are
strongly different for zigzag and armchair edges of RGQD (connecting with
different topologies depending on the edge on the elongated side of the
quantum dot). As it is shown in Figs. 10, 11, with the growth in the number of
carbon atoms the corresponding HHG rates related by the transitions with the
more significant change of the state index occur relating to interband
transitions, moreover, it prevails in the case of armchair edge. Despite this,
as shown in Fig. 10, 11 for (c-f), in both cases of RGQD edges, the energy
difference of transitions tends to saturation with the increasing number of
carbon atoms.

Is in the particular interest the dipole momentum $y$–component. The
transitions, relating with $d_{y}$, appear mainly due to the strong coherent
EM radiation. In Figs. 12, 13 $d_{y}$–components are shown for different edges
and atom numbers of RGQD. As Figs. 12 and 13 show, that although the
$d_{y}$–component by an order of magnitude is smaller than the
$d_{x}$–component, the HHG rates are defined by correspondent transitions
between eigenstates with higher index difference, than in the case of
$d_{x}$–component; and for the RGQD of the armchair edge, the dipole momentum
magnitude is larger than in the case of zigzag edge on the elonged side. In
contradiction with Fig. 13 of the quantum dot of the elonged armchair edge, in
Fig. 12 for the zigzag edge, the middle isolated lines corresponding to the
substantial increase in the magnitude of $d_{y}$–component is clearly visible,
which due to interband transitions (defined by respective state index
difference). The latter has the main contribution to the high-frequency part
corresponding to the interband current, which is connected with
recombination/creation of electron-hole pairs, and it will correspondently
increase the HHG rates. Moreover, Figs. 12, 13 for (d–f) show the satiation of
the dipole momentum magnitude depending on the state energy difference with
the atom number growth, as in Figs. 10, 11 for (d–f).

Figure 12: (Color online) The same as for Fig. 10 but for $d_{y}$–component.
Figure 13: (Color online) The same as for Fig. 11 but for $d_{y}$–component.

In relation to the mentioned, Fig. 14 shows the high harmonic spectrum versus
the angle $\theta$ of the EM wavefield relative to the elonged side, and
harmonic order. (a) and (c) are demonstrated $|a_{x}|$ and $|a_{y}|$,
respectively, in RGQD of zigzag edge on the elonged side. (b) and (d)
correspond $|a_{x}|$ and $|a_{y}|$, respectively, in RGQD of armchair edge.
Fig. 14 is demonstrated for $n=120$ harmonic numbers in accordance with the
energistate index range $(-120,120)$ in Fig. 2 (c) for correspondent number of
atoms. As shown in Fig. 14, in contradistinction to the RGQD with an armchair
edge, in the case of a zigzag edge on the elonged side of RGQD, if for the
$y$–component we have maxima (red and black dots) in the high-frequency part
of HHG corresponding to transitions from occupied states to unoccupied ones
(the electron-hole creation and subsequent recombination), meanwhile for the
$x$–component the preferable angles (in accordance with Fig. 5–8) appear in
the low-frequency part of HHG spectra for the electron/hole transitions
between nearest-neighbor eigenstates (within unoccupied/occupied states).
Former transitions correspond to the interband current, which represents the
recombination/creation of electron-hole pairs, while the latter corresponds to
the intraband current. As shown in Fig. 14, that is defined by respective
state energy difference ${n}\hbar\omega$. The interband current sufficiently
enhances the HHG process rate in the elonged zigzag edge case.

Figure 14: (Color online) The color bars represent the HHG emission rate in
strong field regime in the logarithmic scale via dipole acceleration Fourier
transformation $N^{-1}|a_{i}\left(\Omega\right)|/a_{0}$ (in arbitrary units,
$i=x,y$ ) versus $\theta$ and harmonic order for (a) $N=244$ at zigzag edge
and (b) $N=252$ at armchair edge of RGQD. The wavefield strength is
$E_{0}=0.1\ \mathrm{V/\mathring{A}}$ with frequency $\omega=0.1\
\mathrm{eV}/\hbar$, and the EEI energies are $U=3\ \mathrm{eV}$ and $V\simeq
0.9\ \mathrm{eV}$. The relaxation rate is taken $\hbar\gamma=50\
\mathrm{meV}$.

##  3 Conclusion

We have investigated the generation of infrared laser field harmonics in the
extreme nonlinear strong field regime on the RGQD of different lateral sizes
and edges. We treated the system evolution in the strong laser field on the
basis of the Hartree-Fock approximation. The multiparticle Coulomb interaction
we consider in the extended Hubbard approximation. The numerical results
revealed the typical non-perturbative multiple-plateaus of the spectrum of
HHG, the boundaries of which are determined by the internal excitation lines
of the dot, and the differences related to the RGQD edge. Note, that the
dominant plateau shifted towards higher frequencies with an increase in the
number of atoms. Therefore, by changing the lateral size of the nanostructure,
it is possible to increase the orders of harmonics within the main plateau. In
addition, the cutoff photon energy also shifts toward blue with an increase in
the transverse size of the nanostructure. We have shown that the HHG spectra
have a strong anisotropy depending on the orientation of the electromagnetic
wave field strength relative to the RGQD zigzag edge on the elonged side.
Different polarization angles lead to different maxima in the harmonic spectra
and cutoff energies. The isolated line corresponds to mainly an increase in
the dipole momentum $y$–the component that appears in the HHG process on RGQD
along the zigzag edge, which causes a sufficient increase in the high-order
harmonics yield. The obtained results show that RGQD can serve as a medium for
HHG at interacting with a strong laser field, for the appearance of new energy
states, including that HHG probability increases with growth in the number of
dot atoms. This can be a proper way to enhance the quantum yield of the HHG
and the energy of photons in this process on the graphene-like quantum dots.

##  4 Appendixes: Interaction Hamiltonian

Here we will present the total Hamiltonian by the empirical TB model [91] in
the form:

$\widehat{H}=\widehat{H}_{\mathrm{0}}+\widehat{H}_{\mathrm{int}},$ (8)

where

$\widehat{H}_{0}=\frac{1}{2}\sum_{\left\langle
i,j\right\rangle}V_{ij}n_{i}n_{j}+\frac{U}{2}\sum_{i\sigma}n_{i\sigma}n_{i\overline{\sigma}}-\sum_{\left\langle
i,j\right\rangle\sigma}t_{ij}c_{i\sigma}^{\dagger}c_{j\sigma}.$ (9)

is the free RGQD Hamiltonian. Here $c_{i\sigma}^{\dagger}$ is the operator of
creation of an electron with spin polarization
$\sigma=({\uparrow,\downarrow})$ at site $i$, and $\left\langle
i,j\right\rangle$ is the summation over nearest neighbor sites with the
transfer energy $t_{ij}$ ($\overline{\sigma}$ is the opposite to spin
polarization); and ${n}_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$ is the
electron density operator with total electron density for the site $i$:
${n}_{i}={n}_{i\uparrow}+{n}_{i\downarrow}$. The first and the second terms in
free Hamiltonian (9) correspond to the EEI within the extended Hubbard
approximation ($\widehat{H}_{ee}$) with inter-site $\sim V_{ij}$ and on-site
$\sim U$ Coulomb repulsion energies. The inter-site Coulomb repulsion is
described by the distance $d_{ij}$ between the nearest-neighbour pairs varied
over the system: $\sim V_{ij}=Vd_{min}/d_{ij}$ ($d_{min}$ is the the minimal
nearest-neighbor distance). For the all calculations we have taken $V=0.3U$
[97], [98]. The third term in (9) is the kinetic energy part of the TB
Hamiltonian with tunneling matrix element $t_{ij}$ neighboring sites. The
hopping integral $t_{ij}$ between the nearest-neighbor atoms of GQDs can be
determined experimentally, and is usually taken to be $t_{ij}=2.7\
\mathrm{eV}$ [1]. Note, that we neglected the lattice vibrations in the
Hamiltonian.

The laser-RGQD interaction is described in the length-gauge via the pure
scalar potential:

$\widehat{H}_{\mathrm{int}}=e\sum_{i\sigma}\mathbf{r}_{i}\mathbf{E}\left(t\right)c_{i\sigma}^{\dagger}c_{i\sigma},$

where $\mathbf{r}_{i}$ is the position vector, $e$ is the elementary charge.
We will obtain evolutionary equations for the single-particle density matrix
$\rho_{ij}^{\left(\sigma\right)}=\left\langle
c_{j\sigma}^{\dagger}c_{i\sigma}\right\rangle$ from the Heisenberg equation
$i\hbar\partial\widehat{L}/\partial t=\left[\widehat{L},\widehat{H}\right]$.
Let us the system relaxes at a rate $\gamma$ to the equilibrium
$\rho_{0ij}^{\left(\sigma\right)}$ distribution. The EEI will be considered
under the Hartree-Fock approximation, and to describe a closed set of
equations for the single-particle density matrix
$\rho_{ij}^{\left(\sigma\right)}$, we will assumed the Hamiltonian (9) in the
form:

$\widehat{H}_{0}^{HF}\simeq-\sum_{\left\langle
i,j\right\rangle\sigma}t_{ij}c_{i\sigma}^{\dagger}c_{j\sigma}+U\sum_{i}(\overline{n}_{i\uparrow}-\overline{n}_{0i\uparrow})n_{i\downarrow}$
$+U\sum_{i\sigma}(\overline{n}_{i\downarrow}-\overline{n}_{0i\downarrow})n_{i\uparrow}+\sum_{\left\langle
i,j\right\rangle}V_{ij}(\overline{n}_{j}-\overline{n}_{0j})n_{i}$
$-\sum_{\left\langle
i,j\right\rangle\sigma}V_{ij}c_{i\sigma}^{\dagger}c_{j\sigma}\left({\left\langle{c_{i\sigma}^{\dagger}c_{j\sigma}}\right\rangle}-{\left\langle{c_{i\sigma}^{\dagger}c_{j\sigma}}\right\rangle}_{0}\right),$
(10)

with
$\overline{n}_{i\sigma}={\left\langle{c_{i\sigma}^{\dagger}c_{i\sigma}}\right\rangle}=\rho_{ii}^{\left(\sigma\right)}$.
Thus, the following equation for the density matrix is obtained:

$i\hbar\frac{\partial\rho_{ij}^{\left(\sigma\right)}}{\partial
t}=\sum_{k}\left(\tau_{kj\sigma}\rho_{ik}^{\left(\sigma\right)}-\tau_{ik\sigma}\rho_{kj}^{\left(\sigma\right)}\right)+\left(V_{i\sigma}-V_{j\sigma}\right)\rho_{ij}^{\left(\sigma\right)}$
$+e\mathbf{E}\left(t\right)\left(\mathbf{r}_{i}-\mathbf{r}_{j}\right)\rho_{ij}^{\left(\sigma\right)}-i\hbar\gamma\left(\rho_{ij}^{\left(\sigma\right)}-\rho_{0ij}^{\left(\sigma\right)}\right),$
(11)

and the matrixes $V_{i\sigma}$, $\tau_{ij\sigma}$ are approximated by the
density matrix $\partial\rho_{ij}^{\left(\sigma\right)}$:

$V_{i\sigma}=\sum_{j\alpha}V_{ij}\left(\rho_{jj}^{\left(\alpha\right)}-\rho_{0jj}^{\left(\alpha\right)}\right)+U\left(\rho_{ii}^{\left(\overline{\sigma}\right)}-\rho_{0ii}^{\left(\overline{\sigma}\right)}\right),$
(12)
$\tau_{ij\sigma}=t_{ij}+V_{ij}\left(\rho_{ji}^{\left({\sigma}\right)}-\rho_{0ji}^{\left({\sigma}\right)}\right).$
(13)

In this representation, the initial value of density matrix
${\left\langle{c_{i\sigma}^{\dagger}c_{j\sigma}}\right\rangle}_{0}$ is defined
via TB Hamiltonian ${\widehat{H}_{0}^{TB}}=-\sum_{\left\langle
i,j\right\rangle\sigma}t_{ij}c_{i\sigma}^{\dagger}c_{j\sigma}$. We will
numerically diagonalize the TB Hamiltonian $\widehat{H}_{0}$. That is, in the
static limit the Hartree-Fock Hamiltonian vanishes
$\widehat{H}_{ee}^{HF}\simeq 0$, and the EEI in Hartree-Fock limit is included
in empirical hopping integral between the nearest-neighbor atom $t_{ij}$ which
is chosen to be close to experimental data [59]. Thus, the EEI in the Hartree-
Fock approximation is correspond only for quantum dynamics initiated by the
pump laser field.

The authors are deeply grateful to prof. H. K. Avetissian and Dr. G. F.
Mkrtchian for permanent discussions and valuable recommendations. This work
was supported by the Science Committee of RA in Frames of Project 21AG-1C014.

## References

## References

  * [1] Novoselov K S, Geim A K, Morozov S V at el. 2004 Science 306 666
  * [2] Geim A K 2009 Science 324 1530
  * [3] Castro Neto A H, Guinea F, Peres N M R at el. 2009 Rev. Mod. Phys. 81 109
  * [4] Brabec T and Krausz F 2000 Rev. Mod. Phys. 72 545
  * [5] Corkum P B 1993  Phys. Rev. Lett. 71 1994
  * [6] Agostini P and DiMauro L F 2004  Rep. Prog. Phys. 67 813
  * [7] Kohler M C, Pfeifer T, Hatsagortsyan K Z, Keitel C H 2012 Adv. Atom. Mol. Opt. Phys. 61 159
  * [8] Avetissian H K 2016 Relativistic Nonlinear Electrodynamics: The QED Vacuum and Matter in Super-Strong Radiation Fields (New York: Springer)
  * [9] Ferray M, L’Huillier A, Li X F, Lompre L A at el. 1988 Journal of Physics B 21 L31
  * [10] Ghimire S, DiChiara A D, Sistrunk E at el. 2011 Nature Physics 7 138
  * [11] Schubert O, Hohenleutner M, Langer F at el. 2014 Nature Photonics 8 119
  * [12] Vampa G, Hammond T J, Thire N at el. 2015 Nature 522 462
  * [13] Ndabashimiye G, Ghimire S, Wu M at el. A 2016 Nature 534 520
  * [14] You Y S, Reis D A, and Ghimire S 2017 Nature Physics 13 345349
  * [15] Liu H, Guo C, Vampa G at el. 2018 Nature Physics 14 1006
  * [16] Avetissian H K, Avetisyan V N, Avchyan B R, Mkrtchian G F 2022 Phys. Rev. A 106 033107
  * [17] Mikhailov S A and Ziegler K 2008 J. Phys. Condens. Matter 20 384204
  * [18] Avetissian H K, Avetissian A K, Mkrtchian G F, Sedrakian Kh V 2012 Phys. Rev. B 85 115443
  * [19] Avetissian H K, Avetissian A K, Mkrtchian G F, Sedrakian Kh V 2012 J. Nanophoton. 6 061702
  * [20] Avetissian H K, Mkrtchian G F, Batrakov K G at el. 2013 Phys. Rev. B 88 165411
  * [21] Bowlan P, Martinez-Moreno E, Reimann K at el. 2014 Phys. Rev. B 89 041408
  * [22] Al-Naib I, Sipe J E, and Dignam M M 2015 New J.Phys. 17 113018
  * [23] Chizhova L A, Libisch F, and Burgdorfer J 2016 Phys. Rev. B 94 075412
  * [24] Avetissian H K and Mkrtchian G F 2016 Phys. Rev. B 94 045419
  * [25] Avetissian H K, Ghazaryan A G, Mkrtchian G F, Sedrakian Kh V 2017 J. Nanophoton.  11 016004
  * [26] Avetissian H K, Avchyan B R, Mkrtchian G M, Sargsyan K A 2020 J. Nanophoton.  14 026018
  * [27] Chizhova L A, Libisch F, and Burgdorfer J 2017 Phys. Rev. B 95 085436
  * [28] Dimitrovski D, Madsen L B, and Pedersen T G 2017 Phys. Rev. B 95 035405
  * [29] Yoshikawa N, Tamaya T, and Tanaka K 2017 Science 356 736
  * [30] Golub A, Egger R, Muller C, and Villalba-Chavez S 2020 Phys. Rev. Lett 124 110403
  * [31] Avetissian H K and Mkrtchian G F 2018 Phys. Rev. B 97 115454
  * [32] Avetissian A K, Ghazaryan A G, Sedrakian Kh V 2019  J. Nanophoton. 13(3) 036010
  * [33] Ghazaryan A G, Sedrakian Kh V 2019 J. Nanophoton. 13(4) 046004
  * [34] Ghazaryan A G, Sedrakian Kh V 2019 J. Nanophoton. 13(4) 046008
  * [35] Avetissian A K, Ghazaryan A G, Sedrakian Kh V, Avchyan B R 2017 J. Nanophoton. 11 036004
  * [36] Avetissian A K, Ghazaryan A G, Sedrakian Kh V, Avchyan B R 2018 J. Nanophoton. 12 016006
  * [37] Avetissian H K, Avetissian A K, Ghazaryan A G at el. 2020 J. Nanophoton. 14 026004
  * [38] Avetissian H K, Avetissian A K, Avchian B R, and Mkrtchian G F 2019 Phys. Rev. B 100 035434
  * [39] Bludov Yu, Peres N, M. Vasilevskiy M 2020 Phys. Rev. B 101 075415
  * [40] Liu H, Li Y, You Y S at el. 2017 Nature Physics 13 262
  * [41] Mkrtchian G F, Knorr A, and Selig M, 2020 Phys. Rev. B 100125401
  * [42] Avetissian H K, Mkrtchian G F, Hatsagortsyan K Z 2020 Phys. Rev. Research 2 023072
  * [43] Breton G Le, Rubio A, Tancogne-Dejean N 2018 Phys. Rev. B 98 165308
  * [44] Avetissian H K, Avetissian A K, Avchyan B R, Mkrtchian G F 2018 J. Phys. Condens. Matter 30 185302
  * [45] Avetissian H K, Avetissian A K, Avchyan B R, Mkrtchian G F 2019 Phys. Rev. B 100 035434
  * [46] Pedersen T G 2017 Phys. Rev. B 95 235419
  * [47] Avetissian H K and Mkrtchian G F 2019 Phys. Rev. B 99 085432
  * [48] Almalki S, Parks A M, Bart G at el. 2018 Phys. Rev. B 98 144307
  * [49] Cheng B, Kanda N, Ikeda T N at el. 2020 Phys. Rev. Lett 124 117402
  * [50] Cao T, Li Z, and Louie S G 2015 Phys. Rev. Lett. 114 236602
  * [51] Seixas L, Rodin A S, Carvalho A, and Castro Neto A H 2016 Phys. Rev. Lett. 116 206803
  * [52] Sevinzli H 2017 Nano Lett. 17 258
  * [53] Zhang X, Zhu T, Du H at el. 2022 Phys. Rev. Research 4, 033026
  * [54] Faist J, Capasso F, Sivco D L at el. 1994 Science 264(5158) 553
  * [55] Zhang G P, Bai Y H 2020 Phys. Rev. B 101 081412(R)
  * [56] Avetissian H K, Ghazaryan A G, and Mkrtchian G F 2021 Phys. Rev. B 104 125436
  * [57] Avchyan B R, Ghazaryan A G, Sargsyan K A, and Sedrakian Kh V 2022 JETP 134(2) 125
  * [58] Avchyan B R, Ghazaryan A G, Sargsyan K A, Sedrakian Kh V 2022 Pis’ma ZETF 116(7) 426
  * [59] Guclu A D, Potasz P, Korkusinski M, Hawrylak P 2014 Graphene quantum dot (Berlin: Springer)
  * [60] Avetissian H K, Mkrtchian G F, and Knorr A. 2022 Phys. Rev. B 105 195405
  * [61] Avetissian H K, Mkrtchian G F 2022 Phys. Rev. A 105 063504
  * [62] Fujita M, Wakabayashi K, Nakada K, and Kusakabe K 1996 J. Phys. Soc. Jpn. 65 1920
  * [63] Nakada K, Fujita K, Dresselhaus G, and Dresselhaus M S 1996 Phys. Rev. B 54 17954
  * [64] Son Y-W, Cohen M L, and Louie S G 2006 Phys. Rev. Lett. 97 216803
  * [65] Ezawa M 2008  Physica E 40 1421
  * [66] Ezawa M 2006 Phys. Rev. B 73 045432
  * [67] Chu W, Xie Y, Duan S at el. 2010 Phys. Rev. B 82 125301
  * [68] Yoon H, Park M, Kim J at el. 2021 Chem. Phys. Rev. 2 031303
  * [69] Han M Y, Ozyilmaz B, Zhang Y, and Kim Ph 2007 Phys. Rev. Lett. 98 206805
  * [70] Brey L and Fertig H A 2006 Phys. Rev. B 73 235411
  * [71] Yamijala Sh, Mukhopadhyay M, and Pati S 2015 J. Phys. Chem. C 119, 21 12079
  * [72] Luryi S, Xu J, Zaslavsky A 2013 Future trends in microelectronics: Frontiers and innovations (New York: Wiley)
  * [73] Guclu A D, Potasz P, and Hawrylak P 2010 Phys. Rev. B 82 155445
  * [74] Feng X, Qin Y, and Liu Y 2018 Opt. Express 26 7132
  * [75] Guclu A D, Potasz P, Voznyy O at el. 2009 Phys. Rev. Lett. 103 246805
  * [76] Voznyy O, Guclu A D, Potasz P, and Hawrylak P 2011 Phys. Rev. B 83 165417
  * [77] Wang W L, Meng S, and Kaxiras E 2008 Nano Lett. 8 241
  * [78] Yang L, Cohen M L, and Louie S G 2007 Nano Lett. 7 3112
  * [79] Prezzi D, Varsano D, Ruini A at el. 2008 Phys. Rev. B 77 041404
  * [80] Prezzi D, Varsano D, Ruini A at el. 2008, and Molinari E 2011 Phys. Rev. B 84 041401
  * [81] Murray C B, Kagan C R, and Bawendi M G 2000  Annu. Rev. Mater. Sci. 30 545
  * [82] Bera D, Qian L, Tseng T-K, and Holloway P H 2010 Materials 3 2260
  * [83] Kumar A, Laux S E, and Stern F 1990  Phys. Rev. B 42 5166
  * [84] Ashoori R C, Stormer H L, Weiner J S at el. 1993  Phys. Rev. Lett. 71 613
  * [85] Hogele A, Seidl S, Kroner M at el. 2004  Phys. Rev. Lett. 93 217401
  * [86] Lewenstein M, Balcou Ph, Ivanov M Yu at el. 1994 Phys. Rev. A 49 2117
  * [87] Chen R B, Chang C P, and Lin M F 2010 Physica E 42 2812
  * [88] Chang C P, Huang Y C, Lu C L at el. 2006 Carbon 44 508
  * [89] Qin Y, Feng X, Liu Y 2019 Appl. Sci. 9 325
  * [90] Acharya S, Pashov D, Rudenko A N at el. 2021 npj Computational Materials 7 208
  * [91] Wallace P R 1947 Phys. Rev. 71 622
  * [92] Zurron-Cifuentes O, Boyero-Garcia R, Hernandez-Garcia C at el. 2019 Optics Express 27 7776
  * [93] Avetissian H K, Avchyan B R, and Mkrtchian G F 2012 J. Phys. B 45 025402
  * [94] Avetissian H K, Markossian A G, and Mkrtchian G F 2011 Phys. Rev. A 84 013418
  * [95] Avetissian H K, Markossian A G, and Mkrtchian G F 2011 Phys. Lett. A 375 3699
  * [96] Vampa G, McDonald C R, Orlando G at el. 2014 Phys. Rev. Lett. 113 073901
  * [97] Martin R L and Ritchie J P 1993 Phys. Rev. B  48 4845
  * [98] Zhang G P 2000 Phys. Rev. B  61 4377