# Inhomogeneous Dirac-Bloch Equations for Graphene Interacting with Structured
Light

Yaraslau Tamashevich  Faculty of Engineering and Natural Sciences, Tampere
University, Tampere 33720, Finland  Marco Ornigotti  Faculty of Engineering
and Natural Sciences, Tampere University, Tampere 33720, Finland

###### Abstract

We generalise the usual framework of Dirac-Bloch equations, used to compute
the nonlinear optical response of 2D materials excited by spatially uniform
optical pulses, to the case of structured light pulses. We derive the general
form of Dirac-Bloch equations in the presence of a spatially inhomogeneous
field, for the case of a graphene-like material. Then, as as example of
application of our method, we consider explicitly the case of angular momentum
carrying optical pulses interacting with graphene.

##  I Introduction

Since the discovery of graphene Novoselov  _et al._ (2004) a lot of effort in
theoretical and experimental works has been done in order to describe
different properties of graphene. This material exhibit exotic electronic and
optical properties such as minimal conductivity Katsnelson (2006), universal
absorbance Nair  _et al._ , anomalous Quantum Hall Effect Novoselov  _et al._
(2005); zha , and high nonlinear optical response Wright  _et al._ (2009);
Hafez  _et al._ (2020), to name a few. Inspired by the new physics discovered
in graphene, an intensive research has been carried out in the last decade on
other 2D materials as well, such as transition metal dichalcogenides (TMDs)
Manzeli  _et al._ (2017); Autere  _et al._ (2018),hexagonal boron nitride
(hBN) Neto  _et al._ (2009); Caldwell  _et al._ (2019), and black phosphorous
Geim and Novoselov . The most interesting characteristics of 2D materials,
however, is that they possess high nonlinear response, which can be several
orders of magnitude higher than their bulk counterparts. This enhancement
allows to observe a lot of nonlinear effects such as, for example, second-
Zeng  _et al._ (2013) and third-harmonic generation Kumar  _et al._ (2013),
four-wave mixing Hendry  _et al._ (2010), and optical limiting Dong  _et al._
(2015). While graphene is intrinsically centrosymmetric and therefore lacks
second order response, recent works have demonstrated how it is possible to
induce second-harmonic generation in graphene by inducing stress or strains
Ornigotti  _et al._ (2021) or by exploiting the fact that an ultra-intense
pulse can open its own gap near Dirac points, thus breaking centrosymmetry
Carvalho  _et al._ .

On a seemingly parallel trail, structured light and, in particular, orbital
angular momentum (OAM) carrying beams have been widely studied in the past
decades since the pioneering work of Allen and Woerdman in 1992 Allen  _et
al._ (1992), and represent a versatile tool used in several different
experimental situations, and show wide range of application in different
areas, namely, optical tweezers Franke-Arnold  _et al._ , quantum information
Fickler  _et al._ , microscopy Klar and Hell , spectroscopy Noyan and Kikkawa
, classical Wang  _et al._ and quantum Cozzolino  _et al._ communication, to
name a few.

Despite the considerable amount of research that has been done in both 2D
materials and structured light, investigating the benefits and potential new
physics that would result by combining these two disciplines is still in its
infancy, and it is only recently, that different scenarios, such as the
interaction of OAM beams with ring semiconductors Wätzel and Berakdar (2016),
photoexcitation of graphene Farías  _et al._ (2013) and structured nonlinear
optics in 2D materials have been studied. A comprehensive toolbox for
analysing the effects of structured light on electron dynamics in 2D
materials, however, is still not available.

In this work, we then extend the usual Dirac-Bloch equation approach Ishikawa
to the case of a spatially inhomogeneous, impinging electromagnetic field and
develop a framework that will allow us to calculate nonlinear current and
nonlinear optical response of graphene-like materials when excited by
structured light.

This work is organized as follows: in Sect. II we briefly summarize the main
results of graphene, in terms of low-energy Hamiltonian, its eigenstates and
eigenvalues, which will constitute the building blocks of our approach. Then,
in Sect. III we discuss the explicit form of the interaction Hamiltonian, in
the minimal coupling approximation, for the case of a structured impinging
light pulse, and give the explicit expression for the electric dipole moment
of graphene in $k$-space, and express the interaction Hamiltonian in the
electron-hole basis. Section IV is then dedicated to constructing Dirac-Bloch
equations for the case of impinging structured light pulses, while Sect. V
deals with calculating the Dirac current for this general case. The nonlinear
response of graphene under structured light excitation is then discussed in
Sect. VI. Finally, conclusions and future perspectives of this work are
discussed in Sect. VII.

##  II Low-Energy Hamiltonian for Graphene

To derive the low energy Hamiltonian for graphene in the presence of an
external field, we can start from the minimally coupled crystalline
Hamiltonian

$\hat{H}=\int
d^{3}r\hat{\psi}^{\dagger}(\textbf{r})\left[\frac{\left(\hat{\textbf{p}}-e\textbf{A}\right)^{2}}{2m}+\phi_{c}(\textbf{r})+V_{B}(\textbf{r})\right]\hat{\psi}(\textbf{r}),$
(1)

where $\phi_{c}(\textbf{r})$ is the Coulomb potential, $V_{B}(\textbf{r})$ is
the crystalline (periodic) potential, $\hat{\psi}(\textbf{r})$ is the electron
wavefunction operator, and $\textbf{A}(\textbf{r},t)$ is the vector potential.
From here henceforth, we are assuming to work in the Coulomb gauge, where
$\nabla\cdot\textbf{A}=0$, and we also neglect the Coulomb interaction terms,
since they don’t contribute significantly for the case of graphene Katsnelson
. Since for our purposes it is more convenient to work with the electric
field, rather than the vector potential, we can apply the Power-Zienau-Wooley
transformation Barnett and Radmore  to the above Hamiltonian, to get the usual
light-matter interction in the electric dipole approximation as follows:

$\hat{H}=\hat{H}_{0}+\hat{H}_{I},$ (2)

where

$\displaystyle\hat{H}_{0}$
$\displaystyle=\int\,d^{3}r\,\hat{\psi}^{\dagger}(\textbf{r})\left[\frac{\hat{\textbf{p}}^{2}}{2m}+V_{B}(\textbf{r})\right]\hat{\psi}(\textbf{r})$
(3a)
$\displaystyle\equiv\int\,d^{3}r\,\hat{\psi}^{\dagger}(\textbf{r})\mathcal{H}_{0}\hat{\psi}(\textbf{r}),$
$\displaystyle\hat{H}_{I}$
$\displaystyle=-e\int\,d^{3}r\,\hat{\psi}^{\dagger}(\textbf{r})\,\left[\textbf{r}\cdot\textbf{E}_{T}(\textbf{r},t)\right]\,\hat{\psi}(\textbf{r}),$
(3b)

where $\textbf{E}_{T}(\textbf{r},t)=-\partial_{t}\textbf{A}(\textbf{r},t)$ is
the (transverse) electric field. To calculate the graphene Hamiltonian, we use
the following Bloch representation for the electron wavefunction operator

$\displaystyle\hat{\psi}(\textbf{r})=\frac{1}{\sqrt{N}}\sum_{k}\sum_{n=0}^{N}e^{i\textbf{k}\cdot\textbf{R}_{n}}\Big{[}\phi_{A}(\textbf{r}-\textbf{R}_{n})\hat{a}_{k}$
(4) $\displaystyle+\phi_{B}(\textbf{r}-\textbf{R}_{n})\hat{b}_{k}\Big{]},$

where $\phi_{A,B}(\textbf{r})$ are $2p_{z}$-orbital wavefunctions at site $A$
and $B$ (and they form an orthogonal normalised basis set),
$\textbf{R}_{n}=j\textbf{a}_{1}+m\textbf{a}_{2}$ is a Bravais lattice vector
(see Fig. 1), and $\hat{a}_{k},\hat{b}_{k}$ are the electron site operators,
annihilating an electron on site A or B, respectively.

Figure 1: (a) Lattice structure of graphene. The unit cell, containing two
atoms (one belonging to lattice A, the other to lattice B), is shaded in gray.
The basis vectors used in this work are chosen to be
$\textbf{a}_{1}=(a/2)(3,\sqrt{3})$, and $\textbf{a}_{2}=(a/2)(3,-\sqrt{3})$.
(b) Reciprocal space lattice structure of graphene (the Wigner-Seitz cell is
shaded in grey). The reciprocal lattice basis vectors are defined as
$\textbf{b}_{1}=(2\pi/3a)(1,\sqrt{3})$ and
$\textbf{b}_{2}=(2\pi/3a)(1,-\sqrt{3})$.The two inequivalent Dirac points K
and K’ are defined as $\textbf{K}=(2\pi/3a)(1,-1/\sqrt{3})$, and
$\textbf{K'}=(2\pi/3a)(1,1/\sqrt{3})$,respectively. Figure 2: Band structure
of graphene. You can see 6 touching point of valence (orange) and conduction
(blue) band which are called Dirac points. Locally Dirac point can be
represented as a touching of two identical cones (see figure inside
rectangular) that allows there treating as a massless particles.

By substituting Eq. (4) into the first of Eqs. (3), it is not difficult to
see, that the free Hamiltonian $\hat{H}_{0}$ can be then reduced to

$\hat{H}_{0}=\gamma\sum_{k}\left[f(\textbf{k})\hat{a}^{\dagger}_{k}\hat{b}_{k}+\text{h.c.}\right],$
(5)

where $\gamma$ is the hopping amplitude and
$f(\textbf{k})=1+\exp{(i\textbf{k}\cdot\textbf{a}_{1})}+\exp{(i\textbf{k}\cdot\textbf{a}_{2})}$,
and h.c. stands for Hermitian conjugate. Diagonalising the above free
Hamiltonian leads to the band structure (eigenvalues) of graphene, i.e.,
$E(\textbf{k})=\pm\gamma|f(\textbf{k})|$, and its corresponding eigenvectors

$u_{\pm}(\textbf{k})=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}e^{i\frac{\varphi(\textbf{k})}{2}}\\\
\pm\,e^{-i\frac{\varphi(\textbf{k}}{2})}\end{array}\right),$ (6)

where $\varphi(\textbf{k})=\text{arg}\left[f(\textbf{k})\right]$. A plot of
the full band structure of graphene is shown in Fig. 2 We can transform the
site operators as well into the diagonal basis of the Hamiltonian
$\hat{H}_{0}$ ,obtaining

$\left(\begin{array}[]{c}\hat{e}_{k}\\\
\hat{h}^{\dagger}_{-k}\end{array}\right)=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}e^{-i\frac{\varphi(\textbf{k})}{2}}&e^{i\frac{\varphi(\textbf{k})}{2}}\\\
e^{-i\frac{\varphi(\textbf{k})}{2}}&-e^{i\frac{\varphi(\textbf{k})}{2}}\end{array}\right)\left(\begin{array}[]{c}\hat{a}_{k}\\\
\hat{b}_{k}\end{array}\right).$ (7)

Notice that this is equivalent to the electron-hole representation discussed,
for example, in Ref. Haug and Koch (2009), where $\hat{e}_{k}$ annihilates an
electron with momentum $k$, and $\hat{h}^{\dagger}_{-k}$ creates a hole with
momentum $-k$.

The free Hamiltonian in the electron-hole representation can be then written
as follows

$\hat{H}_{0}=\sum_{k}\,E(\textbf{k})\left(\hat{e}^{\dagger}_{k}\hat{e}_{k}+\hat{h}^{\dagger}_{-k}\hat{h}_{-k}\right).$
(8)

###  II.1 Low-Energy Approximation

We can expand the eigenvalues $E(\textbf{k})$ around the zeros of
$f(\textbf{k})$, which happen to be the two inequivalent Dirac points K and K’
[see Fig. 1 (b)] Katsnelson (2007) to obtain the low-energy approximation of
$\hat{H}_{0}$ around the Dirac points. To do that, we first redefine the
coordinates of the Dirac points as $\textbf{K}=(2\pi/3a)(1,\xi/\sqrt{3})$,
where $\xi=\pm 1$, and then expand $f(\textbf{k})$ around
$\textbf{q}=\textbf{k}+\textbf{K}$ obtaining

$\displaystyle f(\textbf{q})$ $\displaystyle\simeq$ $\displaystyle
f(\textbf{K})+\nabla_{k}f(\textbf{k})\Big{|}_{\textbf{K}}\textbf{q}$
$\displaystyle=$ $\displaystyle\frac{3ai}{2}\left(k_{x}-i\xi\,k_{y}\right)$

Substituting this result in the expression of the Hamiltonian $\hat{H}_{0}$
and defining the Fermi velocity as $\hbar v_{f}3a\gamma/2$, we get the
familiar form of the low-energy graphene Hamiltonian, i.e.,

