{'fraction_non_alphanumeric': 0.04700696055684455, 'fraction_numerical': 0.008213457076566125, 'mean_word_length': 4.600571725571726, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Inhomogeneous Dirac-Bloch Equations for Graphene Interacting with Structured Light

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 [1] 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 [2], universal absorbance [3], anomalous Quantum Hall Effect [4,5], and high nonlinear optical response [6,7], 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) [8,9],hexagonal boron nitride (hBN) [10,11], and black phosphorous [12]. 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- [13] and third-harmonic generation [14], four-wave mixing [15], and optical limiting [16]. While graphene is intrinsically centrosymmetric and therefore lacks second order response, recent works have demonstrated how it is possible to induce secondharmonic generation in graphene by inducing stress or strains [17] or by exploiting the fact that an ultra-intense pulse can open its own gap near Dirac points, thus breaking centrosymmetry [18].
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 [19], and represent a versatile tool used in several different experimental situations, and show wide range of application in different areas, namely, optical tweezers [20], quantum information [21], microscopy [22], spectroscopy [23], classical [24] and quantum [25] 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 [26], photoexcitation of graphene [27] and structured nonlinear optics in 2D materials have been stud-ied. 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 [28] 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 where φ c (r) is the Coulomb potential, V B (r) is the crystalline (periodic) potential,ψ(r) is the electron wavefunction operator, and A(r, t) is the vector potential. From here henceforth, we are assuming to work in the Coulomb gauge, where ∇·A = 0, and we also neglect the Coulomb interaction terms, since they don't contribute significantly for the case of graphene [29]. 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 [30] to the above Hamiltonian, to get the usual light-matter interction in the electric dipole approximation as follows: where E T (r, t) = −∂ t A(r, t) is the (transverse) electric field. To calculate the graphene Hamiltonian, we use the following Bloch representation for the electron wavefunction operator where φ A,B (r) are 2p z -orbital wavefunctions at site A and B (and they form an orthogonal normalised basis set), R n = ja 1 + ma 2 is a Bravais lattice vector (see Fig.  1), andâ k ,b k are the electron site operators, annihilating an electron on site A or B, respectively. By substituting Eq. (4) into the first of Eqs. (3), it is not difficult to see, that the free HamiltonianĤ 0 can be then reduced tô where γ is the hopping amplitude and f (k) = 1 + exp (ik · a 1 ) + exp (ik · 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(k) = ±γ|f (k)|, and its corresponding eigenvectors where ϕ(k) = arg [f (k)]. 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Ĥ 0 ,obtaining Notice that this is equivalent to the electron-hole representation discussed, for example, in Ref. [31], wherê e k annihilates an electron with momentum k, andĥ † −k creates a hole with momentum −k.
The free Hamiltonian in the electron-hole representation can be then written as followŝ

