A comprehensive study of the velocity, momentum and position matrix elements for Bloch states using a local orbital basis

We present a comprehensive study of the velocity operator, $\hat{\boldsymbol{v}}=\frac{i}{\hbar} [\hat{H},\hat{\boldsymbol{r}}]$, when used in crystalline solids calculations. The velocity operator is key to the evaluation of a number of physical properties and its computation, both from a practical and fundamental perspective, has been a long-standing debate for decades. Our work summarizes the different approaches found in the literature, connecting them and filling the gaps in the sometimes non-rigorous derivations. In particular we focus on the use of local orbital basis sets where the velocity operator cannot be approximated by the $k$-derivative of the Bloch Hamiltonian matrix. Among other things, we show how the correct expression can be found without unequivocal mathematical steps, how the Berry connection makes its way in this expression, and how to properly deal with the two popular gauge choices that coexist in the literature. Finally, we explore its use in density functional theory calculations by comparing with its real-space evaluation through the identification with the canonical momentum operator. This comparison offers us, in addition, a glimpse of the importance of non-local corrections, which may invalidate the naive momentum-velocity correspondence.


I. INTRODUCTION
The quantum mechanical velocityv operator plays a central role in the evaluation of macroscopic optoelectronic properties of crystalline solids. The velocity matrix elements (VMEs) are generically needed to determine transitions between band states through several formalisms such as Fermi's golden rule [1] for decay and optical excitation processes or the more general Kubo linear response theory [2]. Closely related tov, the canonical momentum operatorp also plays a key role, but more from a methodological standpoint. The momentum matrix elements (MME) are, for instance, needed to find parameter-free effective models within k ·p perturbation theory [3].
It is common to consider the velocity operator aŝ v = i[Ĥ,r] (in atomic units, which we will use throughout the text), following a classical to quantum mechanics identification through the Heisenberg equations of motion (Heisenberg picture). In the cases whereĤ only contains the kinetic energy and a potential commuting with the position operatorr (i.e., in absence of spinorbit coupling or non-local potentials), to work withv or withp becomes completely equivalent. In coordinate representation this means that i[H, r] and −i∇ r are interchangable.
In the following, we review the state of the art of the uses and misuses of these two operators as well as the position operatorr when evaluating matrix elements between band states. Evaluating the commutator matrix elements presents no problems when dealing with localized states, as in atomic physics, but fundamental difficulties can be found when dealing with Bloch eigenstates due to rψ nk (r) not belonging to the same Hilbert space as that of the states themselves, ψ nk (r). This issue has been addressed by Gu and coworkers [4], in addition to presenting an extensive review of the existing p-r relations (as called in their work) in the literature. Gu et al. find the correct relation between the momentum (or velocity for the case whenv andp are equivalent) and position matrix elements: Here the matrix element is taken between eigenstates normalized to a finite volume v and obeying periodic boundary conditions (PBCs), therefore not decaying at the boundaries even if the volume tends to infinity. The surface term C is calculated on the surface of the solid. With this is mind, one can easily convince oneself that the position (also called dipole) matrix element depends on the origin of coordinates and that C compensates this choice (as the VME cannot depend on the origin). Eq. (1) above presents a fundamental view of i[Ĥ,r] rather than convenient shortcut to evaluate the VME, as a challenging integration in coordinate space is needed on the right side of the equality. It also remarks the difference with the atomic case, where the surface term does not appear.
In practice, Bloch eigenstates are very often represented in a Bloch basis which, in turn, may be expanded in a local orbital basis. In this regard, a good effort has been put in the actual evaluation of the VME over the last decades. Between same-k Bloch eigenstates, the arXiv:2201.12290v2 [cond-mat.mes-hall] 31 Jan 2022 VME can be calculated through the following expression: where the c's are the coefficients of the expansion of eigenstates in a generic |αk v Bloch basis and A αα (k) is the Berry connection associated with such Bloch basis. The first term is sometimes referred to as the Peierls approximation, while the second term is needed to deliver the full matrix element. The suppression of the second term, as we will show later, can lead to large and uncontrolled quantitative errors. This issue was first explored by Pedersen et al. [5], by trying to complete the Peierls approximation within a tight-binding scheme. We also point out the remarkable work by Tomczak et al. [6], introducing an intra-unit cell correction to ∇ k H αα (k) (we will see later that this can be understood as a consequence of the gauge choice), and adding extra terms to this quantity similar to those in Eq. (2). Tomczak et al. contribution was later replicated, perhaps in a clearer and more complete way, in the work of Lee et al. [7], who presented a complete expression similar to that in Eq. (2) for a general, nonorthonormal basis. Actually, as originally reported in Ref. [7], the Berry connection did not appear. The fact that Eq. (2) can be recasted in this form will be shown below in this work, thus generalizing the evaluation of the VME to any basis, not necessarily comprised of local orbitals. In the light of Eq. (1), Eq. (2) presents a somewhat puzzling aspect: first, there is no surface term and, second, no term depends on the placement of the integration volume v. However, Eq. (2) was derived in Ref. [7] essentially in the same manner as Eq. (1) was derived in Ref. [4], namely, by making use ofv = i[Ĥ,r] projected in a chosen basis, the coordinate basis in the former and a local orbital basis in the latter. Addressing this intriguing observation is part of our motivation to carry out the present work, as well as studying the role of ∇ k H αα (k) in the calculation of velocity matrix elements. We will also explain how the the popular expression for position matrix elements given by Blount [8] fits into this comparison.
To this end, we organize this study as follows. In Sec. II we present the main theoretical ingredients by first recalling the differences between periodic boundary conditions versus the infinite volume case when defining Bloch eigenstates. We follow by introducing two ways of treatingv = i[Ĥ,r], one relying on an integration in the whole finite volume of the system, and the other based on using the k representation for operators, involving matrix elements between the cell-periodic part of the Bloch eigenstates. In Sec. III we explain how Eq. (2) rigorously comes about from the second method, while showing the way it has been previously derived in the literature is mathematically inconsistent, to say the least. In Sec. IV we present a numerical study that gives us insight into the quantitative error that one makes when assuming the equality nk|v|n k = nk|p|n k , even in the presence of a local potential, and the trade-off between computational simplicity and accuracy when using the Peierls approximation in a practical situation. Finally, we summarize our main conclusions in Sec. V.

