{'fraction_non_alphanumeric': 0.02632901215958195, 'fraction_numerical': 0.008575352560881653, 'mean_word_length': 4.743362831858407, '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': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear optical responses of crystalline systems: Results from a velocity gauge analysis

In this work, the difficulties inherent to perturbative calculations in the velocity gauge are addressed. In particular, it is shown how calculations of nonlinear optical responses in the independent particle approximation can be done to any order and for any finite band model. The procedure and advantages of the velocity gauge in such calculations are described. The addition of a phenomenological relaxation parameter is also discussed. As an illustration, the nonlinear optical response of monolayer graphene is numerically calculated using the velocity gauge.


I. INTRODUCTION
A common approach to understanding nonlinear optical effects in atomic and condensed matter systems comes from a perturbation theory where the electric current is expanded in powers of an external applied electric field 1 , assumed to be sufficiently weak for the expansion to be physically meaningful. In this framework, one attempts to calculate nonlinear optical response functions, usually second or third order susceptibilities, the different frequency components of which describe a variety of physical phenomena like the Kerr effect, harmonic generation, electro-optic effect, etc.
The theory was first developed for atomic systems and, in the early nineties 2,3 , extended to crystalline systems, characterized by electronic bands and the corresponding Bloch functions. The simpler calculations follow two essential assumptions which will be adopted throughout this work: the independent particle approximation, where explicit electron-electron interactions are disregarded (except possibly in the determination of the electron bands, as in Density Functional Theory), and the electric dipole approximation, where the response functions are taken to be local in space, a consequence of the long wavelength limit of the applied electric field.
Even in this simple approach, however, difficulties were found regarding different ways to describe the perturbation. We can write the complete Hamiltonian of the crystal under the influence of the external field in two ways, either in the length gauge, or in the velocity gauge or minimal coupling Hamiltonian, where H 0 is the unperturbed crystal Hamiltonian, q = −e is the charge of the electron and A(t) the vector potential, chosen so that E(t) = −∂ t A(t). These two descriptions can easily be shown to be equivalent and related by a time-dependent unitary transformation 4 . The first attempts at computing nonlinear optical responses in crystals made use of the minimal coupling method 2,5 , since it has the advantage of retaining the translation symmetry of the Hamiltonian, only coupling states with the same Bloch vector k. However, it presented some serious difficulties as it seemed that response functions computed in the velocity gauge diverged in the limit of low frequencies, even for the case of zero temperature insulators, where such divergences should clearly be absent. A more detailed analysis of the linear response showed that these were only apparent divergences that could be removed by careful manipulations and sum rules 2 . As presented, these procedures were cumbersome and not easily generalizable to higher order response functions.
This led to a preference for the length gauge in nonlinear optical response calculations 3,6-9 . This gauge presents its own difficulties, the most notorious one being that the matrix elements of the position operator in the Bloch basis are ill-defined until the thermodynamic limit is taken and, even then, they can be defined only as a distribution. Inspired by Blount's work 10 , Aversa et al. 3 pointed out that the position operator will appear in the calculations only inside commutators whose matrix elements are well defined and successfully developed the theory in the length gauge. This formalism has since been applied to various systems 6,8,9,[11][12][13] .
More recently, and still within the length gauge approach, simple expressions for the nonlinear conductivity to arbitrary order in the electric field were derived by the present authors, making use of the covariant derivative notation 4 . These expressions are useful for inspection and reduced the calculation procedure to a fairly straightforward expansion of commutators. Although the idea of a "generalized derivative" was already around in the literature 3 , it involved separating the intra and inter-band components of the position operator and response functions. The emphasis on these distinctions, originally motivated by the intent of applying the equations to the special case of clean, cold semiconductors and to make analogies with atomic and free electron systems, made the calculations less transparent, in our view.
In this same work 4 , similar and equivalent expressions were derived using the velocity gauge, for the Hamiltonian where the second term is the periodic lattice potential. It was shown that, in the form derived from the velocity gauge, the expressions for the nonlinear optical coefficients lose their validity if only a finite number of bands around the Fermi level are taken into account. This difficulty was recognized early on 3,14 , and, together with the apparent infrared divergences, led to the velocity gauge being less adopted. The reasons for the failure of these expressions under a truncation in band space were also understood 3,15 and recently subjected to a more quantitative analysis 16 . One of the arguments 3 consisted in noting that the sum rules that connected the expressions in the two gauges seem to rely on commutator identities such as (r and v, position and velocity operators), which require an infinite number of bands to hold. This led to a common misconception that the velocity gauge could only be properly implemented if an infinite number of bands is taken into account 3,8,14,[16][17][18] . In fact, sum rules of greater generality have been constructed which remain true even under truncation of the bands 4 . Nonetheless, various authors have indicated that the velocity gauge, if employed, would lead to different, unphysical, predictions 16,17 . The fundamental difference between the two gauges that seems not to have been properly appreciated concerns the form of the perturbation. In the velocity gauge, the form of the perturbation depends explicitly on H 0 , unlike in the case of the length gauge. Expressions for the nonlinear coefficients derived from the Hamiltonian in Eq. 3, cannot be truncated to a finite set of bands. A truncation of H 0 implies a different form of the perturbation, and leads to different expressions for the nonlinear conductivities in the velocity gauge, but not in the length gauge.
In order to build the theory in its most general form, we will make no assumptions on the form of H 0 , other than it has the periodicity of some Bravais lattice, so that Bloch's theorem applies and there is a well defined First Brillouin Zone (FBZ). The derived forms for the nonlinear conductivities will be completely equivalent to the ones obtained from the length gauge and can be applied to any finite band model; no commutator identities of the kind of Eq. 4 are assumed.
In the next section, the equivalence of the two gauges is revisited, starting from an Hamiltonian with a finite set of bands, and some formal concepts behind the velocity gauge formulation are reviewed. In section III, the equation of motion for the reduced density matrix in the velocity gauge is introduced in a form more general than previously presented. Recursively solving the equation of motion leads to the nonlinear conductivities. The advantages and subtleties of using this gauge are discussed in Section IV, where it is argued that the velocity gauge should prove more efficient for numerical computations. In Section V, the introduction of a phenomenological relaxation parameter is addressed, and proves to be less trivial than expected. As an example, in Section VI the formalism is applied to third harmonic generation in graphene, in a full Brillouin zone calculation, and a comparison is made with already existing results in the literature. Section VII contains some closing remarks.

II. DENSITY MATRIX FORMALISM
In very general terms, which do not exclude finite band models, the single particle unperturbed Hamiltonian H 0 can be written as where |ψ ks are the Bloch band states and [H 0 ] kss ≡ ks δ ss , with ks the band energies.
To represent the scalar potential, −q E · r, in the Bloch basis, we require from the start the infinite volume limit. We define the normalization of the Bloch wave functions as where d is the dimensionality of the system. Using Blount's results for the matrix elements of the position operator 10 , where ξ kss is the Berry connection 19,20 , one obtains for the single particle perturbed Hamiltonian in the length gauge 4 (see Appendix B) where the covariant derivative, D kss , is defined by In our previous paper 4 , we discussed the equivalence of the length and velocity gauges, starting from a theory with an infinite number of bands (Eq. 3); we showed that a suitable truncation of final expressions for the nonlinear optical response functions to a finite set of bands leads to a reasonable approximation only in the length gauge. Therefore, if we want to formulate correctly a velocity gauge calculation with a finite set of bands, we should take Eqs. 8 and 9 as our starting point, and obtain the single particle velocity gauge Hamiltonian, H A , from a time dependent unitary transformation of H E . In this fashion, we preserve the equivalence of both descriptions, even when the sum over band indexes is finite. In appendix B we derive with the time dependent perturbation, V kss (t), being expressed as a power series in the vector potential, where h α1...αn An implicit summation over repeated Cartesian components α i is henceforth assumed. The coefficients in the expansion are written explicitly in Eq. 13 in terms of commutators involving covariant derivatives. If we take the coefficient of the first order term in the expansion of Eq. 13 as an example, These are the unperturbed velocity matrix elements, since v (0) (15) Had we started from the Hamiltonian of Eq. 3, the perturbation expansion of Eq. 12 would be reduced to the linear term in A(t); the second order term would have been a k independent constant, irrelevant for the dynamics, and all higher order terms would have been zero 4 . But to proceed in more general terms and ensure equivalence between length and velocity gauges for finite band models, we must, for now, retain all the terms in the series.
For the analysis of the steady state response of the electric current density, J, it is useful to rewrite the perturbation of Eq. 12 in frequency space, where the connection can already be made with the electric field components by E(ω) = iωA(ω), Having carefully defined the perturbation in the velocity gauge, we can turn to the nonlinear optical response functions, defined in frequency space by a similar expansion, The ensemble average of the electric current density, as for any observable, is given in terms of the density matrix ρ, where the reduced density matrix (RDM) is defined as a the expectation of a product of a creation and a destruction operator in Bloch states: The standard density matrix formalism computes the nonlinear conductivities by expanding the RDM in powers of the electric field. In the velocity gauge, however, the electric current is an explicitly time and field dependent observable. The velocity matrix elements are defined by −i −1 r, H A = −1 D, H A and also have to be expanded in powers of the electric field: The expansion must therefore be done in the density matrix and the velocity matrix elements simultaneously. In the absence of an external field, the current is The first order response is and, in general, The expansion of the velocity matrix elements in the external field is already defined in Eq. 20. The expansion of the reduced density matrix involves solving its equation of motion recursively.

III. RECURSIVE SOLUTIONS TO THE EQUATION OF MOTION
In the absence of scattering, the equation of motion for the reduced density matrix is We can isolate the perturbation on the right hand side, where ∆ kss ≡ ks − ks .
To solve the equation of motion perturbatively, it will be written in frequency space, order by order in the electric field, using the expansion of Eq. 16. For every order n, the reduced density matrix will be expressed in terms of its lower order terms. To alleviate notation and make the recursion relation clearer, we factorize the electric fields and other common factors by defining, The goal is now to express the recursion relation between objects of the form ρ α1...αn kss (ω 1 , . . . , ω n ). In the absence of a perturbation, the reduced density matrix is simply the Fermi-Dirac distribution, which, when replaced in Eq. 21, implies J (0) = 0, as expected.
The first order term is and the second order, The pattern is already becoming clear. As an additional example, the third order term has the form, Finally, to general order n, the perturbative solution to the equation of motion is recursively expressed as ρ α1...αn kss (ω 1 , . . . , ω n ) = 1 This recursion relation can be unfolded into lengthy expressions and its structure analyzed in more detail. However, we shall see that the real value of these expressions lies in their numerical evaluation, for which a recursion relation is sufficient. Applying Eq. 20 and Eq. 26 to the Eq. 23, the general form of the nonlinear optical response functions in the velocity gauge is obtained, This form of the nonlinear optical response functions will still have to undergo the usual symmetrization procedure to ensure intrinsic permutation symmetry 11 . Albeit trivial, this last step is a bit cumbersome to write down and will be left implicit.These expressions are equivalent to the ones we derived in a previous work 4 using the length gauge. Although far more complicated, they have their advantages, which we will discuss in the next section. The equivalence of the results of the two approaches, length and velocity gauges, is guaranteed by their being related by a unitary transformation (see Appendix B). It can also be explicitly shown by using very general sum rules to map the expressions for the nonlinear conductivities in the velocity gauge to those of the length gauge, order by order, as shown in an appendix of our previous work 4 . The proof for first order responses is nonetheless presented in appendix A, as an example of these sum rules, which generalize those derived in earlier works 2 , no longer rely on commutator identities (Eq. 4), and are valid for a model with a finite number of bands.

IV. LENGTH VS VELOCITY GAUGE
As usual, there are advantages and disadvantages associated with any particular choice of gauge. By considering the (exactly equivalent) forms of the nonlinear conductivities derived in the two gauges, the strengths and weaknesses of each can be analyzed.
A first look at Eq. 32 will immediately bring out the usual concerns with infrared divergences in the velocity gauge form, due to all the inverse frequency factors. We emphasize again, however, that this expression is equivalent to the one obtained from the length gauge and therefore these divergences are only apparent. Manipulating the expressions in the velocity gauge and using a series of sum rules, it can be shown that Eq. 32 is divergence free. This approach was the one originally pursued 2 , but this use of sum rules became rather pointless after the length gauge formulation had been developed 3 . If the sum rules are employed in the velocity gauge to remove apparent divergences, one will simply arrive at the same expression as obtained more straightforwardly in the length gauge. This equivalence through sum rules 4 does not demand an infinite number of bands, but it does put a constraint on the use of Eq. 32, namely that the integration must be done over the full FBZ to cancel divergences. An effective continuum Hamiltonian describing a portion of the FBZ-like the Dirac Hamiltonian for graphene-, will not suffice, since these sum rules rely on the periodicity in k space of the quantities involved.
Having clarified this point, it can still be noted that the velocity gauge form is considerably more elaborate; less useful not only for inspection, but in an actual analytical calculation. As an example, the expression of the second harmonic generation with all components along the x axis, in the length gauge 4 , is 2 , while in the velocity gauge, 2 The symbol • stands for the Hadamard product: where we still have to write the density matrix components, This simple example illustrates that there is no advantage in doing the analytical calculations in the velocity gauge, although inspection of previous equations shows an interesting point: there are only simple poles in the velocity gauge form ( ω − ∆ ) −1 , while in the length gauge, by differentiation, higher order poles emerge. Still, for analytical calculations, we would advocate the more transparent and simpler length gauge approach 3,4 .
The strength of the velocity gauge form lies in the different arrangement of the commutators. The covariant derivatives are no longer applied to the density matrix in its recursion relation 3 . Instead, they operate only on the unperturbed Hamiltonian H 0 in the determination of the functions h kss (Eq. 13), which are independent of frequency, temperature and chemical potential.
A careful look at the algorithm of the previous section, shows that for the nonlinear conductivity of order n, there are n + 1 such functions to compute by successively applying a covariant derivative: h α1...αm kss with m = 1, ..., n + 1. In the previous example of second harmonic generation, these would be h x kss , h xx kss and h xxx kss . Further reducing this algorithm to its fundamental ingredients, we note that these calculations demand only a knowledge of two objects, which fully define the system under consideration: the dispersion relation ks and the Berry connection ξ kss .
Once these h α1...αn kss functions are analytically determined, the integrand in Eq. 32 or in Eq. 34 can be numerically evaluated at each k value, independently and quite easily. In fact, the procedure involves evaluating the analytic h α1...αn kss functions and the Fermi-Dirac distribution at the k point and then computing simple commutators and traces of numeric matrices. It has no numerical derivatives at all. This is in contrast with the length gauge, where either the full expression of the response function is analytically calculated or numerical derivatives have to be applied in each step of the density matrix recursion relation. Either way, via the product rule and higher order poles, the number of complicated terms to evaluate grows very fast with n in the length gauge approach.
For this reason, the form of the nonlinear conductivity in Eq. 32, derived from the velocity gauge, should provide a more efficient numerical approach. The authors have implemented numerically the expressions in both gauges and done calculations on the nonlinear conductivity of monolayer graphene and observed that the computation times were indeed considerably smaller when Eq. 32 was used. The velocity gauge results are presented in Section VI.

V. PHENOMENOLOGICAL RELAXATION PARAMETERS
In Eq. 32, like in all previous equations, the addition of the usual infinitesimal imaginary part to the frequencies, ω → ω + i 0 + , is always implicit, as imposed by causality. This provide us with well defined relaxation free expressions. From the numerical standpoint, the imaginary part must always be finite, but can be taken to be smaller than any other energy scale in the problem.
However, one is also interested in modeling relaxation processes due to scattering of electrons with impurities, phonons and other electrons. A simple phenomenological approach is to add one or more relaxation parameters to the frequency poles. A common justification for the addition of a relaxation parameter γ comes from considering a scattering term in equation of motion 21 , If the perturbation is turned off, the density matrix relaxes to the equilibrium distribution ρ eq . In the length gauge approach, ρ eq = ρ (0) . A simple rearrangement 4 , makes clear how this will impact the nonlinear conductivities. The second term on the right hand side of Eq. 38 is not relevant to any order n = 0, implying that the only alteration will be to add an imaginary constant to each pole in frequency space, Again, using the second harmonic generation as an example, the scattering term will induce the following simple modification of the expression in Eq. 33, In the case of the velocity gauge, this change is much more complicated, since the equilibrium distribution depends on the perturbation. Having in mind the connection between the two gauges, it is easy to see that the equivalent ρ eq in the velocity gauge equation of motion should be obtained by an unitary transformation 4 of the Fermi-Dirac distribution 5 . This again translates into an expansion in the vector potential, If this equilibrium distribution is used, the results obtained will, indeed, be equivalent to the ones from the length gauge. Nevertheless, as discussed in the previous section, the main advantage of using the velocity gauge comes from the absence of derivatives acting on the density matrix, and once the term in Eq. 41 is added to the equation of motion, that advantage, and the greater numerical efficiency of this approach, are lost.
As a consequence, we will follow a different phenomenology, more appropriate for our computations. To each frequency ω j we shall add a constant imaginary part 6 : ω j → ω j + iγ. This method comes naturally from considering the adiabatic switching of the external fields. It looks similar to the previous method, but the expression for second harmonic generation (as would be obtained in the length gauge) has now a factor of two in the relaxation parameter that appears in one of the poles. This might seem like a slight difference, but is a distinct phenomenology and, as we shall see, even for low values of γ it can lead to completely different results for very low frequencies ω γ and in a region of width γ around resonances.

VI. NUMERICAL IMPLEMENTATION: THIRD HARMONIC GENERATION IN GRAPHENE
In this section, the velocity gauge approach is tested numerically, by computing the third order conductivity of a material whose nonlinear optical properties have been the subject of intensive research in recent years: monolayer graphene 11,12,18,[22][23][24][25][26][27] .
Monolayer graphene is a two dimensional sheet of carbon atoms, displayed in an honeycomb lattice. It has two atoms per Bravais lattice site. Here, we shall consider only the simplest nearest neighbor tight binding model that describes the electronic properties of graphene 28 , with t being the hopping parameter and a 1 = (1/2, √ 3/2) a and a 2 = (−1/2, √ 3/2) a the basis vectors of the Bravais lattice. From this model Hamiltonian, the dispersion relation and the Berry connection can be computed, ks = s t |φ(k)| (45) with band index s = ±1 (−1 for valence and +1 for conduction band) 7 . The knowledge of the dispersion relation and the Berry connection is sufficient for a calculation of the nonlinear conductivity, independently of the gauge. It can be shown that a single independent component σ (3) xxxx (ω, ω, ω) describes third harmonic generation. We shall confine ourselves to the study of this frequency component of the third order nonlinear conductivity.
Results obtained from the standard length gauge approach are presented in Fig. 1 Anomalous feature of the third order nonlinear conductivity of graphene from Fig. 1(b), in a region near the three photon resonance 3 ω = 2µ.
Dirac point. In this case, the analytical form of the third order conductivity in the Dirac point approximation was calculated and then plotted. Also, γ is introduced via the additional scattering term in the equation of motion (first type of phenomenological treatment described in the previous section), as in Mikhailov's work 21 .
These results are in agreement with analogous calculations already in the literature 12,18,21 , with the strongest feature present at the three photon resonance 3 ω = 2µ. In particular, Fig. 1(a) is identical to Fig. 3(b) and Fig. 5 in refs 18 and 21 , respectively. Fig. 1(b) shows a very anomalous behavior at the resonance, which is always present in a region |3 ω − 2µ| ≤ γ for small but finite γ. Fig. 2 shows a close-up of this region. This strange feature near resonance is analyzed in detail in 21 , where it is regarded as a prominent feature of potential interest, since in practice γ is always finite. However, a more careful analysis shows that despite the unusual shape of this feature, in the limit γ → 0 + , the scattering free result in 18 is indeed obtained (see Appendix C).
To do the same calculations with the velocity gauge approach developed here, we evaluate analytically not the full third order conductivity but only the following functions All these analytical functions follow from direct evaluation of the commutators (Eq. 13) using Eqs. 45 and 46. At this point, the expression in Eq. 32 for the nonlinear conductivity (with n = 3) can be evaluated by numerical integration over the full FBZ. The phenomenology adopted here is the one that follows from adiabatic switching, for the reasons discussed in the previous section.
The velocity gauge results are in Fig. 3. The curves are markedly different than the ones obtained from the length gauge, especially for large γ. This should came with no surprise, since a different phenomenological treatment is adopted. To prove that the difference between the two curves is only due to the way the relaxation parameters are introduced, the length gauge calculations were performed again, now with this second type phenomenological treatment (using the third order analogue of Eq. 42) and included in the plots of Fig. 3. The results obtained in the two gauges are identical.
Of course, as we take the relaxation free limit γ → 0 + the results of both phenomenological treatments also coincide. However, for finite γ it may lead to some considerably different behavior of the nonlinear conductivity. In particular, we note that in the phenomenology adopted here ( ω i → ω i + iγ), those anomalous features seen in Fig. 1 and discussed in 21 are not present. Instead, a more plausible smooth curve is seen at the resonances.
As a final remark, we emphasize that the results presented here for the velocity gauge involved a complete two-band tight binding model and an integration over the full FBZ, not an effective Dirac Hamiltonian, as in the case of the length gauge results. This has some consequences. First, it had been suggested that a possible source for the two order of magnitude discrepancy between theoretical and experimental results in graphene could be the use of an effective Hamiltonian 18 . The agreement displayed in Fig. 3 establishes that the Dirac point approximation is valid for the range of frequencies adopted in previous studies 12,18,21 . Second, the use of complete bands will allow us to proceed towards higher frequencies and study spectral regions where the Dirac Hamiltonian does not give an accurate description.

VII. CONCLUSIONS
In summary, a velocity gauge approach to calculations of nonlinear optical conductivities was developed in this work, within the density matrix formalism. It was shown that, contrary to common belief, the results are the same as those from the length gauge, as demanded by gauge invariance. No difficulties come from applying this velocity gauge formalism to finite band models.
The velocity gauge provides an efficient algorithm for computing nonlinear conductivities. The dispersion relation and the Berry connection of a crystal are the prerequisites and define the crystalline system under study. A series of commutators (Eq. 13) can then be analytically computed. From these, the conductivity can be numerically evaluated for any frequencies, temperature, chemical potential, relaxation parameter and to any order, without using any numerical derivatives.
Previously, studies of nonlinear conductivities often involved writing down analytical expressions for the full conductivity. For third order, this already becomes rather cumbersome. Analytically, it is only really tractable for simple effective continuum Hamiltonians (such as the Dirac Hamiltonian). The velocity gauge approach developed here does not have such complications, is easily generalizable to higher orders and is always applied to complete bands. This gauge provide us with a simple yet powerful method to compute the nonlinear optical response functions of crystalline systems.

Appendix A: Equivalence of linear responses
As an illustration of the general sum rules implicit in a velocity gauge treatment, the expression for the linear conductivity in the velocity gauge (A1) will be shown to be equivalent to the one obtained from a length gauge approach, On first inspection, this equation might suggest that this operator is diagonal in Bloch momentum; it is not because of the presence of the gradient with respect to k. Any observable that can be written as a differential operator acting on the wave function in a continuous basis will have a similar representation 8 .
The full single particle Hamiltonian in the length gauge is, therefore, The time evolution in this gauge is The second term is seen to cancel the scalar potential term in the Hamiltonian, The velocity gauge Hamiltonian therefore becomes