$\hat{H}(\textbf{k})=\hbar v_{f}(\sigma_{x}k_{x}+\xi\sigma_{y}k_{y}),$ (10)

where $\xi=\pm 1$ is the valley index, and $\sigma_{i}$ are Pauli matrices.
The eigenvalues of $\hat{H}(\textbf{k})$ are now simply $E(\textbf{k})=\hbar
v_{f}|\textbf{k}|$. Notice, moreover, that in this approximation

$\displaystyle\varphi(\textbf{k})$ $\displaystyle=$
$\displaystyle\text{arg}\left[f(\textbf{k})\right]\simeq-\frac{\pi}{2}+\xi\arctan\left(\frac{k_{y}}{k_{x}}\right)$
$\displaystyle\equiv$ $\displaystyle-\frac{\pi}{2}+\varphi_{\xi}(\textbf{k}).$

And, that the eigenvectors $u_{\pm}(\textbf{k})$ acquire a valley dependence
through $\varphi_{\xi}(\textbf{k})$, i.e., $u_{\pm}(\textbf{k})\rightarrow
u_{\pm}^{\xi}(\textbf{k})$.

##  III Interaction Hamiltonian

We now substitute the Bloch Ansatz given by Eq. (4) into the expression of the
interaction Hamiltonian, i.e., the second of Eqs. (3). To calculate the
integrals, we proceed as before, first by introducing the change of variables
$\textbf{x}=\textbf{r}-\textbf{R}_{n}$, assuming only nearest neighbour
coupling, such that $\textbf{R}_{\ell}=\textbf{R}_{m}-\textbf{R}_{n}$ only
allows for $\ell=\\{0,\pm 1\\}$, and we make the extra assumption that the
electric field is constant over a single unit cell, i.e.,

$\textbf{E}_{T}(\textbf{x}+\textbf{R}_{n},t)\simeq\,\textbf{E}_{T}(\textbf{R}_{n},t).$
(12)

Notice that this approximation is consistent with previous work dealing with
semiconductor Bloch equations in the presence of inhomogeneous fields
Kuklinski and Mukamel (1991)

Since we are doing these calculation in the site basis, we will have to deal
with both diagonal terms, proportional to $\hat{a}^{\dagger}_{k}\hat{a}_{k}$
and $\hat{b}^{\dagger}_{k}\hat{b}_{k}$, and nondiagonal terms, proportional to
$\hat{a}^{\dagger}_{k}\hat{b}_{k}$. The former are related to the presence of
permanent electric dipoles in the system, and without loss of generality we
can neglect their contribution, since graphene (and most 2D materials, in
general) do not possess permanent dipoles. Therefore, we focus our attention
on the nondiagonal interaction terms, which describe light-induced transitions
and are regulated by the usual electric dipole moment. The general form of the
(nondiagonal) interaction Hamiltonian is then given by

$\displaystyle\hat{H}_{I}(\textbf{k})$ $\displaystyle=$
$\displaystyle\sum_{n}e^{i(\textbf{k}-\textbf{k'})\cdot\textbf{R}_{n}}\,\Bigg{[}\textbf{E}_{T}(\textbf{R}_{n},t)\cdot\textbf{d}_{CV}(\textbf{k'})$
(13) $\displaystyle-$ $\displaystyle e$
$\displaystyle\,f(\textbf{k'})\textbf{R}_{n}\cdot\textbf{E}_{T}(\textbf{R}_{n},t)\,\int\,\frac{d^{3}x}{N}\phi_{A}^{*}(\textbf{x})\phi_{B}(\textbf{x})\Bigg{]}\hat{a}^{\dagger}_{k}\hat{b}_{k^{\prime}}$

where the integral in the second line is zero because $\phi_{A}(x)$ and
$\phi_{B}(x)$ are orthogonal,

$\textbf{d}_{CV}(\textbf{k})=-e\sum_{\ell}e^{i\textbf{k}\cdot\textbf{R}_{\ell}}\int\,d^{3}x\phi_{A}^{*}(\textbf{x})\,\textbf{x}\,\phi_{B}(\textbf{x}-\textbf{R}_{\ell})\\\
$

is the (conduction-to-valence-band) transition dipole moment in reciprocal
space. Notice, that while normally, with spatially homogeneous fields, the
first term in the equation above would be readily calculated because
$\textbf{E}_{T}(\textbf{R}_{n},t)=\textbf{E}_{T}(t)$ constant in space, this
is not the case here, as we are considering a general, spatially inhomogeneous
field. In this case, we can introduce the angular spectrum representation for
the impinging field (i.e., its Fourier transform) as

$\textbf{E}_{T}(\textbf{R}_{n},t)=\sum_{q}\mathcal{E}_{T}(\textbf{q},t)e^{-i\textbf{q}\cdot\textbf{R}_{n}},$
(15)

and then introduce the generalised Rabi frequency as follows

$\Omega(\textbf{k},\textbf{q},t)=\frac{1}{\hbar}\mathcal{E}_{T}(\textbf{q},t)\cdot\textbf{d}_{CV}(\textbf{k}-\textbf{q}),$
(16)

so that the interaction Hamiltonian above can be written as

$\hat{H}_{I}(\textbf{k})=\hbar\sum_{q}\,\Omega(\textbf{k},\textbf{q},t)\hat{a}^{\dagger}_{k}\hat{b}_{k-q}.$
(17)

Notice, that in order for $\Omega(\textbf{k},\textbf{q},t)$ to be a bona-fide
Rabi frequency, it must be a dimensionless quantity, which implicitly means
that the dimension of the electric field must be $Vs/m$, which is consistent
with the fact that $|\textbf{E}|=|\partial_{t}\textbf{A}|\simeq\tau A_{0}$,
with $\tau$ being the characteristic timescale of the impinging electric field
(for an impinging pulse, $\tau$ is the duration of the pulse).

We are now in the position to write down the total interaction Hamiltonian
$\hat{H}_{I}$ in the site-basis as follows

$\hat{H}_{I}=\hbar\sum_{k,q}\left[\Omega(\textbf{k},\textbf{q},t)\hat{a}^{\dagger}_{k}\hat{b}_{k-q}+\Omega^{*}(\textbf{k},\textbf{q},t)\hat{b}^{\dagger}_{k}\hat{a}_{k-q}\right].$
(18)

This is the first result of our work. Notice, how accounting for the spatial
variation of the impinging electromagnetic pulse means that we are forced to
introduce a nonlocal coupling between the various states in $k$ space, through
the operators $\hat{a}^{\dagger}_{k}\hat{b}_{k-q}$, with the (transverse)
momentum of the light pulse acting as the connection between the various $k$
states available in the first Brillouin zone.

The interaction Hamiltonian in Eq. (18) can be put in a simpler and more
straightforward form, by representing the (complex) Rabi frequency in polar
coordinates, such that we get

$\displaystyle\hat{H}_{I}=\hbar\sum_{k,q}\left|\Omega(\textbf{k},\textbf{q},t)\right|\Big{[}e^{i\omega(\textbf{k},\textbf{q},t)}\hat{a}^{\dagger}_{k}\hat{b}_{k-q}$
(19)
$\displaystyle+e^{-i\omega(\textbf{k},\textbf{q},t)}\hat{b}^{\dagger}_{k}\hat{a}_{k-q}\Big{]},$

where
$\left|\Omega(\textbf{k},\textbf{q},t)\right|=\sqrt{\Omega_{x}^{2}(\textbf{k},\textbf{q},t)+\Omega_{y}^{2}(\textbf{k},\textbf{q},t)}$
is the modulus of the Rabi frequency, and

$\omega(\textbf{k},\textbf{q},t)=\text{arg}\left[\Omega(\textbf{k},\textbf{q},t)\right]=\arctan\left[\frac{\Omega_{y}(\textbf{k},\textbf{q},t)}{\Omega_{x}(\textbf{k},\textbf{q},t)}\right],$
(20)

is its phase.

###  III.1 Conduction-to-Valence band Dipole Moment in Reciprocal Space

The explicit expression of the conduction-to-valence band dipole moment can be
calculated analytically, in the electron-hole representation, using the
eigenstates (in reciprocal space) of the free Hamiltonian
$\hat{H}(\textbf{k})$, as detailed in Ref. Kelardeh  _et al._ (2015) as

$\textbf{d}_{CV}(\textbf{k})=\langle\psi_{c}(\textbf{k})|e\textbf{r}|\psi_{v}(\textbf{k})\rangle,$
(21)

where the conduction and valence band eigenstates are given as
$\ket{\psi_{c,v}(\textbf{k})}=e^{i\textbf{k}\cdot\textbf{r}}\ket{u_{\pm}(\textbf{k})},$
where the $+$ ($-$) sign refers to conduction (valence) band. Substitution
gives

$\textbf{d}_{CV}(\textbf{k})=\frac{e}{2}\nabla_{k}\varphi(\textbf{k}).$ (22)

Notice that in the vicinity of a Dirac point
$\varphi(\textbf{k})\simeq\frac{\pi}{2}+\xi\arctan(k_{y}/k_{x})$, and the
dipole moment becomes

$\displaystyle\textbf{d}_{CV}(\textbf{k})$ $\displaystyle=$
$\displaystyle\frac{e\xi}{2}\left(-\frac{k_{y}}{k_{x}^{2}+k_{y}^{2}}\hat{\textbf{x}}+\frac{k_{x}}{k_{x}^{2}+k_{y}^{2}}\hat{\textbf{y}}\right)$
$\displaystyle=$
$\displaystyle-\frac{e\xi}{2|\textbf{k}|}\left[\sin\varphi_{\xi}(\textbf{k})\hat{\textbf{x}}-\cos\varphi_{\xi}(\textbf{k})\hat{\textbf{y}}\right],$

or, by introducing the helicity basis
$\hat{\textbf{h}}_{\xi}=\frac{1}{\sqrt{2}}\left(\hat{\textbf{x}}+i\xi\hat{\textbf{y}}\right)$
Mandel and Wolf , we get

$\textbf{d}_{CV}(\textbf{K})=-\frac{ie\xi}{2\sqrt{2}|\textbf{k}|}\left(e^{-i\varphi_{\xi}(\textbf{k})}\hat{\textbf{h}}_{+}-e^{i\varphi_{\xi}(\textbf{k})}\,\hat{\textbf{h}}_{-}\right)$
(24)

###  III.2 Interaction Hamiltonian in the Electron-Hole Representation

We can transform the interaction Hamiltonian in the electron-hole
representation using the inverse of the change of basis given by Eq. (7),
i.e.,

$\left(\begin{array}[]{c}\hat{a}_{k}\\\
\hat{b}_{k}\end{array}\right)=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}e^{i\frac{\varphi(\textbf{k})}{2}}&e^{i\frac{\varphi(\textbf{k})}{2}}\\\
e^{-i\frac{\varphi(\textbf{k})}{2}}&-e^{-i\frac{\varphi(\textbf{k})}{2}}\end{array}\right)\left(\begin{array}[]{c}\hat{e}_{k}\\\
\hat{h}^{\dagger}_{-k}\end{array}\right).$ (25)

Substituting this into Eq. (18), and using the anti-commutation relations
$\\{\hat{h}^{\dagger}_{k},\hat{h}_{q}\\}=\delta_{k,q}$ we get

$\displaystyle\hat{H}_{I}=\hbar\sum_{k,q}\left|\Omega(\textbf{k},\textbf{q},t)\right|$
$\displaystyle\Big{\\{}\cos\Phi(\textbf{k},\textbf{q},t)\Big{(}\hat{e}^{\dagger}_{k}\hat{e}_{k-q}+\hat{h}^{\dagger}_{q-k}\hat{h}_{-k}$
$\displaystyle-\delta_{q,0}\Big{)}+i\sin\Phi(\textbf{k},\textbf{q},t)$
$\displaystyle\left(\hat{e}^{\dagger}_{k}\hat{h}^{\dagger}_{q-k}-\hat{h}_{-k}\hat{e}_{k-q}\right)\Big{\\}},$
(26)

where

$\Phi(\textbf{k},\textbf{q},t)=\frac{\Delta\varphi(\textbf{k},\textbf{q})}{2}-\omega(\textbf{k},\textbf{q},t),$
(27)

with
$\Delta\varphi(\textbf{k},\textbf{q})=\varphi(\textbf{k})+\varphi(\textbf{k}-\textbf{q})]$,
and $\omega(\textbf{k},\textbf{q},t)$ defined in Eq.(20). The form of the
interaction Hamiltonian in the electron-hole representation is quite
complicated, because it contains not only interband transitions
($e^{\dagger}e$ and $h^{\dagger}h$ terms) but also interband ones
($e^{\dagger}h^{\dagger}$ and $he$ terms). Moreover, written in this form it
is equivalent to the linear minimal coupling Hamiltonian in Ref. Stroucken
_et al._ (2011).

##  IV Dirac-Bloch Equations