A. Preliminary definitions
We start by recalling the normalization choices for eigenstates in a crystal. This turns out to be a key point to understand the relation between position and velocity operators. Two options are compatible with Bloch theorem: one can assume a finite volume v normalization or let the eigenstates extend to all space. In order to distinguish the two cases, we write being the normalization conditions for every case case Note that this convention assumes |u nk to be normalized to one in the finite volume case and to the unit cell volume (denoted with Ω in the following) in the distribution case. The v[x0] means that the integration domain (the volume v) is defined by the parallelepiped determined by N i a i with N i being the number of cells in each direction given by the primitive vectors a i . Its origin vertex is located at the x 0 point, which have to be selected in the first place. Finally N = N 1 N 2 N 3 is the total number of cells of the crystal. Born-von Karman boundary conditions ψ nk (r + N i a i ) = ψ nk (r) are applied, leading to a quantization of the crystal momentum according to k = l1 N1 G 1 + l2 N2 G 2 + l3 N3 G 3 with G being reciprocal lattice vectors and l i = −M i , . . . , M i such as N i = 2M i + 1. In general, matrix elements of operators whose application on eigenstates breaks periodicity may depend on x 0 . On the other hand, in the case of infinite volume normalization, the integrals run over the whole unbounded space, including infinities. In this case, k vectors become a dense set inside the Brillouin zone (BZ).
The representation of a given operatorÔ in both cases becomes (we will assume that is equivalent to ∞ −∞ in the following). If the operatorÔ is such that O(r)ψ n k (r) still satisfies Bloch theorem, the matrix elements can be reduced to an integration within the unit cell involving the periodic part of eigenstates whereÔ k ≡ e −ik·rÔ e ik·r , sometimes called the "k representation of an operator". In the infinite volume case the expression is similar but taking δ kk → δ(k − k )/Ω. In what follows we particularize to the velocity operator v and its relation to other quantities.