A. Low-Energy Approximation
We can expand the eigenvalues E(k) around the zeros of f (k), which happen to be the two inequivalent Dirac points K and K' [see Fig. 1 [32] to obtain the lowenergy approximation ofĤ 0 around the Dirac points. To do that, we first redefine the coordinates of the Dirac points as K = (2π/3a Substituting this result in the expression of the Hamilto-nianĤ 0 and defining the Fermi velocity as v f 3aγ/2, we get the familiar form of the low-energy graphene Hamiltonian, i.e.,Ĥ where ξ = ±1 is the valley index, and σ i are Pauli matrices. The eigenvalues ofĤ(k) are now simply E(k) = v f |k|. Notice, moreover, that in this approximation And, that the eigenvectors u ± (k) acquire a valley depen-

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 x = r − R n , assuming only nearest neighbour coupling, such that R = R m − R n only allows for = {0, ±1}, and we make the extra assumption that the electric field is constant over a single unit cell, i.e., Notice that this approximation is consistent with previous work dealing with semiconductor Bloch equations in the presence of inhomogeneous fields [33] Since we are doing these calculation in the site basis, we will have to deal with both diagonal terms, proportional toâ † kâ k andb † kb k , and nondiagonal terms, proportional toâ † kb 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 bŷ where the integral in the second line is zero because φ A (x) and φ B (x) are orthogonal, 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 E T (R n , t) = 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 and then introduce the generalised Rabi frequency as follows so that the interaction Hamiltonian above can be written asĤ Notice, that in order for Ω(k, 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 V s/m, which is consistent with the fact that |E| = |∂ t A| τ A 0 , with τ being the characteristic timescale of the impinging electric field (for an impinging pulse, τ is the duration of the pulse).
We are now in the position to write down the total interaction HamiltonianĤ I in the site-basis as followŝ (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â † kb 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 where |Ω(k, q, t)| = Ω 2 x (k, q, t) + Ω 2 y (k, q, t) is the modulus of the Rabi frequency, and is its phase.

A. 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Ĥ(k), as detailed in Ref. [34] as where the conduction and valence band eigenstates are given as |ψ c,v (k) = e ik·r |u ± (k) , where the + (−) sign refers to conduction (valence) band. Substitution gives Notice that in the vicinity of a Dirac point ϕ(k) π 2 + ξ arctan(k y /k x ), and the dipole moment becomes or, by introducing the helicity basisĥ ξ = 1

B. 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., Substituting this into Eq. (18), and using the anticommutation relations {ĥ † k ,ĥ q } = δ k,q we get where with ∆ϕ(k, q) = ϕ(k) + ϕ(k − q)], and ω(k, 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 † e and h † h terms) but also interband ones (e † h † and he terms). Moreover, written in this form it is equivalent to the linear minimal coupling Hamiltonian in Ref. [36].

IV. DIRAC-BLOCH EQUATIONS
To calculate the Dirac-Bloch equations, we first introduce the following two-body populations and polarisation and calculate their time derivative with the help of the Heisenberg equation. Using the anti-commutation rules {ê † k ,ê k } = δ k,k and {ĥ † −k ,ĥ −k } = δ k,k we get, after some long but straightforward algebra, the following twobody Dirac-Bloch equations where E(k) [E(k )] are the energy eigenvalues of valence and conduction bands in state k (k ), and we have defined 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 (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., where Ω = Ω x (k, q, t)x + Ω y (k, q, t)ŷ, σ = σ xx + σ yŷ , ψ k = (ê kĥ † −k ) T , and U (k) being the matrix whose columns are the eigenstates ofĤ 0 calculated before. If we then define the following 2 × 2 matrices and call m x,y ij (k, k , q) their matrix elements, whose explicit expressions are given as follows , the current operatorĵ =ĵ xx +ĵ yŷ in Eq. (31) can be then written as followsĵ Taking the expectation values we finally get where we have used the fact that m x,y 11 = −m x,y 22 and m x,y 12 = im x,y 12 = −i(m x,y 21 ) * . Notice, that the current above has the same formal expression that the one for the homogeneous case, once the replacements n e,h (k, k ) → n e,h (k), p(k, k ) → p(k), and m x,y ij (k , k, q, t) → m x,y ij (k, k, 0, t) are made.

VI. NONLINEAR OPTICAL RESPONSE
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., where τ is the pulse duration, t 0 an arbitrary offset, R = {R, φ}, and with w(z) = w 0 1 + (z/z R ) 2 is the Gaussian beam waist (with z R being the Rayleigh range of the LG-beam), C ,p = w 0 2p!/[π(p + | |)!], and L | | p (x) are the associated Laguerre polynomials [37], and represents the OAM index, while p represents the radial index of the LG-beam. We then solve Eqs. (29) for different values of the OAM parameter (and set p = 0, for simplicity), and then use Eqs. (36) to compute the correspondent nonlinear current j(t) by integrating over the transverse shape of the pulse, and we then calculate the nonlinear signal as I(ω) ∼ |ω J(ω)| 2 , where J(ω) is the Fourier transform of the nonlinear current.
The resulting nonlinear spectrum for different values of the OAM parameter 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 [32].
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 [18] 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.

A. 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 4dimensional reciprocal space {k, k }, 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 spatiallydependent nonlinear current where J x,y (k, k , 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 J x,y (k, k , 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 k 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 2dimensional one by setting k = k in Eqs. (29). This, essentially, means, that we restrict the domain of definition of k 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 k = k essentially erases any information about the even harmonics. This, again, is not surprising, since the limit k = k 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., j y (R, t) = k,q J y (k, q, t)e ik·R .
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 spatiallyvarying 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.