To calculate the Dirac-Bloch equations, we first introduce the following two-
body populations and polarisation

$\displaystyle n_{e}(\textbf{k},\textbf{q})$
$\displaystyle=\langle\hat{e}^{\dagger}_{k}\hat{e}_{q}\rangle,$ (28a)
$\displaystyle n_{h}(\textbf{k},\textbf{q})$
$\displaystyle=\langle\hat{h}^{\dagger}_{-k}\hat{h}_{-q}\rangle,$ (28b)
$\displaystyle p(\textbf{k},\textbf{q})$
$\displaystyle=\langle\hat{h}_{-k}\hat{e}_{q}\rangle,$ (28c)

and calculate their time derivative with the help of the Heisenberg equation.
Using the anti-commutation rules
$\\{\hat{e}^{\dagger}_{k},\hat{e}_{k^{\prime}}\\}=\delta_{k,k^{\prime}}$ and
$\\{\hat{h}^{\dagger}_{-k},\hat{h}_{-k^{\prime}}\\}=\delta_{k,k^{\prime}}$ we
get, after some long but straightforward algebra, the following two-body
Dirac-Bloch equations

$\displaystyle-i\frac{d}{dt}n_{e}(\textbf{k},\textbf{k}^{\prime})$
$\displaystyle=\frac{E(\textbf{k})-E(\textbf{k}^{\prime})}{\hbar}n_{e}(\textbf{k},\textbf{k}^{\prime})+\sum_{q}\Bigg{\\{}G_{c}(k+q,q,t)n_{e}(\textbf{k}+\textbf{q},\textbf{k}^{\prime})-G_{c}(\textbf{k}^{\prime},\textbf{q},t)n_{e}(\textbf{k},\textbf{k}^{\prime}-\textbf{q})$
$\displaystyle-i\left[G_{s}(\textbf{k}^{\prime},\textbf{q},t)p^{*}(\textbf{k}^{\prime}-\textbf{q},\textbf{k})-G_{s}(\textbf{k},\textbf{q},t)p(\textbf{k}^{\prime},\textbf{k}-\textbf{q})\right]\Bigg{\\}},$
(29a) $\displaystyle-i\frac{d}{dt}n_{h}(\textbf{k},\textbf{k}^{\prime})$
$\displaystyle=\frac{E(\textbf{k})-E(\textbf{k}^{\prime})}{\hbar}n_{h}(\textbf{k},\textbf{k}^{\prime})+\sum_{q}\Bigg{\\{}G_{c}(\textbf{k},\textbf{q},t)n_{h}(\textbf{k}-\textbf{q},\textbf{k}^{\prime})-G_{c}(\textbf{k}^{\prime}+\textbf{q},\textbf{q},t)n_{h}(\textbf{k},\textbf{k}^{\prime}+\textbf{q})$
$\displaystyle+\left[G_{s}^{*}(\textbf{k}^{\prime}+\textbf{q},\textbf{q},t)p^{*}(\textbf{k},\textbf{k}^{\prime}+\textbf{q})+G_{s}(\textbf{k},\textbf{q},t)p(\textbf{k}-\textbf{q},\textbf{k}^{\prime})\right]\Bigg{\\}},$
(29b) $\displaystyle-i\frac{d}{dt}p(\textbf{k},\textbf{k}^{\prime})$
$\displaystyle=-\frac{E(\textbf{k}^{\prime})+E(\textbf{k})}{\hbar}p(\textbf{k},\textbf{k}^{\prime})+\sum_{q}\Bigg{\\{}-G_{c}(\textbf{k}^{\prime},\textbf{q},t)p(\textbf{k},\textbf{k}^{\prime}-\textbf{q})-G_{c}(\textbf{k}+\textbf{q},\textbf{q},t)p(\textbf{k}+\textbf{q},\textbf{k}^{\prime})$
$\displaystyle-
iG_{s}(\textbf{k}+\textbf{q},\textbf{q},t)n_{e}(\textbf{k}+\textbf{q},\textbf{k}^{\prime})+iG_{s}(\textbf{k}^{\prime},\textbf{q},t)\left[\delta_{k,k^{\prime}-q}-n_{h}(\textbf{k}^{\prime}-\textbf{q},\textbf{k})\right]\Bigg{\\}}$
(29c)

where $E(\textbf{k})$ [$E(\textbf{k}^{\prime})$] are the energy eigenvalues of
valence and conduction bands in state k ($\textbf{k}^{\prime}$), and we have
defined

$\displaystyle G_{c}(\textbf{k},\textbf{q},t)$
$\displaystyle=\left|\Omega(\textbf{k},\textbf{q},t)\right|\cos\Phi(\textbf{k},\textbf{q},t),$
(30a) $\displaystyle G_{s}(\textbf{k},\textbf{q},t)$
$\displaystyle=\left|\Omega(\textbf{k},\textbf{q},t)\right|\sin\Phi(\textbf{k},\textbf{q},t).$
(30b)

The equations above constitute a generalisation of the usual Dirac-Bloch
equations, to account for the presence of a spatially varying electromagnetic
pulse. This is the second result of our work. Notice, moreover, that written
in the above form, the generalised Dirac-Bloch equations have a similar form
that the semiconductor Bloch equations for an inhomogeneous electric field, as
given in Ref. 33, with the difference, that we don’t consider the effect of
the Coulomb potential in the specific case of graphene. Accounting for the
effect of Coulomb interactions on Dirac-Bloch equations, is of importance for
certain categories of 2D materials, such as TMDs, for example, and will be the
subject of a future publication.

##  V Current

As a final step towards building the nonlinear optical response of graphene,
we need to calculate the induced (nonlinear) current. To do so, we use the
usual definition of current, in terms of population and polarisation

$\hat{\textbf{j}}(t)=\sum_{k^{\prime},k,q}\hat{\psi}^{\dagger}_{k^{\prime}}U^{\dagger}(\textbf{k}^{\prime})\,\nabla_{k}\,h_{I}(\textbf{k},\textbf{q})\,U(\textbf{k}-\textbf{q})\hat{\psi}_{k-q},$
(31)

where, again, we work in the interaction picture, so that the $k$-dependent
interaction Hamiltonian above, in the band representation, is derived directly
from Eq. (19), i.e.,

$\displaystyle h_{I}(\textbf{k},\textbf{q})$ $\displaystyle=$
$\displaystyle\hbar\,\bm{\sigma}\cdot\bm{\Omega},$ (32)

where
$\bm{\Omega}=\Omega_{x}(\textbf{k},\textbf{q},t)\hat{\textbf{x}}+\Omega_{y}(\textbf{k},\textbf{q},t)\hat{\textbf{y}}$,
$\bm{\sigma}=\sigma_{x}\hat{\textbf{x}}+\sigma_{y}\hat{\textbf{y}}$,
$\hat{\psi}_{k}=(\hat{e}_{k}\hskip 2.84526pt\,\hat{h}^{\dagger}_{-k})^{T}$,
and $U(\textbf{k})$ being the matrix whose columns are the eigenstates of
$\hat{H}_{0}$ calculated before. If we then define the following $2\times 2$
matrices

$\displaystyle M_{x}$ $\displaystyle=$ $\displaystyle
U^{\dagger}\nabla_{k}h_{I}\,U\Big{|}_{x},$ (33) $\displaystyle M_{y}$
$\displaystyle=$ $\displaystyle U^{\dagger}\nabla_{k}h_{I}\,U\Big{|}_{y},$

and call $m^{x,y}_{ij}(\textbf{k},\textbf{k}^{\prime},\textbf{q})$ their
matrix elements, whose explicit expressions are given as follows

$\displaystyle m_{11}^{x}$
$\displaystyle=|\Omega_{0}(\textbf{q},t)|\Big{[}-E_{x}\sin
2\varphi(\textbf{k}-\textbf{q})\cos\delta_{x}(\textbf{k}^{\prime},\textbf{k}-\textbf{q})$
$\displaystyle+E_{y}\cos
2\varphi(\textbf{k}-\textbf{q})\cos\delta_{y}(\textbf{k}^{\prime},\textbf{k}-\textbf{q})\Big{]},$
(34a) $\displaystyle m_{12}^{x}$
$\displaystyle=-i|\Omega_{0}(\textbf{q},t)|\Big{[}E_{x}\sin
2\varphi(\textbf{k}-\textbf{q})\sin\delta_{x}(\textbf{k}^{\prime},\textbf{k}-\textbf{q})$
$\displaystyle-E_{y}\cos
2\varphi(\textbf{k}-\textbf{q})\sin\delta_{y}(\textbf{k}^{\prime},\textbf{k}-\textbf{q})\Big{]},$
(34b) $\displaystyle m_{11}^{y}$
$\displaystyle=|\Omega_{0}(\textbf{q},t)|\Big{[}E_{x}\cos
2\varphi(\textbf{k}-\textbf{q})\cos\delta_{x}(\textbf{k}^{\prime},\textbf{k}-\textbf{q})$
$\displaystyle+E_{y}\sin
2\varphi(\textbf{k}-\textbf{q})\cos\delta_{y}(\textbf{k}^{\prime},\textbf{k}-\textbf{q})\Big{]},$
(34c) $\displaystyle m_{12}^{y}$
$\displaystyle=i|\Omega_{0}(\textbf{q},t)|\Big{[}E_{x}\cos
2\varphi(\textbf{k}-\textbf{q})\sin\delta_{x}(\textbf{k}^{\prime},\textbf{k}-\textbf{q})$
$\displaystyle+E_{y}\sin
2\varphi(\textbf{k}-\textbf{q})\sin\delta_{y}(\textbf{k}^{\prime},\textbf{k}-\textbf{q})\Big{]},$
(34d)

where
$\delta_{k}(\textbf{k}^{\prime},\textbf{k}-\textbf{q})=\alpha_{k}+\frac{\varphi(\textbf{k}^{\prime})+\varphi(\textbf{k}-\textbf{q})}{2}$,
the current operator
$\hat{\textbf{j}}=\hat{j}_{x}\hat{\textbf{x}}+\hat{j}_{y}\hat{\textbf{y}}$ in
Eq. (31) can be then written as follows

$\displaystyle\hat{j}_{x}(t)$ $\displaystyle=$
$\displaystyle\sum_{k,q}(m_{11}^{x}\hat{e}^{\dagger}_{k^{\prime}}\hat{e}_{k-q}+m_{12}^{x}\hat{e}^{\dagger}_{k^{\prime}}\hat{h}^{\dagger}_{-k+q}$
$\displaystyle+$ $\displaystyle
m_{21}^{x}\hat{h}_{-k^{\prime}}\hat{e}_{k-q}+m_{22}^{x}\hat{h}_{-k^{\prime}}\hat{h}^{\dagger}_{-k+q}),$
$\displaystyle\hat{j}_{y}(t)$ $\displaystyle=$
$\displaystyle\sum_{k,q}(m_{11}^{y}\hat{e}^{\dagger}_{k^{\prime}}\hat{e}_{k-q}+m_{12}^{y}\hat{e}^{\dagger}_{k^{\prime}}\hat{h}^{\dagger}_{-k+q}$
$\displaystyle+$ $\displaystyle
m_{21}^{y}\hat{h}_{-k^{\prime}}\hat{e}_{k-q}+m_{22}^{y}\hat{h}_{-k^{\prime}}\hat{h}^{\dagger}_{-k+q}).$

Taking the expectation values we finally get

$\displaystyle j_{x}(t)$
$\displaystyle=\sum_{k^{\prime},k,q}\Bigg{\\{}-m_{11}^{x}(\textbf{k}^{\prime},\textbf{k},\textbf{q},t)\left[1-n_{e}(\textbf{k}^{\prime},\textbf{k}-\textbf{q})-n_{h}(\textbf{k}^{\prime},\textbf{k}-\textbf{q})\right]-2\operatorname{Im}\Big{\\{}\tilde{m}_{12}^{x}(\textbf{k}^{\prime},\textbf{k},\textbf{q},t)\,p(\textbf{k}^{\prime},\textbf{k}-\textbf{q})\Big{\\}}\Bigg{\\}},$
(36a) $\displaystyle j_{y}(t)$
$\displaystyle=\sum_{k^{\prime},k,q}\Bigg{\\{}-m_{11}^{y}(\textbf{k}^{\prime},\textbf{k},\textbf{q},t)\left[1-n_{e}(\textbf{k}^{\prime},\textbf{k}-\textbf{q})-n_{h}(\textbf{k}^{\prime},\textbf{k}-\textbf{q})\right]+2\operatorname{Im}\Big{\\{}\tilde{m}_{12}^{y}(\textbf{k}^{\prime},\textbf{k},\textbf{q},t)\,p(\textbf{k}^{\prime},\textbf{k}-\textbf{q})\Big{\\}}\Bigg{\\}},$
(36b)