B. Relation between the velocity and momentum matrix elements
As discussed in the introduction, the velocity and momentum operators can only be interchanged if spin-orbit coupling is neglected and the periodic potential in the crystal is assumed to be local. Let the Hamiltonian be separated intoĤ =Ĥ L +Ĥ , whereĤ L =p 2 /2 + V (r) with V (r) the local periodic part of the lattice potential and where no spin-orbit coupling has been included. Then, for the first term one can writep = i[Ĥ L ,r], and the projection of the full velocity operator on a subspace of band states can be written as nk,n k .
nk,n k is an extra term that one encounters in practice and which is a consequence of the finite nature of the Hilbert space in actual calculations. This can be ultimately traced back to the fact that the canonical commutation relation [r α ,p β ] = iδ αβ can never exactly hold for a finite-matrix representation since in such cases Tr(r αpβ ) = Tr(p αrβ ). Therefore, one can only expect the quantity ∆ (v) nk,n k to be negligible if the physical states are sufficiently well represented in the working Hilbert space. We will give below a few examples of this practical limitation.
The presence of the second term in Eq. (7) is challenging from a practical standpoint and only whenĤ = 0, Eq. (7) becomes the theoretical velocity-momentum equality. In any case, evaluating the VME seems to require, in principle, the evaluation of the MME through its representation -i∇ r . In the following we explore two routes that can be followed to by-pass the evaluation of the MME and, at the same time, of the second term if present.

C. Relation between velocity matrix elements and the Berry connection
We first write the VME in the k representation. It is easy to see thatv k = e −ik·r i[Ĥ,r]e ik·r = ∇ kĤk , so one can write By applying chain rule it is straightforward to find where ω nk,n k ≡ n (k) − n (k) and A nn (k) ≡ i u nk |∇ k u n k Ω , this last quantity being the Berry connection [9]. Eq. (9) can be found in the literature, see e.g. Ref. [10]. The equation above replaces Eq. (7) by introducing the evaluation of the Berry connection associated with the Bloch eigenstates. This, however, can be a cumbersome task since k derivatives of eigenstates are not known in numerical diagonalization procedures. While this problem can be circumvented through perturbation theory [11], in Sec. III we show how Eq. (9) can be recasted in a more convenient and familiar form. Incidentally, note that Eq. (9) provides a way to compute the non-diagonal Berry connection elements if the VMEs are known.

D. Relation between velocity and position matrix elements
Alternatively, we can directly perform the integrals that appear in Eq. (7) when representing in coordinate space. AssumingĤ = 0 and, therefore, being able to writev ≡ i[Ĥ,r] =p, one is free to use −i∇r or i[Ĥ,r]. In both cases the explicit knowledge of the real-space wavefunction of the eigenstates is required. In the former case derivatives need to be carried out, which depending on the orbital basis can be more or less cumbersome to implement. In the latter, the use of the commutator entails further steps, where one needs to pay attention to the correct use of the hermiticity ofĤ in theĤr product. This procedure, which has been followed by Gu and coworkers in Ref. [4], only applies to eigenstates in the framework of finite volume normalization, where integrals for matrix elements can be converged. One starts n k (r). After performing the derivatives and using Gauss's theorem, one obtains (same comment [9] applies here). The appearance of this last surface term arises from the periodicity of the wave functions at the surface of the material volume, which we denote with ∂v[x 0 ]. It is important to note that the wavefunctions do not decay even in the limit of an infinite volume and this term is always present.
As noticed in Ref. [4], the hermiticity property cannot be applied as usual in nk|v|n k v = i nk|[Ĥr −rĤ]|n k v , which results in the surface term above. Secondly, both matrix elements on the right hand side (RHS) in Eq. (11) depend on the origin of the integration volume and are not k-diagonal, while the sum does not depend on this arbitrary choice of origin and is diagonal in the wave vector as the VME actually is. The relative weight of nk|r|n k v versus C nk,n k with respect to the full VME is also explored in Ref. [4] showing that, in general, one cannot find a point x 0 that makes the surface term to vanish, even for certain analytical models. Therefore, Eq. (11) presents no advantage versus directly computing −i nk|∇r|n k to find the VME, as one still has to perform nontrivial integrations for position and surface matrix elements. It provides us, however, with the conclusion that momentum and dipole matrix elements (multiplied by the frequency) should never be interchanged when dealing with Bloch eigenstates in a finite volume.

