Theoretical calculations of nonlinear optical calculations of 2D materials

One important feature of two dimensional (2D) materials is that they possess an exceptional nonlinear optical (NLO) response to light, with conductivities that are several orders of magnitude larger than their 3D counterparts. The theoretical descriptions of these NLO responses in crystalline systems involve two different representations of the perturbation: the length and velocity gauges. The former has been the formalism of choice for the past two decades; the latter was implemented only recently, due to concerns that it could not be pratically implemented without breaking sum rules – a set of identities that ensure the equivalence between the two formalisms – which would then render the results unphysical. In this work, we shall review and summarize our contributions to the study of the two formalisms and of their relationship by means of the aforementioned sum rules. Nonlinear optical responses The description of the nonlinear optical response of electrons in a crystal to an external, spatially constant and monochromatic, electric field was first formulated – in terms of the density matrix formalism – in Claudio Aversa and J. E. Sipe’s work of 1995 [1]. There, it is considered that, in the absence of the electric field perturbation, an electron in a crystal can described by the periodic Hamiltonian, H0(k), of eigenvalues (bands), ks, and eigenvectors (Bloch states), Ψks, both of which are labelled by a crystal momentum k, a vector in the first Brillouin zone (FBZ) and a band index, s. Electrons are deemed to be non-interacting among themselves – that is, they are subjected only to the electric field of light that is shining on the material – and so the many-particle ground state is the Fermi sea, a state whose occupation numbers follow from the Fermi-Dirac distribution function, 〈ckscks′ 〉 = δss′ fks. In rather simplistic terms, the effect that the external electric field has on this system of electrons in a set of bands is to change the populations within each band and to introduce transitions between them so that the matrix elements of the density matrix (DM), ρkss′ (t) = 〈cks(t)cks′ (t)〉, are now dynamical. The equations of motion for the DM matrix elements are given by, i ρ̇kss′ (t) = ( ks − ks′ ) ρkss′ (t) + [V(t), ρ(t)]kss′ . (1) A perturbative treatment of ρkss′ (t) then allows us to study the crystal’s response functions – here we shall consider the conductivity – in its different orders of the electric field [2]. To ∗e-mail: gbventura@fc.up.pt © The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0 (http://creativecommons.org/licenses/by/4.0/). EPJ Web of Conferences 233, 03001 (2020) https://doi.org/10.1051/epjconf/202023303001 CMPNC 2019


Nonlinear optical responses
The description of the nonlinear optical response of electrons in a crystal to an external, spatially constant and monochromatic, electric field was first formulated -in terms of the density matrix formalism -in Claudio Aversa and J. E. Sipe's work of 1995 [1]. There, it is considered that, in the absence of the electric field perturbation, an electron in a crystal can described by the periodic Hamiltonian, H 0 (k), of eigenvalues (bands), ks , and eigenvectors (Bloch states), Ψ ks , both of which are labelled by a crystal momentum k, a vector in the first Brillouin zone (FBZ) and a band index, s. Electrons are deemed to be non-interacting among themselves -that is, they are subjected only to the electric field of light that is shining on the material -and so the many-particle ground state is the Fermi sea, a state whose occupation numbers follow from the Fermi-Dirac distribution function, c † ks c ks = δ ss f ks . In rather simplistic terms, the effect that the external electric field has on this system of electrons in a set of bands is to change the populations within each band and to introduce transitions between them so that the matrix elements of the density matrix (DM), ρ kss (t) = c † ks (t)c ks (t) , are now dynamical. The equations of motion for the DM matrix elements are given by, A perturbative treatment of ρ kss (t) then allows us to study the crystal's response functionshere we shall consider the conductivity -in its different orders of the electric field [2]. To do so, one must first specify the form of the perturbation, V(t). Due to the inherent gauge freedom involved in expressing the electric field in terms of the scalar and vector potentials, there is no unique way of writing the perturbation; for simplicity, either the vector potential or the scalar potential is set to zero, and obtains the length and velocity gauges, respectively: Similarly to electromagnetism, the choice of a gauge -we shall presently see -has its own particular advantages and disadvantages. We also note that the equivalence between the descriptions associated to each gauge, though following directly from gauge invariance, has to be reconsidered in light of practical calculations.

Length gauge
In the length gauge, the perturbation is expressed in terms of the position operator and, as such, it is not invariant under lattice translations. Some attention is thus needed when expressing this object in terms of Bloch states [3], for ξ kss , the Berry connection [1]. Two equivalent formulations using this gauge have been put forth. In [1], the authors separate the perturbation into intraband (i : s = s) and interband (e : s s) portions so that the response at a given order is obtained by considering all relevant possible pairings of the two parts of the perturbation, e.g., in second order, one of the terms reads as ρ ie . In our first work on this subject [4], no such distinction is made: the position operator in the Bloch representation is treated as a covariant derivative, D kss = δ ss ∇ k −iξ kss . This formulation makes it clear that the Bloch representation of the position matrix element is gauge covariant under phase shifts of the Bloch states. Following [4], the equation of motion for the DM reads as, By expanding ρ kss (t) in powers of the external field, one can see the two main features of this gauge: the complexity of the solutions increases for high order responses, as an nth-order solution involves the k-space derivative of the (n − 1)th-order one, which is (in principle) a rather complicated function in that variable; in the frequency domain, one can see that energy dependences, ∆ kss = ks − ks , appear only in the denominators. Practical calculations in this gauge can then use a truncation of the Hilbert space, both in the number of bands selected and the region of k-space that is taken into consideration, which makes length gauge the most suitable gauge for analytical work.

Velocity gauge
In the case of the velocity gauge, the interaction with the external electric field is introduced via minimum coupling. This means that the form of the associated perturbation, V(t), depends on the Hamiltonian describing the unperturbed system, Eq.(3). So, while the perturbation in the length gauge remains unchanged -whether one is considering the continuum Hamiltonian or a tight-binding one -the same is not true in the velocity gauge, where the perturbation has to be considered on a case-by-case basis. To illustrate a second point concerning the velocity gauge, we consider the continuum Hamiltonian, H 0 (k) = ( 2 /2m)(−i∇ + k) 2 + U(r), for U(r) the periodic lattice potential. The corresponding perturbation, V A (t), written in the Bloch basis, reads as, The relevant portion of the perturbation is proportional to the velocity operator, which being lattice periodical, can be easily expressed in terms of Bloch states, v kss = δ ss ∇ k ks − i∆ kss ξ kss .
One can easily see that for such a perturbation, the solutions to the equations of motion for the DM will have an energy dependence on ∆ kss in the numerators (as well as in the denominators) since the interband matrix elements of the velocity carry such energy factors. The velocity gauge, therefore, does not allow for a truncation of the Hilbert space along the same lines as the length gauge and so practical calculations in the former gauge have to be implemented differently than those in the latter. This sticking point was only fully solved in our second work [5] and required a closer study on how equivalence is ensured in these calculations.

Equivalence between procedures and sum rules
Gauge invariance tells us that an observable computed in the length gauge must be necessarily equal to its counterpart in the velocity gauge, In the present case, the operators are related by a unitary transformation that describes timedependent translations in momentum space. By expressing these operators in the Bloch basis, one can see that Eq.(8) holds if and only if a set of sum rules are fulfilled; these sum rules are a set of integrals that must be shown to be trivial -a condition that is satisfied when the original operators are defined in the full FBZ and if the set of bands considered remains untruncated [4]. By enforcing these criteria from the start, it is possible to formulate a procedure for practical calculations in the velocity gauge whose results are equivalent to those in the length gauge [5]. We must emphasize that this seemingly technical point has important consequences for the calculations of NLO since the velocity gauge -being fully diagonal in k-space -is better suited for numerical implementations.

Numerical results in the velocity gauge: THG in the monolayer of graphene
A proper calculation in the velocity gauge requires the use of an Hamiltonian that is defined in the full FBZ and whose set of bands cannot be truncated. In [5], we compared the results for the third harmonic generation [2] using the velocity gauge description for the standard two band tight-binding model of graphene with those following from a length gauge one. These results are shown here in Figure 1. One can see that the perfect agreement between the two procedures.

Conclusion
For the past twenty years, the length gauge has been the formalism of choice for computing the nonlinear optical response of crystals while its counterpart formalism -the velocity gauge Figure 1. Frequency dependence of the third-order nonlinear conductivity of graphene, normalized to σ 0 at zero temperature. The parameters used here were (a) µ/t = 0.1, γ/t = 0.011; (b) µ/t = 0.1, γ/t = 0.001, for γ and t the scattering and hopping parameters, respectively. The solid curves were obtained from a length gauge approach and the dashed ones from a velocity gauge calculation. Image and caption were taken from [5].
-was considered to be of little use for practical calculations. In our recent works [4,5], we have described a more workable representation of the perturbation of the length gauge; related the observables in the length and velocity gauges; understood how the sum rules that ensure the equivalence between procedures are broken; and more importantly, we developed a working formalism for calculations in the velocity gauge.