where we have used the fact that $m_{11}^{x,y}=-m_{22}^{x,y}$ and
$m_{12}^{x,y}=i\tilde{m}_{12}^{x,y}=-i(\tilde{m}_{21}^{x,y})^{*}$.

Notice, that the current above has the same formal expression that the one for
the homogeneous case, once the replacements
$n_{e,h}(\textbf{k},\textbf{k}^{\prime})\rightarrow n_{e,h}(\textbf{k})$,
$p(\textbf{k},\textbf{k}^{\prime})\rightarrow p(\textbf{k})$, and
$m_{ij}^{x,y}(\textbf{k}^{\prime},\textbf{k},\textbf{q},t)\rightarrow\,m_{ij}^{x,y}(\textbf{k},\textbf{k},0,t)$
are made.

##  VI Nonlinear optical response

Figure 3: Nonlinear Optical Response $I(\omega)$ of graphene, as a function of
the different harmonics of the carrier frequency $\omega_{L}$ of the impinging
Laguerre-Gaussian electromagnetic pulse, for different values of the OAM
parameter $\ell$, namely $\ell=0$ (dotted, blue line), $\ell=1$ (dashed,
orange line), $\ell=2$ (dotted, red line), and $\ell=3$ (dashed, green line).
To generate this figure, the following parameters have been used:
$E_{0}=10^{7}$ V/m, $\tau=50$ fs, $\omega_{L}=1$ THz. Moreover, $z=0$ and
$w_{0}=1$ (corresponding to scaled, transverse coordinates) have been assumed.

To check the validity of our approach, in this section we calculate the
nonlinear optical response of graphene, when excited by a structured light
pulse. To this aim, and without loss of generality, we assume the impinging
pulse to be polarised along the $x$-direction, and we choose its temporal
profile of the pulse to be Gaussian, while its spatial intensity distribution
matches that of a Laguerre-Gaussian (LG) beam,i.e.,

$\textbf{E}(\textbf{r},t)=E_{0}e^{-\frac{(t-t_{0})^{2}}{\tau^{2}}-\frac{R^{2}}{w^{2}(z)}}u(\textbf{R},z)\cos\omega_{L}t\,\hat{\textbf{x}},$
(37)

where $\tau$ is the pulse duration, $t_{0}$ an arbitrary offset,
$\textbf{R}=\\{R,\phi\\}$, and

$u(x,y,z)=\frac{C_{\ell,p}}{w(z)}\left(\frac{R\sqrt{2}}{w(z)}\right)^{|\ell|}\text{L}_{p}^{|\ell|}\left(\frac{2R^{2}}{w^{2}(z)}\right)e^{i\ell\phi},$
(38)

with $w(z)=w_{0}\sqrt{1+(z/z_{R})^{2}}$ is the Gaussian beam waist (with
$z_{R}$ being the Rayleigh range of the LG-beam),
$C_{\ell,p}=w_{0}\sqrt{2p!/[\pi(p+|\ell|)!]}$, and $\text{L}_{p}^{|\ell|}(x)$
are the associated Laguerre polynomials DLMF , and $\ell$ represents the OAM
index, while $p$ represents the radial index of the LG-beam.

Figure 4: Nonlinear Optical Response $I(\omega)$ of graphene as a function of
the different harmonics of the carrier frequency $\omega_{L}$ of the impinging
Laguerre-Gaussian electromagnetic pulse, for different values of the OAM
parameter $\ell$, namely $\ell=0$ (dashed, blue line), $\ell=1$ (dotted,
orange line), $\ell=2$, (dashed, red line), and $\ell=3$ (dotted, green line).
The parameters for this simulations are the same as those used for generating
Fig. 3. Moreover, this figure has been calculated under the assumption
$\textbf{k}=\textbf{k}^{\prime}$. Notice, how this assumption essentially
removed any dependence on the transverse intensity distribution, levelling all
the different OAM modes to the same result obtained for $\ell=0$.

We then solve Eqs. (29) for different values of the OAM parameter $\ell$ (and
set $p=0$, for simplicity), and then use Eqs. (36) to compute the
correspondent nonlinear current $\textbf{j}(t)$ by integrating over the
transverse shape of the pulse, and we then calculate the nonlinear signal as
$I(\omega)\thicksim|\omega\,\mathbf{J}(\omega)|^{2}$, where
$\mathbf{J}(\omega)$ is the Fourier transform of the nonlinear current.

The resulting nonlinear spectrum for different values of the OAM parameter
$\ell$ is shown in Fig. 3. At a first glance, this spectrum shares the same
characteristics of the usual harmonic spectrum obtained for graphene, i.e.,
the odd harmonics are predominant in magnitude than the even one, since the
even ones are suppressed because of symmetry reasons Katsnelson (2007).

If we take a closer look at the even harmonics, however, we can see that they
are not completely suppressed, as it happens for the case of graphene [see
Fig. 4]. This is an indirect result of the fact, that the intensity of the
impinging pulse was not uniform overall the graphene sample, but had the
characteristic doughnut shape of LG modes. The emergence of even harmonics in
the nonlinear spectrum of a monolayer of graphene has been reported before
Carvalho  _et al._ and it has been explained in terms of a pulse-induced
oscillation of the Dirac cone, that ultimately amplifies even harmonics. In
the context of this manuscript, we observe the same effect, but triggered by
the non-uniform spatial intensity distribution of the pulse, rather than its
temporal features. To corroborate this statement, we can compare the results
in Fig. 3, obtained by solving the full Dirac-Bloch system [Eqs. (29)], with
Fig. 4, which essentially corresponds to the uniform intensity case (see
below). As it can be seen, while the even harmonics are completely suppressed
in the uniform illumination case, for the case of structured illumination
reported in Fig. 3, the actual non-uniform spatial intensity distribution has
the overall effect to slightly break graphene’s centrosymmetry, thus allowing
the generation of even harmonics.

###  VI.1 A Remark on the numerical model

It is worth noticing, that to obtain the result shown in Fig. 3, we have used
the expression of the current given by Eq. (31), which is a time-dependent
only quantity, since it is the result of the integration, over the
4-dimensional reciprocal space $\\{\textbf{k},\textbf{k}^{\prime}\\}$, of the
$k$-dependent current density. If we want to have information on the spatial
variation of the current or, analogously, of the nonlinear spectrum, we would
instead need to compute the spatial Fourier transform of the current density,
instead of simply integrating. By doing so, we will obtain the following,
nonlocal, expression for the spatially-dependent nonlinear current

$\displaystyle j_{x}(\textbf{R},\textbf{R}^{\prime},t)$
$\displaystyle=\sum_{k,k^{\prime},q}\mathcal{J}_{x}(\textbf{k},\textbf{k}^{\prime},\textbf{q},t)e^{i(\textbf{k}\cdot\textbf{R}+\textbf{k}^{\prime}\cdot\textbf{R}^{\prime})},$
(39a) $\displaystyle j_{y}(\textbf{R},\textbf{R}^{\prime},t)$
$\displaystyle=\sum_{k,k^{\prime},q}\mathcal{J}_{y}(\textbf{k},\textbf{k}^{\prime},\textbf{q},t)e^{i(\textbf{k}\cdot\textbf{R}+\textbf{k}^{\prime}\cdot\textbf{R}^{\prime})},$
(39b)

where $\mathcal{J}_{x,y}(\textbf{k},\textbf{k}^{\prime},t)$ are the
expressions inside the curly brackets in Eqs. (36). Calculating the Fourier
transforms above, however, is computationally very demanding, since the
current densities $\mathcal{J}_{x,y}(\textbf{k},\textbf{k}^{\prime},t)$ are
defined on an 8-dimensional grid (since they are functions of three
independent $k$-vectors, namely k, k’, and q), which evolves in time. To
efficiently perform these calculations, one would probably need to solve Eqs.
(29) in a more efficient way, optimising the code and using, for example, a
supercomputer, and then be very careful in defining a big enough grid in both
k and $\textbf{k}^{\prime}$ to ensure enough resolution for the Fourier
transform. Doing this, however, would require a considerable amount of time
and it is outside the scope of this work.

To partly solve this issue, however, we suggest a simplification of the
problem, which results in significantly reducing the computational complexity
of the problem, and could be employed to make qualitative predictions on the
spatial distribution of the current. To do that, we restrict the 4-dimensional
$k$-space to an effective 2-dimensional one by setting
$\textbf{k}=\textbf{k}^{\prime}$ in Eqs. (29). This, essentially, means, that
we restrict the domain of definition of $\textbf{k}^{\prime}$ to a point,
rather than a plane. To test the validity of this approximation, in Fig. 4 we
again show the harmonic spectrum and notice, that we get a pretty similar, at
least qualitatively, result, than the one obtained in fig. 3 using the full
model. Notice, in particular, how setting $\textbf{k}=\textbf{k}^{\prime}$
essentially erases any information about the even harmonics. This, again, is
not surprising, since the limit $\textbf{k}=\textbf{k}^{\prime}$ can be
interpreted, to some extent, as the homogeneous field limit.

Under this approximation, however, we would be able to look at the qualitative
spatial evolution of the current, as now the components of the spatially
dependent nonlinear current are much easier to calculate, as their explicit
expressions now only amount to a regular spatial Fourier transform, i.e.,

$\displaystyle j_{x}(\textbf{R},t)$
$\displaystyle=\sum_{k,q}\mathcal{J}_{x}(\textbf{k},\textbf{q},t)e^{i\textbf{k}\cdot\textbf{R}},$
(40a) $\displaystyle j_{y}(\textbf{R},t)$
$\displaystyle=\sum_{k,q}\mathcal{J}_{y}(\textbf{k},\textbf{q},t)e^{i\textbf{k}\cdot\textbf{R}}.$
(40b)

Despite the equations above are greatly simplified, in terms of computational
challenge, with respect to Eqs. (39), they still need quite a high resolution
on the $k$-space grid, in order to be able to generate meaningful images in
direct space. Although the implementation of such optimised code is out of the
scope of this manuscript, we intend to work on it for a future publication,
where we aim at showing the different current dynamics corresponding to
different spatial profiles of the impinging pulse.

##  VII Conclusion

In this work, we introduced a general framework to study the optical response
of graphene under the action of spatially-varying electromagnetic pulses,
i.e., when graphene interacts with structured light. In particular, we have
provided detailed derivation of the spatially-varying interaction Hamiltonian,
as well as the two-body population and polarisation, constituting the building
blocks for the Dirac-Bloch equations. We have then calculated the expression
of the induced nonlinear current and used it to test our model against known
results, by calculating the nonlinear spectrum for the case of an impinging
Laguerre-Gaussian optical pulse. Surprisingly, although the overall structure
of the nonlinear spectrum is not changed by a spatially varying
electromagnetic field, we have reported how a spatially inhomogeneous pulse
can locally break the centrosymmetry of graphene and allow the generation of
even harmonics. Finally, we have shown how to obtain qualitative results to
study the spatial evolution of the current density.

Although in this work we have focused our attention to graphene, our results
can be easily adapted to any 2D material and provide a general framework to
include the effects of structured light in both the linear and nonlinear
optical response of such materials.

## Acknowledgements

We acknowledge the financial support of the Academy of Finland Flagship
Programme, Photonics research and innovation (PREIN), decision 320165. Y.T.
also acknowledges support from the Finnish Cultural Foundation, decision
00221008.