E. Relation between velocity and position matrix elements with a distribution basis
If Bloch eigenstates are normalized as distributions [recall Eqs. (3) and (4)], then one can still use them as a basis to represent general physical quantum states in the crystal. This was originally referred to as the crystal momentum representation (CMR) [8], where one writes with g n (k) being the envelope function for the n band. The matrix elements between two physical states is written n (k) nk|Ô|n k .
(14) Now one needs the find the matrix element nk|Ô|n k that enters the calculation above. As only the full n sums and k integrations are relevant, we can take into account the boundary properties of the state |φ . This is the case of the matrix elements for the position operator, for which Blount [8] noticed that nk|r|n k is ill-defined by itself but that a distribution form can be given ifr is assumed to act on a state |φ belonging to its domain. Specifically, Blount showed that The effect of boundary conditions on eigenstates is highlighted here, as nk|r|n k fundamentally differs from nk|r|n k v , addressed in the previous section. Position operator matrix elements here do not depend on any arbitrary origin, but only make sense within Eq. (14). As far as the velocity operator is concerned, Eq. (9) is still perfectly valid in the infinite volume case: expression to be used, again, only in the context of Eq. (14). Alternatively, in Appendix A we also show that projectingv on general physical states φ|v|φ = i φ|[Ĥ,r]|φ , along with Eq. (15), also leads to Eq. (16). Finally, to complete the connection between the different matrix element expressions, it is straightforward to show that Again, one simply needs to project φ|v|φ = i φ|[Ĥ,r]|φ and proceed in the same manner as explained in previous section. Here, however, the surface term vanishes due to |φ and |φ being square-integrable over all space. Notice that this momentum and position relation matches that in atomic physics. This is also shown in Ref. [4] by using narrow k envelope functions in the limit of zero width.

III. VELOCITY MATRIX ELEMENTS WHEN REPRESENTING IN A BLOCH BASIS.
Having established a comprehensive overview of the available recipes to evaluate the VME and their proper use, we proceed now with their actual computation when a generic and possibly non-orthonormal Bloch basis is used to expand the Bloch eigenstates: We stress again that |αk v is a generic basis state satisfying Bloch's theorem in a finite volume, with α being a generic quantum number. The coefficients c (n) α (k) are found by solving the generalized eigenvalue problem (19) where H Eq. (9) can now be properly converted into more familiar expression. First the Berry connection reads In this expression one has to perform derivatives in k space of the coefficients c (n) α (k). In most cases these coefficients are obtained by numerical diagonalization of Eq. (19) so that they are not continuous and, therefore, differentiable. However, this can be avoided by directly employing the chain rule after inserting Eq. (20) into Eq. (9), leading to We show the complete derivation in Appendix B. Eq. (21) is one important result of this work: it allows to compute VME from the Hamiltonian matrix elements and the Berry connection in whichever Bloch basis. We have differentiated two contributions, A and B, to the VME. The first one is evokes the exact expressionv k = ∇ kĤk , but the second one is equally important, as we will show below. Eq. (21) clearly shows that ∇ kĤk is not, in general, αα (k). In many practical cases the Bloch basis is expanded, in turn, in a local orbital basis. Regarding this, two different types of basis can be found in the literature: where |αR v is an orbital with d α position vector inside the unit cell of site R. A finite size crystal containing N cells is assumed throughout. Bloch eigenstates are now given by with both expansions being related by c αk . We will use the former basis in this work by default, which we refer to as the cell gauge, and make considerations related to the other one, the atom gauge, when appropriate. In this basis the matrices needed in Eq. (19) become Now we can recast the Berry connection terms in Eq. (21) into an explicit form involving the position operator so that Eq. (21) becomes identical to that reported in Ref. [7] (see Appendix B for all the details.) It is important to note that neither of the two terms in Eq. (21) is gauge independent. One can easily check how the two terms change when switching to the atom gauge, according to Eq. (22). For instance, the first term becomes v (A, atom) nn while the correction for the B term is the same with opposite sign, showing that the absolute value of the sum v nn is gauge invariant, as it should be for a physical operator.
It is worth obtaining the form of Eq. (21) in the maximally localized (Wannier orbitals [12]) or tight-binding limit, denoted here with ν. In this case, only the intraatomic dipolar matrix elements in Eq. (25) survive, leading to which is the same as Eq. (26). This tells us that, under a maximal localization condition, computing the VME using only the gradient term (the A term) in the atom gauge, therefore neglecting v (B,atom) nn , is equivalent to computing both terms (the full VME) in the cell gauge. This means that v (B,atom) nn = 0, as can be easily checked. The second line of Eq. (26) [or equivalently Eq. (27)] was presented in Ref. [13] as a "Peierls substitution approach to the case of multiatomic unit cells". Based on our previous discussion, we see that it appears naturally within the atom gauge. In the more general case of a non-orthonormal basis, both terms of Eq. (21) must be evaluated regardless of the gauge choice. We examine this more in depth in Sec. IV.
It is worth ending this section by briefly discussing the work of Lee et al. [7]. They present an expression for the VME which is, in fact, a particular case of our general expression Eq. (21) (we reproduce it in Appendix B). However, we believe that in order to reach their expression for the VME, they have inadvertently mixed Hilbert spaces. Their derivation starts fromv = i[Ĥ,r] and, briefly, they follow by projecting nk|v|n k v = i nk|[Ĥ,r]|n k v , expanding eigenstates in a non-orthonormal local orbital basis, and inserting the closure relationÎ = αR,α R S αR,α R between the product of operators. We note that their procedure is equivalent to start by writing and proceeding in the same manner. The problem of starting with Eq. (28) is that, as explained in Sec. II, this equality only holds in the case of open boundary conditions (infinite systems). This is not the case when using local orbital basis sets, where the band eigenstates obey a finite volume normalization [see Eq. (5) and Eq. (22)]. The correct result found in Ref. [7] can only be explained by the unjustified identification of nk|r|n k v + C nk,n k in Eq. (11) n k (r), which leads to Eq. (28). By doing this, the dependence on an arbitrary origin of integration is effectively removed by the new integration limits, but the integral is ill-defined. This subtle issue, which can be easily missed, is stressed by our notation in Eq. (28), where we have put the subscript v is on the left hand side but not on the right hand side of the equality.
In the next section we present some numerical examples in order to explore the details of the VME formulas in a practical situation.