## References

  * Novoselov  _et al._ (2004) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
  * Katsnelson (2006) M. Katsnelson, The European Physical Journal B-Condensed Matter and Complex Systems 51, 157 (2006).
  * (3) R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. R. Peres, and A. K. Geim,  320, 1308.
  * Novoselov  _et al._ (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005).
  * (5) .
  * Wright  _et al._ (2009) A. R. Wright, X. G. Xu, J. C. Cao, and C. Zhang, Applied Physics Letters 95, 072101 (2009).
  * Hafez  _et al._ (2020) H. A. Hafez, S. Kovalev, K. Tielrooij, M. Bonn, M. Gensch, and D. Turchinovich, Advanced Optical Materials 8, 1900771 (2020).
  * Manzeli  _et al._ (2017) S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, Nature Reviews Materials 2, 17033 (2017).
  * Autere  _et al._ (2018) A. Autere, H. Jussila, Y. Dai, Y. Wang, H. Lipsanen, and Z. Sun, Advanced Materials 30, 1705963 (2018).
  * Neto  _et al._ (2009) A. C. Neto, F. Guinea, N. M. Peres, K. S. Novoselov, and A. K. Geim, Reviews of modern physics 81, 109 (2009).
  * Caldwell  _et al._ (2019) J. D. Caldwell, I. Aharonovich, G. Cassabois, J. H. Edgar, B. Gil, and D. N. Basov, Nature Reviews Materials 4, 552 (2019).
  * (12) A. K. Geim and K. S. Novoselov,  6, 183.
  * Zeng  _et al._ (2013) H. Zeng, G.-B. Liu, J. Dai, Y. Yan, B. Zhu, R. He, L. Xie, S. Xu, X. Chen, W. Yao, _et al._ , Scientific reports 3, 1 (2013).
  * Kumar  _et al._ (2013) N. Kumar, J. Kumar, C. Gerstenkorn, R. Wang, H.-Y. Chiu, A. L. Smirl, and H. Zhao, Physical Review B 87, 121406 (2013).
  * Hendry  _et al._ (2010) E. Hendry, P. J. Hale, J. Moger, A. K. Savchenko, and S. A. Mikhailov, Physical Review Letters 105, 097401 (2010).
  * Dong  _et al._ (2015) N. Dong, Y. Li, Y. Feng, S. Zhang, X. Zhang, C. Chang, J. Fan, L. Zhang, and J. Wang, Scientific Reports 5, 14646 (2015).
  * Ornigotti  _et al._ (2021) M. Ornigotti, L. Ornigotti, and F. Biancalana, APL Photonics 6, 060801 (2021).
  * (18) D. N. Carvalho, F. Biancalana, and A. Marini,  3, 10.1515/odps-2017-0006.
  * Allen  _et al._ (1992) L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, Physical Review A 45, 8185 (1992).
  * (20) S. Franke-Arnold, L. Allen, and M. Padgett,  2, 299.
  * (21) R. Fickler, R. Lapkiewicz, W. N. Plick, M. Krenn, C. Schaeff, S. Ramelow, and A. Zeilinger,  338, 640.
  * (22) T. A. Klar and S. W. Hell,  24, 954.
  * (23) M. A. Noyan and J. M. Kikkawa,  107, 032406.
  * (24) J. Wang, J.-Y. Yang, I. M. Fazal, N. Ahmed, Y. Yan, H. Huang, Y. Ren, Y. Yue, S. Dolinar, M. Tur, and A. E. Willner,  6, 488.
  * (25) D. Cozzolino, D. Bacco, B. Da Lio, K. Ingerslev, Y. Ding, K. Dalgaard, P. Kristensen, M. Galili, K. Rottwitt, S. Ramachandran, and L. K. Oxenløwe,  11, 064058.
  * Wätzel and Berakdar (2016) J. Wätzel and J. Berakdar, Scientific Reports 6, 21475 (2016).
  * Farías  _et al._ (2013) M. B. Farías, G. F. Quinteiro, and P. I. Tamborenea, The European Physical Journal B 86, 432 (2013).
  * (28) K. L. Ishikawa,  82, 201402.
  * (29) M. I. Katsnelson,  _Graphene: Carbon in Two Dimensions_ (Cambridge University Press).
  * (30) S. Barnett and P. Radmore,  _Methods in Theoretical Quantum Optics_ (Oxford University Press).
  * Haug and Koch (2009) H. Haug and S. W. Koch,  _Quantum theory of the optical and electronic properties of semiconductors_ (World Scientific Publishing Company, 2009).
  * Katsnelson (2007) M. I. Katsnelson, Materials today 10, 20 (2007).
  * Kuklinski and Mukamel (1991) J. Kuklinski and S. Mukamel, Physical Review B 44, 11253 (1991).
  * Kelardeh  _et al._ (2015) H. K. Kelardeh, V. Apalkov, and M. I. Stockman, Physical Review B 91, 045439 (2015).
  * (35) L. Mandel and E. Wolf,  _Optical Coherence and Quantum Optics_, 1st ed. (Cambridge University Press).
  * Stroucken  _et al._ (2011) T. Stroucken, J. Grönqvist, and S. Koch, Physical Review B 84, 205445 (2011).
  * (37) DLMF, “NIST Digital Library of Mathematical Functions,” http://dlmf.nist.gov/, Release 1.1.5 of 2022-03-15, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.

.


Some application of difference equations in Cryptography and Coding Theory

Cristina FLAUT

Faculty of Mathematics and Computer Science, Ovidius University,

Bd. Mamaia 124, 900527, CONSTANTA, ROMANIA

http://www.univ-ovidius.ro/math/

e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
Abstract. In this paper, we present some applications of a difference equation
of degree $k$ in Cryptography and Coding Theory.

Key Words. Fibonacci numbers; difference equations; encrypting and decrypting;
coding and decoding.

2000 AMS Subject Classification: 11B39, 94Bxx.

1\. Introduction

There are many papers devoted to the study of the properties and the
applications of some particular integer sequences, as for example: Fibonacci
sequences, $p-$Fibonacci sequences, Tribonacci sequences, etc. (see [1], [2],
[4], [5], [8], [9], [10], [11]). In this paper, we generalize these results,
by considering the general case of a difference equation of degree $k$, we
associate a matrix to such an equation and, using some properties of these
matrices, we give some applications of them in Cryptography and Coding Theory.
We generalize the notion of complete positive integers sequence, given in [6],
and a result given in [13], result which states the representation of a
natural number as a sum of nonconsecutive Fibonacci numbers. With these
results, we give an algorithm for messages encryption and decryption and, in
Section 3, we give an application in Coding Theory. In Appendix, we present
some MAPLE procedures. These procedures are very helpful in the encrypting and
decrypting processes.

2\. Some properties of a difference equation of degree $k,k\geq 2$

Let $n$ be an arbitrary positive integer and $a,b$ be arbitrary integers,
$b\neq 0$. We consider the following difference equation of degree two

$d_{n}=ad_{n-1}+bd_{n-2},d_{0}=0,d_{1}=1$ (2.1.)

and the attached matrix

$D_{2}=\left(\begin{array}[]{cc}a&b\\\
1&0\end{array}\right)=\left(\begin{array}[]{cc}d_{2}&bd_{1}\\\
d_{1}&bd_{0}\end{array}\right).$

It results

$D_{2}^{2}=\left(\begin{array}[]{cc}b+a^{2}&ab\\\
a&b\end{array}\right)=\left(\begin{array}[]{cc}d_{3}&bd_{2}\\\
d_{2}&bd_{1}\end{array}\right).$

Therefore, we obtain that

$D_{2}^{n}=\left(\begin{array}[]{cc}d_{n+1}&bd_{n}\\\
d_{n}&bd_{n-1}\end{array}\right).$ (2.2.)

Since det$D_{2}=-b$ and
det$D_{2}^{n}=bd_{n-1}d_{n+1}-bd_{n}^{2}=\left(-b\right)^{n},$ the following
relation is true:

$d_{n-1}d_{n+1}-d_{n}^{2}=\left(-b\right)^{n-1}.$ (2.3.)

The inverse of the matrix $D_{2}^{n}$ is

$D_{2}^{-n}==\frac{1}{\left(-1\right)^{n+1}b^{n}}\left(\begin{array}[]{cc}-bd_{n-1}&bd_{n}\\\
d_{n}&-d_{n+1}\end{array}\right).$ (2.4.)

If we consider the following recurrence relation of degree three

$d_{n}=ad_{n-1}+bd_{n-2}+cd_{n-2},d_{-1}=d_{0}=d_{1}=0,d_{2}=1,c\neq 0,$
(2.5.)

we have the attached matrix

$D_{3}=\left(\begin{array}[]{ccc}d_{3}&bd_{2}+cd_{1}&cd_{2}\\\
d_{2}&bd_{1}+cd_{0}&cd_{1}\\\
d_{1}&bd_{0}+cd_{-1}&cd_{0}\end{array}\right)=\left(\begin{array}[]{ccc}a&b&c\\\
1&0&0\\\ 0&1&0\end{array}\right),$

since we can take $bd_{0}+cd_{-1}=d_{2}-ad_{1}=1.$ Therefore, using relation
$d_{4}=a^{2}+b,$ we get

$D_{3}^{2}=\allowbreak\left(\begin{array}[]{ccc}b+a^{2}&c+ab&ac\\\ a&b&c\\\
1&0&0\end{array}\right)=\allowbreak\left(\begin{array}[]{ccc}d_{4}&bd_{3}+cd_{2}&cd_{3}\\\
d_{3}&bd_{2}+cd_{1}&cd_{2}\\\ d_{2}&bd_{1}+cd_{0}&cd_{1}\end{array}\right).$

From here, it results

$D_{3}^{n}=\left(\begin{array}[]{ccc}d_{n+2}&bd_{n+1}+cd_{n}&cd_{n+1}\\\
d_{n+1}&bd_{n}+cd_{n-1}&cd_{n}\\\
d_{n}&bd_{n-1}+cd_{n-2}&cd_{n-1}\end{array}\right).$ (2.6.)

Since det$D_{3}=c$ and det$D_{3}^{n}=$
$=c^{2}\left[d_{n}\left(d_{n}^{2}-d_{n-1}d_{n+1}\right)+d_{n-2}\left(d_{n+1}^{2}-d_{n}d_{n+2}\right)+d_{n-1}\left(\allowbreak
d_{n-1}d_{n+2}-d_{n}d_{n+1}\right)\right],$ using the fact that
det$D_{3}^{n}=c^{n},n\geq 2$, we obtain the following relation

$d_{n}\left(d_{n}^{2}\text{-}d_{n-1}d_{n+1}\right)\text{+}d_{n-2}\left(d_{n+1}^{2}\text{-}d_{n}d_{n+2}\right)\text{+}d_{n-1}\left(\allowbreak
d_{n-1}d_{n+2}\text{-}d_{n}d_{n+1}\right)\text{=}c^{n-2}.$ (2.7.)

The inverse of the matrix $D_{3}^{n}$ is

$D_{3}^{-n}=\frac{1}{c^{n-2}}\left(\begin{array}[]{ccc}g_{11}&g_{12}&g_{13}\\\
g_{21}&g_{22}&g_{23}\\\ g_{31}&g_{32}&g_{33}\end{array}\right),$ (2.8.)

where
$g_{11}=$ $c^{2}(-d_{n}d_{n-2}+d_{n-1}^{2})\allowbreak$,
$g_{21}=-c(-d_{n}^{2}+d_{n-1}d_{n+1}),$
$g_{31}=-(bd_{n}^{2}+cd_{n}d_{n-1}-bd_{n-1}d_{n+1}-cd_{n+1}d_{n-2}),$
$g_{12}=-c^{2}(d_{n}d_{n-1}-d_{n+1}d_{n-2}),g_{22}=c(-d_{n}d_{n+1}+d_{n-1}d_{n+2}),$
$g_{32}=cd_{n}^{2}+bd_{n}d_{n+1}-bd_{n-1}d_{n+2}-cd_{n-2}d_{n+2},$
$g_{13}=c^{2}(d_{n}^{2}-d_{n-1}d_{n+1}),g_{23}=c(-d_{n}d_{n+2}+d_{n+1}^{2}),$
$g_{33}=bd_{n}d_{n+2}-cd_{n}d_{n+1}-bd_{n+1}^{2}+cd_{n-1}d_{n+2}.$

Now, we consider the general $k-$terms recurrence, $n,k\in\mathbb{N},k\geq
2,n\geq k,$

$d_{n}\text{=}a_{1}d_{n-1}\text{+}a_{2}d_{n-2}\text{+}...\text{+}a_{k}d_{n-k},d_{0}\text{=}d_{1}\text{=...=}d_{k-2}\text{=}0,d_{k-1}\text{=}1,a_{k}\neq
0$ (2.9.)

and the matrix $D_{k}\in\mathcal{M}_{k}\left(\mathbb{R}\right),$

$D_{k}=\left(\begin{array}[]{ccccc}a_{1}&a_{2}&a_{3}&...&a_{k}\\\
1&0&0&...&0\\\ 0&1&0&...&0\\\ ...&...&...&...&...\\\
0&0&...&1&0\end{array}\right),$ (2.10.)

(see [7]).

Proposition 2.1. With the above notations, the following relations are true:
1)

$D_{k}=\
\left(\begin{array}[]{ccccc}d_{k}&\overset{k-1}{\underset{i=1}{\sum}}a_{i+1}d_{k-i}&\overset{k-2}{\underset{i=1}{\sum}}a_{i+2}d_{k-i}&...&a_{k}d_{k-1}\\\
d_{k-1}&\overset{k-1}{\underset{i=1}{\sum}}a_{i+1}d_{k-i-1}&\overset{k-2}{\underset{i=1}{\sum}}a_{i+2}d_{k-i-1}&...&a_{k}d_{k-2}\\\
d_{k-2}&\overset{k-1}{\underset{i=1}{\sum}}a_{i+1}d_{k-i-2}&\overset{k-2}{\underset{i=1}{\sum}}a_{i+2}d_{k-i-2}&...&a_{k}d_{k-3}\\\
...&...&...&...&...\\\
d_{1}&\overset{k-1}{\underset{i=1}{\sum}}a_{i+1}d_{k-i-k+1}&\overset{k-2}{\underset{i=1}{\sum}}a_{i+2}d_{-i+1}&...&a_{k}d_{0}\end{array}\right).$
(2.11.)

2) For $n\in\mathbb{Z},n\geq 1,$ we have that

$D_{k}^{n}\text{=}\left(\begin{array}[]{ccccc}d_{n+k-1}&\overset{k-1}{\underset{i=1}{\sum}}a_{i+1}d_{n+k-i-1}&\overset{k-2}{\underset{i=1}{\sum}}a_{i+2}d_{n+k-i-1}&...&a_{k}d_{n+k-2}\\\
d_{n+k-2}&\overset{k-1}{\underset{i=1}{\sum}}a_{i+1}d_{n+k-i-2}&\overset{k-2}{\underset{i=1}{\sum}}a_{i+2}d_{n+k-i-2}&...&a_{k}d_{n+k-3}\\\
d_{n+k-3}&\overset{k-1}{\underset{i=1}{\sum}}a_{i+1}d_{n+k-i-3}&\overset{k-2}{\underset{i=1}{\sum}}a_{i+2}d_{n+k-i-3}&...&a_{k}d_{n+k-4}\\\
...&...&...&...&...\\\
d_{n}&\overset{k-1}{\underset{i=1}{\sum}}a_{i+1}d_{n+k-i-k}&\overset{k-2}{\underset{i=1}{\sum}}a_{i+2}d_{n-i}&...&a_{k}d_{n-1}\end{array}\right).$
(2.12.)

Proof. 1) It is obvious, using relation $\left(2.9\right)$. Indeed,
$d_{k}=a_{1}.$ Making computations, we obtain

$\overset{k-1}{\underset{i=1}{\sum}}a_{i+1}d_{k-i}=a_{2}d_{k-1}+a_{3}d_{k-2}+...+a_{k}d_{1}=a_{2},$
$...$
$\overset{k-1}{\underset{i=1}{\sum}}a_{i+1}d_{k-i-2}\text{=}a_{2}d_{k-3}+a_{3}d_{k-4}+...+a_{k}d_{-1}\text{=}d_{k-1}-a_{1}d_{k-2}=1,$

and so on.

2) We use induction. For $n=1$, the relation is true. Let
$D_{k}^{n+1}=\left(h_{ij}\right)_{i,j\in\\{1,2,...,k\\}}$. Since
$D_{k}^{n+1}=D_{k}^{n}D_{k}$, assuming that the statement is true for $n$ and
using relation $\left(2.10\right),~{}$it results the following elements for
the matrix $D_{k}^{n+1}:$

$h_{11}\text{=}a_{1}d_{n+k-1}\text{+}\overset{k-1}{\underset{i=1}{\sum}}a_{i+1}d_{n+k-i-1}\text{=}a_{1}d_{n+k-1}+a_{2}d_{n+k-2}+...+a_{k}d_{n}\text{=}d_{n+k}\newline
,$
$h_{12}\text{=}a_{2}d_{n+k-1}+\overset{k-2}{\underset{i=1}{\sum}}a_{i+2}d_{n+k-i-1}\text{=}a_{2}d_{n+k-1}\text{+}a_{3}d_{n+k-2}\text{+...+}a_{k}d_{n-1}\text{=}\overset{k-1}{\underset{i=1}{\sum}}a_{i+1}d_{n+k-i},$

$...$ $h_{1k}=a_{k}d_{n+k-1},$
$h_{21}\text{=}a_{1}d_{n+k-2}\text{+}\overset{k-1}{\underset{i=1}{\sum}}a_{i+1}d_{n+k-i-2}\text{=}a_{1}d_{n+k-2}+a_{2}d_{n+k-3}+...+a_{k}d_{n-1}\text{=}d_{n+k-1},$
$h_{22}\text{=}a_{2}d_{n+k-2}\text{+}\overset{k-2}{\underset{i=1}{\sum}}a_{i+2}d_{n+k-i-2}\text{=}a_{2}d_{n+k-2}\text{+}a_{3}d_{n+k-3}\text{+...+}a_{k}d_{n}\text{=}\overset{k-1}{\underset{i=1}{\sum}}a_{i+2}d_{n+k-i-1},$

etc.$~{}\Box\vskip 6.0pt plus 2.0pt minus 2.0pt$

Remark 2.2. ([7]) The following statements are true:
i) det$D_{k}=\left(-1\right)^{k+1}a_{k}$ and

$D_{k}^{-1}=\frac{1}{a_{k}}\left(\begin{array}[]{ccccc}0&a_{k}&0&...&0\\\
1&0&a_{k}&...&0\\\ ...&...&...&...&...\\\ 0&0&0&...&1\\\
1&-a_{1}&...&-a_{k-2}&-a_{k-1}\end{array}\right).$

ii) $D_{k}^{n}=a_{1}D_{k}^{n-1}+a_{2}D_{k}^{n-2}+...+a_{k}D_{k}^{n-k}.\vskip
6.0pt plus 2.0pt minus 2.0pt$

Remark 2.3. From the above, we have that

$\text{{det}}D_{k}^{n}=\left(-1\right)^{\left(k+1\right)n}a_{k}^{n}.$ (2.13.)

Relations 2.3, 2.7 and 2.13 are Cassini’s type relations.

3\. Applications in cryptography and coding theory

A sequence of positive integers, $\left(d_{n}\right),n\geq 0,$ is complete if
and only if each natural number $n$ can be written under the form
$\underset{i=1}{\overset{m}{\sum}}c_{i}d_{i},$ where $c_{i}$ is either zero or
one (see [6]). In [3], the author proved that a nondecreasing sequence of
natural numbers, with $d_{1}=1,$ is complete if and only if the following
relation

$d_{k+1}\leq 1+\underset{i=1}{\overset{k}{\sum}}d_{i},$

is true for $k\in\\{1,2,3,.....\\}.$

There are sequences of positive integers which are complete and sequences
which are not complete. An example of complete sequence is the sequence of
Fibonacci numbers. In [13], the author proved that each natural number can be
written as a unique sum of non-consecutive Fibonacci numbers, therefore the
sequence of Fibonacci numbers is complete. An example of not complete sequence
was done in [12]. The sequence $\left(d_{n}\right),n\geq 0,$ where
$d_{n}=4d_{n-1}+3d_{n-2},d_{0}=0,d_{1}=1,$ is not complete.

In the following, we will generalize the notion of complete sequence and we
will give applications of this new notion.

Definition 3.1. A sequence of positive integers, $\left(d_{n}\right),n\geq 0,$
is called general complete (or g-complete) if and only if for each natural
number $n\geq 0,$ there is a natural number $q,$ $q\geq 0,$ such that $n$ can
be written, in a unique way, under the form

$n=d_{0}a_{0}+a_{1}d_{1}+...+a_{q}d_{q},a_{0},...a_{q}\in\mathbb{N},a_{q}\neq
0.$

Theorem 3.2. The sequence of positive integers $\left(d_{n}\right),n\geq 0,$
generated by the difference equation $\left(2.9\right)$ is g-complete.

Proof. We consider the sequence of positive integers
$\left(d_{n}\right)_{n\geq 0}$ generated by the difference equation given in
$\left(2.9\right)$. Let $n,q$ be the natural numbers, with $q\geq k\geq
2,~{}$such that the term $d_{q},$ of the sequence $\left(d_{n}\right)_{n\geq
0},$ satisfies the condition $d_{q}\leq n<d_{q+1}.$ We use the Quotient
Remainder Theorem. Therefore, we obtain the natural numbers
$c_{q},c_{q-1},....c_{k-1}$ and $r_{q},r_{q-1},...r_{k}$ such that

$n=d_{q}c_{q}+r_{q},0\leq r_{q}<d_{q},\newline $
$r_{q}=d_{q-1}c_{q-1}+r_{q-1},0\leq r_{q-1}<d_{q-1},\newline $
$r_{q-1}=d_{q-2}c_{q-2}+r_{q-2},0\leq r_{q-2}<d_{q-2},\newline $
$........\newline $
$r_{k+1}=d_{k}c_{k}+r_{k},0\leq r_{k}<d_{k},\newline $ $r_{k}=d_{k-1}c_{k-1},$

with $c_{k-1}=r_{k},$ since $d_{k-1}=1.$
It results that the number $n$ can be written, in a unique way, under the form

$n=d_{q}c_{q}+d_{q-1}c_{q-1}+d_{q-2}c_{q-2}+...+d_{k}c_{k}+d_{k-1}c_{k-1},$
(3.1.)

Indeed, we can’t have $k$ terms in the relation $\left(3.1\right),$
$d_{s},d_{s+1},...d_{s+k-1},$ with $s+k-1\leq q$ such that
$c_{s}d_{s}+c_{s+1}d_{s+1}+...+c_{s+k-1}d_{s+k-1}=d_{s+k}.$ In relation
$\left(3.1\right)$, if $s+k<q,$ then the coefficient of $d_{s+k}$ is
$c_{s+k}+1$. From here, we get

$r_{s+k+1}=\left(c_{s+k}+1\right)d_{s+k}+r_{s+k}-d_{s+k}.$

Since $r_{s+k}-$ $d_{s+k}<0,$ we obtain a contradiction with the Quotient
Remainder Theorem.$~{}\Box\vskip 6.0pt plus 2.0pt minus 2.0pt$

The above Theorem extend the result obtained in [13] for Fibonacci numbers to
a difference equation of degree $k\,$, equation defined by $\left(2.9\right)$.

3.1. An application in Cryptography

In [8], [10], [11], [12], were described applications of Fibonacci numbers,
$k-$Fibonacci numbers or Tribonacci numbers in Cryptography. In the following,
by using Theorem 3.2 and the terms $\left(d_{n}\right)_{n\geq 0}$ of the
difference equation given by the relation $\left(2.9\right)$, we will present
a new method for encrypting and dencrypting messages. This new method give us
a lot of possibilities for finding the encryption and decryption keys, having
the advantage that each natural number $n$ has a unique representation by
using the terms of such a sequence. It results that a number cannot have the
same encrypted value as another number, therefore the obtained encrypted texts
are hard to break.

The Algorithm

Let $\mathcal{A}$ be an alphabet with $N$ letters, labeled from $0$ to $N-1,$
$m$ be a plain text and $n$ be a number obtained by using the label of the
letters from the plain text $m.$ We split the text $m$ in blocks with the same
length, $m_{1},m_{2},...,m_{r}.$ To these blocks correspond the numbers
$n_{1},n_{2},...,n_{r}.$ The numbers $n_{i},i\in\\{1,2,...,r\\},$ will be
encrypted using the following procedure.

The encrypting

1\. We consider the natural numbers $k\geq 2,a_{1},a_{2},...,a_{k};$

2\. We consider the difference equation of degree $k,$ given by the relation
$\left(2.9\right),$ with the above coefficients $a_{1},a_{2},...,a_{k}.$

3\. We compute the elements $d_{k},d_{k+1},...,d_{q},~{}$such that $d_{q}\leq
n_{i}<d_{q+1},i\in\\{1,2,...,r\\},$ as in the proof of Theorem 3.2. We denote
$s=q-k+2.$

4\. We obtain the secret encryption and decryption key
$\left(a_{1},a_{2},...,a_{k},s\right).$

5\. We compute the numbers $c_{iq},c_{i(q-1)},...,c_{i(k-1)}$ such that

$n_{i}=d_{q}c_{iq}+d_{q-1}c_{i(q-1)}+d_{q-2}c_{i(q-2)}+...+d_{k}c_{ik}+d_{k-1}c_{i(k-1)},$

as in relation $\left(3.1\right).$

6\. We obtain $\left(c_{iq},c_{i(q-1)},...,c_{i(k-1)}\right),$ the labels for
the encrypted text. If all $c_{ij}\leq N-1,$ then we will put in the cipher
text the letter $L\,_{j},$ from the alphabet $\mathcal{A},$ which has the
label $c_{ij}.$ If there are $c_{ij}>N-1,$ then we will put $c_{ij}~{}$in the
cipher text. If there are more than one $c_{ij}$ such that $c_{ij}>N-1,$
namely $\\{c_{ij_{1}},...c_{ij_{t}}\\},~{}$we will count the biggest number of
decimals of these elements. Assuming that $c_{ij_{v}}$ has the biggest number
of decimals, $0\leq v\leq t,$ we will put zeros to the left side of the
$c_{ij_{p}},p\neq v,$ such that all new obtained $c_{ij_{l}}~{}$have the same
number of decimals. In this way, we obtain the encrypted text $T_{i}$.

7\. Using the above steps for all blocks, we obtain the cipher text, by joint
the cipher texts $T_{1}T_{2}...T_{r}$.

The decrypting

1\. We use the key $\left(a_{1},a_{2},...,a_{k},s\right).$ The last number $s$
is the length of encoded blocks.