IV. PRACTICAL CASES: HEXAGONAL BORON NITRIDE AND GRAPHENE
The first goal of this section is to gauge the importance of the different terms in Eq.
, by comparing between independent evaluations of the VME and MME. Particularizing to a local orbital basis case, the former can be evaluated from Eq.
(21), while the latter becomes nk|p|n k = −i αα R c (n) * α (k)c (n ) α (k)e ik·R α0|∇r|α R . Our second goal is to explore the relative importance of the two terms in Eq. (21).

A. Detailed numerical analysis of VME and MMe
We start by computing the band structure of a benchmark material. We choose a monolayer of hexagonal boron nitride (hBN), which is a sufficiently complex system to our purposes. In Fig. 1, we show: (i) a tightbinding (TB) two-band calculation for the lower (upper) valence (conduction) bands, including only first neighbour interactions between the p z orbitals of B and N atoms, (ii) a DFT calculation employing a small-core pseudopotential basis set [14] to replace the 1s 2 electrons in every atom (labelled as CREN) and, (iii) an allelectron calculation with the 6-31G* basis set [15]. The DFT calculations were performed using CRYSTAL17 [16] with the local von Barth-Hedin exchange-correlation functional [17]. The use of Gaussian type basis sets allows us to perform real-space integrations in an analytical fashion. We are not concerned here with the accuracy of the obtained gap so we have excluded the use of hybrid functionals and their possible extra non-local contributions.
Both DFT band structures are essentially similar up to the conduction band. The agreement is particularly good for the valence and conduction bands, both with a band gap of 4.55 eV, except maybe for a noticeable difference at the M point of 0.5 eV. As expected, only the more accurate all-electron calculation with a large basis can reproduce results in the literature [18]. The tight-binding parameters can be fitted to resemble one of these calculations. It is easier to obtain a better overall fit to the CREN band structure with t = 2.15 eV, as only two p z orbitals are present to reproduce the energy dispersion.
We now explore in some detail Eqs. (7) and (21). To this purpose, in Fig. 2 we show the magnitude of several quantities relevant to the band-gap optical transition along the Γ − K − M path. Fig. 2(a) shows the absolute value of the x component of the VME and the MME for the three cases shown in Fig. 1. Looking at the VME, the TB result deviates quantitatively from the other two, but not qualitatively. When comparing the VME and the MME, we observe that for the CREN basis the difference is significant, particularly near M, while that for the large basis this difference is negligible. The latter result proves that Eq. (21) is properly implemented since the evaluation of the MME is essentially analytical due to the use of Gaussian orbitals. The difference found in the former calculation can be attributed, as reflected in Eq. (7), to both the presence of the non-local pseudopotential and the difference in the size of the Hilbert space (8 bands versus 36 bands). In order to isolate the effect of the non-locality from that of the basis size, we have repeated the calculation with the all-electron basis using the nonlocal functional HSE06 [19] (not shown in the figures). In this case, the deviation between the VME and MME curves becomes appreciable, but not larger than the one for the CREN case. In summary, these results explicitly show that the VME and MME cannot always be taken as the same quantity.
In Fig. 2(b), we compare the magnitude of the kgradient term [the A term in Eq. (21)] calculated in the atom gauge for the three different cases. In the TB case, this term gives the full value for the VME. In DFT the results deviate significantly from the exact value, showing the importance of the B term in Eq. (21). The CREN basis presents a larger deviation, despite being smaller in size [see Eq. (25)]. As mentioned in Sec. III, only the maximal localization condition for the basis orbitals ensures that v (A, atom) gives the exact VME. This condition is not met in neither of the two DFT basis sets used in our calculations. We also show the result obtained in the cell gauge in Fig. 2(c). Now, not only quantitative differences appear, but also selection rules break when approaching the Γ point (here the VME must be zero according to the irreducible representations of the wave functions). Therefore, identifying the VME simply as a k-derivative of the Bloch Hamiltonian in the cell gauge can lead, not only to quantitative errors, but also to incorrect physical interpretations.

B. Optical conductivity
The calculation of an experimentally measurable quantity such as the optical conductivity can be affected by an incorrect evaluation of the VME. To show this we make use of the Kubo-Greenwood [20] expression (we do not where f nk is the Fermi-distribution occupation number and N k is the number of k points in the discretized Brillouin zone. In Fig. 3 we show the longitudinal optical conductivity, computing the VMEs within the different approximations considered in previous section. We have separated the results obtained with the small-core basis from those with the all-electron basis, as shown in Fig. 3(a) and Fig.  3(b), respectively, where we have also added the calculation with the TB model in both panels. At the bandgap frequency, the DFT and TB calculations involving the exact VME are able to reproduce the quasiuniversal behaviour [21] for a parabolic noninteracting semiconductor, yielding σ = e 2 /2 . The use of MMEs, instead of the VMES, fails for the pseudopotential and small basis case [black dashed line in Fig. 3(b)], as expected from the discussion in previous subsection. At higher frequencies the TB model underestimates the response, which is similar in magnitude for both DFT cases, the only difference being the position of the Van Hove singularity which originates in the bands at the M point (see Fig. 1). In both DFT calculations the k-gradient approximation overestimates the exact result. A calculation with the k-gradient term in the cell gaugev → v (A, cell) (not shown) gives an even larger discrepancy at all frequencies, but worse, also removes the isotropic behaviour of the conductivity tensor with σ xx = σ yy . This erroneous behaviour has been already discussed for graphene in Ref. [22] and highlights the importance of taking the k-derivative approximation for VME using the appropriate gauge. It is also worth mentioning here the work by Wissgott et al. [13]. There, the Peierls approximation in the atom gauge is tested versus the complete VME also through a conductivity analysis of transition-metal oxides. Our conclusion about the gauge choice, not explored in their work, could give a better insight about the discrepancies that are found in Ref. [13].
A direct comparison with experiments can be made by analyzing the optical response of graphene. It is known that monolayer graphene shows a quasi-constant absorbance of ∼ 2.3%, corresponding to σ = e 2 /4 , over the energy region that goes from the far-infrared to the visible spectrum where excitonic effects are negligible (< 2 eV) [23][24][25]. Therefore, in this energy range, Kubo-Greenwood DFT-based calculations are expected to give a faithful optical response. In Fig. 4 we show the optical conductivity calculated with two different basis sets, equivalent to those used for hBN. We present results for the exact VMEs and their approximated values using the MMEs. Experimental results from Ref. [23] are also shown. We can see that both basis sets give results in very good agreement with the experimental ones when employing VME. For the case of MMEs, the CREN basis set gives ∼ 0.175e 2 / , which translates in a 30 % error when comparing to the experimental curve. This result complements our previous study of hBN, showing the significant effect of non-local operators and finite basis sets when trying to replace the VMEs by the MMEs. Fig. 3 for the case of graphene. Experimentals results from Ref. [23] are shown. VME and MMe has been used to represent the velocity operator in the Kubo-Greenwood formula for the two DFT calculations (other approximations are not shown in this case).

V. CONCLUSION
We have presented a comprehensive study of the evaluation of VME in crystalline solids, as obtained from the fundamental relationv = i[Ĥ,r]. We have scrutinized all available expressions in the literature, filling the gaps and connecting them in a coherent story. We have seen that, when working in coordinate representation, one is bound to deal with a very inconvenient surface term which can be avoided by going first into the k-representation. We have obtained a general expression which contains a familiar k-derivative term plus a correction term which involves the Berry connection of the Bloch basis elements. When using local orbitals as a basis, this can be rewritten in a more familiar form (see, e.g., Ref. [7]), but whose previous derivations contain unjustified mathematical steps. We have also shown several equivalences which involve the momentum and position operators, including well-known expressions in the crystal momentum representation (nonphysical distribution basis).
We have numerically tested the validity of different approximations to the VME by computing the optical conductivity of monolayer hBN and graphene through the Kubo-Greenwood formula. In particular, we have shown that approximating the VME by ∇ k H αα (k) in a non-orthonormal basis produces significant quantitative errors and may also give rise to qualitative ones if one is not careful with the choice of gauge. We have also made emphasis on the fact that the velocity and momentum matrix elements can only be safely interchanged if the Hamiltonian is free of non-local terms and eigenstates are well-represented in the working Hilbert space.
In summary, this work may well serve as a complete as well as rigorous guide to the intricate relations behind the evaluation of the velocity, momentum, and position matrix elements in crystalline solids. In this appendix we prove that projectingv = i[Ĥ,r] in CMR along with the corresponding expression for the position operator, Eq. (15), allows to obtain Eq. (9). Let |mk be a general Bloch basis (orthonormal for simplicity) following a distribution normalization. An identification with the eigenstates basis will be made in the end. We have Now we insert the closure relation between the two operators: We have splitted the full matrix elements into two terms according to the two parts in Eq. (15). First we work out φ 1 |v (1) |φ 2 , We have taken into account that a crystal Hamiltonian is diagonal in the k vector, this is H mk,m k ≡ and applying the chain rule, The first term is zero following the conditions required by Blount [8]. Now we look at the Berry connection term We now find the expression in the eigenstates basis. For clarity we rename m = n, and use H nn (k) = n (k)δ nn Ω, obtaining n (k ) nk|v|n k , where nk|v|n k = δ(k − k ) ∇ k n (k)δ nn + iΩ −1 ω nk,n k A nn (k) , which is precisely Eq. (16).

Appendix B DERIVATION OF EQ. (21)
We start from Eq. (9) for the case k = k , nk|v|n k v = iω nk,n k A nn (k) + ∇ k n (k)δ nn Recall that the Berry connection is defined A nn (k) ≡ i u nk |∇ k u n k Ω . Expanding the periodic part in a Bloch Basis, |nk state is |u nk = α c so we have Now we can introduce the Hamiltonian matrix elements in the first two terms according to the eigenvalue equation, yielding The first term cancels the gradient of the energy band in Eq. (41). In order to write the final form of the expression, we note that A αα (k) ≡ i u αk |∇ k u α k Ω = i∇ k S as presented in the main text. Finally this expression is particularized for Bloch states expanded in a local orbital basis, where |αk v = 1/ √ N R e ik·R |αR v , leading to the Berry connection The expression above is a generalization for that of a Wannier basis, see e.g. Ref. [11,12]. Here, an extra term arises accounting from the nonorthonormal character of atomic states, differently from the Wannier orbitals, which are orthonormal by construction. This is also reflected by the appearance of the overlap matrix in the first line of Eq. (20). In Eq. (44), dipole matrix elements between the basis set are integrated in all space and not in the unit cell, different than in the original definition for the Berry connection. This change is done by passing from cell to lim N →∞ 1 N ∞ −∞ , that is well-defined for a periodic integrand. Finally Eq. (43) can be written which is the formula given in Ref. [7].