2\. We split the cipher text in blocks of length $s$ and we consider their
labels $\left(c_{iq},c_{i(q-1)},...,c_{i(k-1)}\right).~{}$We compute the
elements $d_{k},d_{k+1},...d_{q}$, therefore the corresponding plain text for
each block is
$n_{i}=d_{q}c_{iq}+d_{q-1}c_{i(q-1)}+d_{q-2}c_{i(q-2)}+...+d_{k}c_{ik}+d_{k-1}c_{i(k-1)}.$

3\. If there are numbers $\\{c_{ij_{1}},...c_{ij_{t}}\\}$ which are greater
than$~{}N-1,$ these numbers appear with the same number of digits or with new
labels, therefore it is easy to find the real message.

Example 3.3. 1) We consider an alphabet with 27 letters: A, B, C,…..Z, x,
where ”x” represent the blank space, labeled with $0,1,2,...,26.$ We want to
encrypt the message ”JOHNxHASxAxDOG”. We split this message in the following
blocks: JOHN, xHAS, xAxD, OGxx. We added two characters at the end of the last
block to obtain a block with $4$ letters. We consider a difference equation of
degree four, therefore $k=4,$ and $a_{1}=18,a_{2}=10,a_{3}=13,a_{4}=3$.
Therefore, the encoding key is SKND, or 18101303\. If, for example
$a_{1}=182,$ the key will be 182010013003.

To the first block, JOHN, will correspond the number $n_{1}=9140713.$ We have
$d_{0}=d_{1}=d_{2}=0,d_{3}=1$. We obtain the following terms
$d_{4}=18,d_{5}=334,d_{6}=6205,$ $d_{7}=115267,d_{8}=2141252.$ Since
$d_{8}\leq 9140713\leq d_{9},$ we obtain the
key$~{}\left(18,10,13,3,6\right)$.
We have
$n_{1}=\mathbf{4}\cdot d_{8}+575705,$
$575705=\mathbf{4}\cdot d_{7}+114637,\newline 114637=\mathbf{18}\cdot
d_{6}+2947,\newline 2947=\mathbf{8}\cdot d_{5}+275,\newline
275=\mathbf{15}\cdot d_{4}+5,\newline 5=\mathbf{5}\cdot d_{3},$ therefore, we
obtain the labels for the encrypted text $\left(4,4,18,8,15,5\right)$. It
results that the first block is encrypted in the word EESIPF. To do this
faster, we can use the below MAPLE procedures:

    > crypt(JOHN);
     9140713
    >lineq(4, 18, 10, 13, 3, 9140713);
    6
    40418081505
    > decrypt(40418081505);
    EESIPF

For the second block, xHAS, we get the number $n_{2}=26070018.$ Using the same
algorithm, we obtain the following labels for the encrypted text
$\left(12,3,4,13,1,13\right).$ The second block is encrypted in the word
MDENBN. To ease calculations, we can use the below MAPLE procedures:

    > crypt(xHAS);
     26070018
    >lineq(4, 18, 10, 13, 3, 26070018);
    6
    120304130113
    > decrypt(120304130113);
    MDENBN

To the third block, xAxD, corresponds the number $n_{3}=26002603.$ We obtain
the labels for the encrypted text $\left(12,2,12,7,13,13\right).$ It results
that the third block is encrypted in the word MCMHNN, by using the below MAPLE
procedures:

    > crypt(xAxD);
     26002603
    >lineq(4, 18, 10, 13, 3, 26002603);
    6
    120212071313
    > decrypt(120212071313);
    MCMHNN

To the last block, OGxx, corresponds the number $n_{4}=14062626$. We get
the labels for the encrypted text $\left(6,10,10,1,3,6\right).$ The fourth
block is encrypted in the word GKKBDG, as we can see by using the below MAPLE
procedures:

    > crypt(OGxx);
     14062626
    >lineq(4, 18, 10, 13, 3, 14062626);
    6
    61010010306
    > decrypt(61010010306);
    GKKBDG

Therefore, the encrypted text is EESIPFMDENBNMCMHNNGKKBDG.

For decoding, we use the key $\left(18,10,13,3,6\right).~{}$We split the
message in blocks of length $6$: EXKLFC, MDENBN, MCMHNN, GKKBDG. For this, we
can use the below MAPLE procedure:

    > split(EESIPFMDENBNMCMHNNGKKBDG, 6);
    [EESIPF MDENBN MCMHNN GKKBDG]

We compute $d_{4}=18,d_{5}=334,d_{6}=6208,$ $d_{7}=115267,d_{8}=2141252$. For
EESIPF, we obtain $n_{1}=2141252\cdot 4+115267\cdot 4+6205\cdot 18+334\cdot
8+18\cdot 15+5\cdot 1=9140713,$ that means we get the word JOHN, as we can see
below:

    > delineq(EESIPF,18, 10, 13, 3);
    9140713

2) If we consider $k=4,$ $a_{1}=18,a_{2}=4,a_{3}=13,a_{4}=3$, then, the
encoding key is SEND, or 18041303. In this situation, we have the following
terms $d_{4}=18,d_{5}=328,d_{6}=5989,$ $d_{7}=109351,d_{8}=1996592,$ and the
encrypted message is E28OQAPNBAPCPNAHPDKHAOIAM.

3) If the alphabet has 27 letters and we received X032121MBD, as an encrypted
block, and we know that this block has ”length” $6,$ we obtain the labels
$\left(23,32,121,12,1,3\right),$ that means $23032121012001003.$ Indeed, $23$
is label for X, the next two labels are $032$ and $121,$ since we know that
the labels have the same number of digits, in our case $3,$ since the block
has level $6$. For this situation, the procedure decrypt was modified in
decryptM by inserting new characters y for $032$ and z for $121.$

    > decryptM(23032121012001003);
    XyzMBD

If the above block has length $7,$ then the labels are
$\left(23,3,21,21,12,1,3\right)$.

3.2. An application in Coding Theory

In [1], [4], [9], were presented some applications of Fibonacci numbers and
Fibonacci $p-$numbers in Coding Theory. With these elements, were defined some
special matrices used for sending messages. In the following, we generalized
these results, using the matrices given in relations $\left(2.11\right)\,\
$and $\left(2.12\right)$ as a coding matrices for sending messages. These
matrices are associated to difference equation of order $k,$ $k$ a natural
number, $k\geq 2,$ given by the relation $\left(2.9\right),$ with
$a_{1},...a_{k}$ positive integers$.~{}$We attach to this equation the
following equation of degree $k$

$x^{k}-a_{1}x^{k-1}-a_{2}x^{k-2}-...-a_{k}=0$ (3.2.)

and we assume that equation $\left(3.2\right)$ has $k$ distinct real roots
$\alpha_{1},\alpha_{2},...,\alpha_{k},$ such that $\alpha_{1}>1$ and
$\alpha_{1}>\left|\alpha_{i}\right|,i\in\\{2,...,k\\}.$

We assume that the initial message is represented as a square matrix
$M=\left(m_{ij}\right)_{i,j\in\\{1,2,...,k\\}}~{}$of order $k,$ with
$m_{ij}\geq 0,i,j\in\\{1,2,...,k\\}.$ Considering that
$D_{k}^{n}=\left(h_{ij}\right)_{i,j\in\\{1,2,...,k\\}}$ is the coding matrix
and $D_{k}^{-n}$ is the decoding matrix, we will use the relation

$D_{k}^{n}\cdot M=E,$ (3.3.)

as encoding transformation and the relation

$M=D_{k}^{-n}\cdot E,$ (3.4.)

as decoding transformation. We define $E$ as code-message matrix.

Theorem 3.4. Withe the above notations, the following statement is true

$\frac{e_{ij}}{e_{(i+r)j}}\thickapprox\alpha_{1}^{r},r\in\\{1,2,...,k-1\\},$
(3.5.)

for $i+r\leq k.\vskip 6.0pt plus 2.0pt minus 2.0pt$

Proof. Since $a_{1},\alpha_{2},...,\alpha_{k}$ are roots of the equation
$\left(3.2\right),$ we have

$d_{n}=A_{1}\alpha_{1}^{n}+A_{2}\alpha_{2}^{n}+...+A_{k}\alpha_{k}^{n},$
(3.6.)

with coefficients $A_{i}\in\mathbb{R},i\in\\{1,2,...,k\\},$ obtained from the
initial conditions of difference equation $\left(2.9\right).$ Using relation
$\left(3.3\right),$ we get
$e_{ij}=\overset{k}{\underset{t=1}{\sum}}h_{it}m_{tj}$ . From relation
$\left(2.12\right),$ it results that the elements $h_{ij}$ are non-zero
elements and are represented as a linear combination of the terms
$d_{n+k-i},d_{n+k-i-1},...,d_{n-i+1},$

$h_{ij}=\overset{k}{\underset{t=1}{\sum}}B_{tj}d_{n+k-i-t+1},$

$B_{tj}\in\\{0,1,a_{1},a_{2},...a_{k}\\}.$ For example,

$h_{i1}=d_{n+k-i}=\overset{k}{\underset{t=1}{\sum}}B_{t1}d_{n+k-i-t+1},$

where $B_{11}=1,B_{t1}=0,t\in\\{2,...k\\},$

$h_{i2}=\overset{k-1}{\underset{t=1}{\sum}}a_{t+1}d_{n+k-t-i}=\overset{k}{\underset{t=1}{\sum}}B_{t2}d_{n+k-i-t+1},$

where $B_{12}=0,B_{22}=a_{2},B_{32}=a_{3},...B_{k2}=a_{k},~{}$etc.

Therefore, we obtain

$e_{ij}=\overset{k}{\underset{t=1}{\sum}}\left(\overset{k}{\underset{q=1}{\sum}}B_{qt}d_{n+k-i-q+1}\right)m_{tj}$

and

$e_{(i+r)j}=\overset{k}{\underset{t=1}{\sum}}\left(\overset{k}{\underset{q=1}{\sum}}B_{qt}d_{n+k-i-
r-q+1}\right)m_{tj}.$

Using relation $\left(3.7\right),$ it results

$\underset{i\rightarrow\infty}{\lim}\frac{e_{ij}}{e_{(i+r)j}}=\underset{i\rightarrow\infty}{\lim}\frac{\overset{k}{\underset{t=1}{\sum}}\left(\overset{k}{\underset{q=1}{\sum}}B_{qt}d_{n+k-i-q+1}\right)m_{tj}}{\overset{k}{\underset{t=1}{\sum}}\left(\overset{k}{\underset{q=1}{\sum}}B_{qt}d_{n+k-i-
r-q+1}\right)m_{tj}}=$
$=\underset{i\rightarrow\infty}{\lim}\frac{\overset{k}{\underset{t=1}{\sum}}\left(\overset{k}{\underset{q=1}{\sum}}B_{qt}\left(\overset{k}{\underset{s=1}{\sum}}A_{s}\alpha_{s}^{n+k-i-q+1}\right)\right)m_{tj}}{\overset{k}{\underset{t=1}{\sum}}\left(\overset{k}{\underset{q=1}{\sum}}B_{qt}\left(\overset{k}{\underset{s=1}{\sum}}A_{s}\alpha_{s}^{n+k-i-
r-q+1}\right)\right)m_{tj}}=$
$=\underset{i\rightarrow\infty}{\lim}\frac{\overset{k}{\underset{q,s=1}{\sum}}\Gamma_{sq}\alpha_{s}^{n+k-i-q+1}}{\overset{k}{\underset{q,s=1}{\sum}}\Gamma_{sq}\alpha_{s}^{n+k-i-
r-q+1}}=$
$\text{=}\underset{i\rightarrow\infty}{\lim}\frac{\Gamma_{11}\alpha_{1}^{n+k-i}\text{+}\Gamma_{12}\alpha_{1}^{n+k-i-1}\text{+...+}\Gamma_{1k}\alpha_{1}^{n-i+1}\text{+}\Gamma_{21}\alpha_{2}^{n+k-i}\text{+}\Gamma_{22}\alpha_{2}^{n+k-i-1}\text{+...}}{\Gamma_{11}\alpha_{1}^{n\text{+}k\text{-}i\text{-}r}\text{+}\Gamma_{12}\alpha_{1}^{n\text{+}k\text{-}i\text{-}r\text{-}1}\text{+...+}\Gamma_{1k}\alpha_{1}^{n\text{-}i\text{-}r\text{+}1}\text{+}\Gamma_{21}\alpha_{2}^{n\text{+}k\text{-}i\text{-}r}\text{+}\Gamma_{22}\alpha_{2}^{n\text{+}k\text{-}i\text{-}r\text{-}1}\text{+...}}\text{=}L.$

where
$\Gamma_{sq}=A_{s}\left(\overset{k}{\underset{t=1}{\sum}}B_{qt}m_{tj}\right).$
We remark that $\Gamma_{sq}\neq 0,~{}$due to the choice of the elements
$A_{s},B_{qt}$ and $m_{tj}$. Since $\alpha_{1}>1$ and
$\alpha_{1}>\left|\alpha_{i}\right|,i\in\\{2,...,k\\},$ we obtain that
$L=\alpha_{1}^{r}.~{}$From here, if we consider $i+r\leq k,$ it results
$\frac{e_{ij}}{e_{(i+r)j}}\thickapprox\alpha_{1}^{r},r\in\\{1,2,...,k-1\\}.\vskip
6.0pt plus 2.0pt minus 2.0pt$

Remark 3.5. From relation $\left(3.3\right),~{}$we have $\det M\cdot\det
D_{k}^{n}=\det E,$ then

$\left(-1\right)^{\left(k+1\right)n}a_{k}^{n}\det M=\det E,$ (3.7.)

by using $\left(2.13\right).$

The matrices $D_{k}^{n}$ are used for sending messages. In this way, this
method ensures infinite variants for the chioce of transformation of the
initial matrix $M$.

The method developed above has the property to detect and correct errors in
the submitted message $E$. For this purpose, we used Theorem 3.4 and relation
$\left(3.7\right)$. From condition $\left(3.7\right),$ we have the relation
between $\det M$ and $\det E$. The determinant of the matrix $M,$ $\det M,$ is
used as a checking element for the code-message matrix $E,$ when we received
it using a communication channel. After the matrix $E$ and $\det M~{}$were
received, we compute $\det E~{}$and we check if the relation
$\left(3.7\right)$ is satisfied. If the answer is positive, it results that
the matrix $E$ and $\det M~{}$were transmitted without errors. If the answer
is negative, we have that the matrix $E$ or $\det M$ were received with
errors. For correct the errors, we will use relation $\left(3.5\right)$.
First, we suppose that we have single error in the received matrix $E$. Since
such an error can appears on the position
$\left(i,j\right),i,j\in\\{1,2,...,k\\}$, it results that we can have $k^{2}$
possibilities. If we check for all $k^{2}$ possibilities and we don’t obtain
natural numbers as solution for the received error, it is possible to have
double errors, triple errors,…, $k^{2}-1$ errors. If in all these situations
we do not obtain positive integer solutions, it results that $\det M$ was sent
with errors or the matrix $E$ has $k^{2}-$fold errors and this matrix is not
correctable. Therefore, the matrix $E$ must be rejected.

This method, in which we use the matrices $D_{k}^{n}$ as encoding matrices and
$D_{k}^{-n}$ as decoding matrices, generalizes the methods developed in [9]
for Fibonacci numbers, in [1] and [2] for Fibonacci $p-$numbers and in [11]
for Fibonacci $k-$numbers to a difference equation of order $k,$
$k\in\mathbb{N}$, $k\geq 2,$ given by the relation $\left(2.9\right),$ with
$a_{1},...a_{k}$ positive integers.

The error correcting codes are used widely in modern communications networks.
There are many parameters associated to error correcting codes which
determines the ability of a code to detect and correct errors (Hamming
distance, the rate of a code, etc.) One of these parameters is the potential
error correction coefficient $S$, which is the ratio between all correctable
errors and all detectable errors.

We can remark that the code message matrix $E$ can contain single, double,….,
$k^{2}-$fold errors. Therefore, we have

$\complement_{k^{2}}^{1}+\complement_{k^{2}}^{2}+...+\complement_{k^{2}}^{k^{2}}=2^{k^{2}}-1$

possible errors. Since $k^{2}-$fold errors from the code-message matrix $E$
are not correctable, we can correct $2^{k^{2}}-2$ errors. Therefore we get

$S=\frac{2^{k^{2}}-2}{2^{k^{2}}-1}\thickapprox 1,$

that means the correctable possibility of the method is about $100\%.$

This remark generalized the results obtained in [1] and [2] for particular
case of Fibonacci $p-$numbers to a difference equation of degree $k\,$,
defined by $\left(2.9\right)$.

Conclusions. In this paper, we presented some applications of a difference
equation of degree $k$ in Cryptography and Coding Theory.

The algorithm for encryption/decryption messages has some advantages:

i) a block of length $s$ is transformed into a block of different length $r;$

ii) each natural number $n$ has a unique representation using the terms of
this sequence, that means a number cannot have the same encrypted value as
another number;

iii) this method give us a high versatility, ensuring infinite variants for
the chioce of the encrypted keys.

Moreover, using the matrix associated to a difference equation of degree $k$
we have a lot of possibilities to chose a matrix for sending messages,
defining error correcting codes with a very good potential error correction
coefficient.

The above results encourages us to study these equations for finding other
interesting applications of them.

Appendix

In the following, we present some MAPLE procedures used in the encryption and
decryption processes.

    crypt:=proc(st)
    local ll,nn,ss,ii,num,n;
    num := table([’A’=0, ’B’=1, ’C’=2, ’D’=3, ’E’=4, ’F’=5,
    ’G’=6, ’H’=7,’I’=8, ’J’=9, K’=10, ’L’=11, ’M’=12, ’N’=13,
    ’O’=14,’P’=15, ’Q’=16, ’R’=17, ’S’=18, ’T’=19, ’U’=20,
    ’V’=21, ’W’=22, ’X’=23, ’Y’=24, ’Z’=25, ’x’=26]):
    ll := length(st): nn := 1: for ii from 1 to ll do
    ss := num[substring(st, ii .. ii)]:
    nn := 100* nn+ss: od:
    n:=nn-10^(2* ll):print(n):end:
    save crypt, ‘crypt.m‘;
    crypt(JOHN)
    9140713

    lineq := proc (k, a1, a2, a3, a4, n)
    local AA, i, jj, cc, nn, NN, s,j;
    j:=100; AA := array(0 .. j); AA[0] := 0; AA[1] := 0; AA[2] := 0; AA[3] := 1;
    for i from 0 to j-4 do AA[i+4] := AA[i+3]*a1+AA[i+2]*a2+AA[i+1]*a3+AA[i]*a4;
    if n < AA[i+4] and AA[i+3] < n then jj := i+3 end if;
    s := jj-k+2 end do; cc := array(1 .. s); cc[1] := trunc(n/AA[jj]);
    nn := n-cc[1]*AA[jj]; for i from 2 to s do cc[i] := trunc(nn/AA[jj-i+1]);
    nn := nn-cc[i]*AA[jj-i+1] end do; NN := 10^(2*s-2)*cc[1];
    for i from 2 to s do NN := NN+10^(2*s-2*i)*cc[i] end do;print(s); print(NN);
    end proc; save lineq, ‘lineq.m‘; lineq(4, 18, 10, 13, 3, 9140713):
    6
    40418081505

    decrypt:=proc(nn)
    local alpha, ss,mm,rr,ii, ans,A,B,C,D,E,F,G,H,II,J,K,L,M,
    O,P,Q,R,S,T,U,V,W,X,Y,Z,x;
    alpha:= table([0 = A, 1 = B, 2 = C, 3 = D, 4 = E,
    5 = F, 6 = G, 7 = H, 8 = I, 9 = J, 10 = K, 11 = L,
    12 = M, 13 = N, 14 = O,15=P,16=Q,17=R,18=S,19=T,
    20=U,21=V,22=W,23=X,24=Y,25=Z,26=x]):
    mm := nn: rr:=floor(trunc(evalf(log10(mm)))/2)+1: ans:=‘ ‘:
    for ii from 1 to rr do mm:=mm/100: ss:=alpha[frac(mm)*100]:
    ans:=cat(ss,ans): mm:=trunc(mm) od: ans;
    end: save decrypt, ‘decrypt.m‘;
    decrypt(40418081505);
    EESIPF

    decryptM:=proc(nn)
    local alpha, ss,mm,rr,ii, ans,A,B,C,D,E,F,G,H,II,J,K,L,M,
    O,P,Q,R,S,T,U,V,W,X,Y,Z,x;
    alpha:= table([0 = A, 1 = B, 2 = C, 3 = D, 4 = E,
    5 = F, 6 = G, 7 = H, 8 = I, 9 = J, 10 = K, 11 = L,
    12 = M, 13 = N, 14 = O,15=P,16=Q,17=R,18=S,19=T,
    20=U,21=V,22=W,23=X,24=Y,25=Z,26=x,32=y,121=z]):
    mm := nn: rr:=floor(trunc(evalf(log10(mm)))/3)+1: ans:=‘ ‘:
    for ii from 1 to rr do mm:=mm/1000: ss:=alpha[frac(mm)*1000]:
    ans:=cat(ss,ans): mm:=trunc(mm) od: ans;
    end: save decryptM, ‘decryptM.m‘;
    decryptM(23032121012001003);
    XyzMBD

    delineq := proc (st, a1, a2, a3, a4)
    local AA, BB, cc, ll,nn, num, s,ii;
    num := table([’A’=0, ’B’=1, ’C’=2, ’D’=3, ’E’=4, ’F’=5,
    ’G’=6, ’H’=7,’I’=8, ’J’=9, K’=10, ’L’=11, ’M’=12, ’N’=13,
    ’O’=14,’P’=15, ’Q’=16, ’R’=17, ’S’=18, ’T’=19, ’U’=20,
    ’V’=21, ’W’=22, ’X’=23, ’Y’=24, ’Z’=25, ’x’=26]):
    ll := length(st):  BB:=array(1..ll):AA:=array(0..ll+2):
    for ii from 1 to ll do BB[ii]:= num[substring(st, ii .. ii)]  od:
    AA[0] := 0: AA[1] := 0: AA[2] := 0: AA[3] := 1:
    for ii from 0 to ll-2  do
    AA[ii+4] := AA[ii+3]*a1+AA[ii+2]*a2+AA[ii+1]*a3+AA[ii]*a4;od:
    nn[]:=0: for ii from 1 to ll do nn:=nn+BB[ii]*AA[ll+3- ii]: od:
    print(nn);end proc; save delineq, ‘delineq.m‘;
    delineq(EESIPF,18, 10, 13, 3):
    9140713

    split := proc (st, k)
    local ll, ii, rr, AA;
    ll := length(st); rr := 0; AA := array(1 .. ll/k);
    for ii to ll/k do AA[ii] := substring(st, ii+rr .. ii+k-1+rr);
    rr := rr+k-1 end do; print(AA);
    end proc; save split, ‘split.m‘;
    split(EESIPFMDENBNMCMHNNGKKBDG, 6);
    [EESIPF MDENBN MCMHNN GKKBDG]

Acknowledgments. The author thanks the referees for their suggestions and
remarks which helped me to improve this paper.

References

[1] Basu, M., Prasad, B., The generalized relations among the code elements
for Fibonacci coding theory, Chaos, Solitons and Fractals, 41(2009),
2517-2525.

[2] Basu, M., Prasad, B., Coding theory on the m-extension of the Fibonacci
p-numbers, Chaos, Solitons and Fractals, 42(2009), 2522–2530.

[3] Brown, Jr.J. L., Note on Complete Sequences of Integers, The American
Mathematical Monthly, (68)(6)(1961), 557-560.

[4] Esmaeili, M., Esmaeili, Mo., A Fibonacci-polynomial based coding method
with error detection and correction, Computers and Mathematics with
Applications, 60(2010), 2738-2752.

[5] Sergio Falcon, S., Plaza, A., On k-Fibonacci numbers of arithmetic
indexes, Applied Mathematics and Computation, 208(2009), 180–185.

[6] Hoggatt, Jr.V. E., King, C., Problem E1424, Amer. Math. Monthly, Vol.
67(1960), p. 593.

[7] Johnson, R.C., Fibonacci numbers and matrices,
http://maths.dur.ac.uk/~dma0rcj/PED/fib.pdf

[8] Klein, S.T., Ben-Nissan, M.K., On the Usefulness of Fibonacci Compression
Codes, The Computer Journal, 2005, doi:10.1093/comjnl/bxh000

[9] Stakhov, A.P., Fibonacci matrices, a generalization of the “Cassini
formula”, and a new coding theory, Chaos, Solitons and Fractals, 30(2006),
56-66.

[10] Stakhov, A.P., The “golden” matrices and a new kind of cryptography,
Chaos, Solitons and Fractals, 32(2007), 1138–1146.

[11] Tahghighi, M., Jafaar, A., Mahmod, R., Said, M. R., Generalization of
Golden Cryptography based on $k$-Fibonacci Numbers, International Conference
on Intelligent Network and Computing (ICINC 2010), 2010.

[12] Yeates, A. C., Application of Linear Sequences to Cryptography, Honors
Theses, Paper 191, 2013.

[13] Zeckendorf, E., Reprćsentation des nombres naturels par une somme de
nombres de Fibonacci ou de nombres de Lucas, Bull. Soc. R. Sci. Lićge,
41(1972), 79-182.