Assessing the role of interatomic position matrix elements in tight-binding calculations of optical properties

We study the role of hopping matrix elements of the position operator $\mathbf{\hat{r}}$ in tight-binding calculations of linear and nonlinear optical properties of solids. Our analysis relies on a Wannier-interpolation scheme based on \textit{ab initio} calculations, which automatically includes matrix elements of $\mathbf{\hat{r}}$ between different Wannier orbitals. A common approximation, both in empirical tight-binding and in Wannier-interpolation calculations, is to discard those matrix elements, in which case the optical response only depends on the on-site energies, Hamiltonian hoppings, and orbital centers. We find that interatomic $\mathbf{\hat{r}}$-hopping terms make a sizeable contribution to the shift photocurrent in monolayer BC$_2$N, a covalent acentric crystal. If a minimal basis of $p_z$ orbitals on the carbon atoms is used to model the band-edge response, even the dielectric function becomes strongly dependent on those terms.


Introduction
Empirical tight-binding (TB) is the method of choice for obtaining a simple and intuitive description of the electronic structure of solids [1]. In this method, the basis is defined implicitly through the on-site energies and the Hamiltonian hopping matrix elements, without an explicit real-space representation of the basis orbitals. 1 While that information is already sufficient to evaluate many physical quantities (energy bands, elastic constants, phonon spectra, etc.), the calculation of optical responses requires, in addition, the matrix elements of the position (dipole) operatorr in the TB basis.
In TB calculations of optical properties, it is customary to make the simplest possible approximation for the position matrix; namely, to discard all matrix elements ofr between different Wannier orbitals, which we will refer to as the "hopping matrix elements ofr", or simply "r hoppings." When doing so, the only spatial information that is retained in the model are the orbital centers, where τ m is the center of the mth Wannier orbital in the home cell. This minimal spatial embedding of a TB model already allows to incorporate electromagnetic fields in a gaugeinvariant manner [2]. It is, nonetheless, a rather uncontrolled approximation: symmetryallowed intra-atomic matrix elements such as s|x|p x are discarded along with interatomic matrix elements, all of which can in principle contribute to the optical response. The above approximation has been discussed in the literature [3][4][5], and the general prescription for incorporatingr-hopping terms into both orthogonal [6] and non-orthogonal [7] TB models has been described. However, the impact of those terms on calculations of optical responses has not been thoroughly examined. An important step was taken in Refs. [4,8], which examined the corrections from intra-atomicr hoppings to the linear optical response of a toy model.
In this work, we revisit the problem from an ab initio perspective. We restore the hopping terms that are missing from Eq. (2), and employ a first-principles-based WF method to assess their contribution to the linear and quadratic optical responses of a single layer of the graphitic material BC 2 N [9][10][11].
Our calculations use the Wannier interpolation method [6] and they proceed as follows. After performing an ab initio calculation of the electronic structure, we construct, in a postprocessing step, well-localized WFs spanning the relevant bands. We then take those WFs and use them as an orthogonal TB basis to evaluate the band structure and the optical matrix elements. Since the Wannier orbitals are constructed explicitly, ther hoppings can be tabulated and included, along with the on-site energies, orbital centers, andĤ hoppings, in the calculation of optical responses; by selectively discarding some or all of ther hoppings, we are able to gauge their contributions.
The optical responses analyzed in this work are the dielectric function and the shift photoconductivity. The latter is a quadratic response associated with a shift in the center of mass of an electron as it is optically excited from a valence band to a conduction band in a piezoelectric crystal [12][13][14], and is known to be quite sensitive to the spatial embedding of the TB model [15,16]. (The same happens with the ground-state electric polarization [3], which is associated with the center of mass of valence electrons [17].) Our test system, a single layer of BC 2 N, was chosen for the following reasons. (i) Its structure is noncentrosymmetric, a necessary condition for observing a quadratic response. (ii) As it is a covalent crystal with strong orbital overlap between different sites, interatomic terms can be expected to play a significant role. (iii) Its simple band structure near the fundamental gap allows for a complementary study based on a two-band k · p model constructed from the Wannier Hamiltonian, which further highlights the impact of approximation (2).
The paper is organized as follows. In Sec. 2 we review the expressions for the dielectric function and shift photoconductivity in the independent-particle approximation. The connection between Wannier interpolation and TB theory, and the contribution ofr hoppings to the optical matrix elements, is discussed in Sec. 3. In Sec. 4 we describe the ab initio and Wannier-interpolation calculations, and in Sec. 5 we present and analyze the numerical results. We conclude in Sec. 6 with a summary and discussion, and in Appendix A we describe how the k · p model is constructed.

Dielectric function and shift photoconductivity
In a nonmagnetic material such as BC 2 N, the dielectric function is a symmetric tensor whose imaginary part is absorptive, with an interband contribution given by [14] ab (ω) = Here e > 0 is the elementary charge, f knm = f kn −f km and ω knm = E kn −E km are differences between occupation factors and between band energies, [dk] = d d k/(2π) d in d dimensions, and the integral is over the first Brillouin zone (BZ). Finally, is the interband dipole given by the off-diagonal part of the Berry connection matrix where |u km denotes the cell-periodic part of a Bloch state |ψ km . Noncentrosymmetric crystals display a nonlinear optical effect known as the bulk photovoltaic (or photogalvanic) effect, which can be divided phenomenologically into "linear" and "circular" parts [18][19][20]. The linear bulk photovoltaic effect occurs in piezoelectric crystals such as BC 2 N, and is described by the relation The shift current corresponds to the interband part of this response and is given by [14] σ abc (0; ω, −ω) = πe 3 where denotes the gauge-covariant "generalized derivative" of the interband dipole with respect to k (∂ a = ∇ ka is the ordinary k derivative).

Wannier interpolation 3.1 Energy bands
After the ab initio total-energy calculation, we construct in a post-processing step a set of well-localized WFs. These are chosen to span a group of bands that includes the initial and final states involved in interband absorption processes up to some desired frequency. From the Wannier orbitals we then define a set of Bloch basis states as take matrix elements of the ab initio Hamiltonian between those states, and diagonalize the resulting matrix, With a proper choice of WFs, the eigenvalues {E km } provide a smooth interpolation across the BZ of the selected group of density functional theory (DFT) bands [21].

Optical matrix elements
The same interpolation strategy can be applied to other k-space quantities. In particular, the transformed cell-periodic states interpolate the ab initio cell-periodic eigenstates, allowing to treat wavefunction-derived quantities such as the Berry connection [6]. Inserting Eq. (13) in Eq. (6) yields where and Eq. (3) was used in the last step. We will refer to A k and a k as the "internal" and "external" parts of the Berry connection matrix A k . 2 The convention adopted in Eq. (10), with the phase factor e ik·τm included in the Bloch sum, is the most natural one for dealing with geometric quantities in k space [23]. With that convention, when allr hoppings d nm (R) are discarded Eq. (15) for A (W) k vanishes identically so that A k reduces to the A k term given by [6] A knm = where the numerator is an "effective" TB velocity matrix element obtained by . The A k term is "internal" in that it only depends on the minimal TB ingredients contained in Eqs. (1) and (2): on-site energies, (interatomic)Ĥ hoppings, and Wannier centers. Let us now restore the "external"r-hopping terms, and examine their contributions to A k . Intra-atomicr hoppings make a k-independent contribution to Eq. (15), while the contribution from the interatomic ones varies with k. 3 In the limit where the overlap between orbitals on different sites is negligibly small, all interatomic terms vanish and Eq. (14) reduces to denotes the intra-atomic part of ther-hopping matrix d(R). Away from that limit, interatomicĤ-andr-hopping terms give competing contributions to A k .
Thus, the standard treatment of optical properties in TB can be systematically improved by progressively addingr-hopping terms to the model in a sequence of steps: 1. Only include on-site energies, orbital centers, andĤ hoppings.
Step 1 corresponds to the uncorrected TB model, step 2 adds the intra-atomic corrections, and steps 3 and higher incorporate interatomic corrections. Intra-atomic corrections are expected to be important for optical transitions involving localized d orbitals [4]. The role of interatomic corrections remains less clear, as they were not considered in some of the previous works; they will be included below in our study of monolayer BC 2 N.

Structural and computational details
The structure of monolayer BC 2 N, depicted in Fig. 1, is formed by alternating zigzag chains of carbon and boron nitride atoms. The unit cell contains four atoms; the two carbon atoms are inequivalent, and we use the symbols C N and C B to label the ones with a nitrogen and a boron atom among their nearest neighbors (NNs), respectively. The structure is polar alongŷ, and has mirror symmetries M x and M z . The BZ is also shown in Fig. 1, with the high-symmetry points indicated.
The DFT calculations are performed using the Quantum ESPRESSO code package [24]. We treat the core-valence interaction using scalar-relativistic projector augmented-wave pseudopotentials taken from the Quantum ESPRESSO website. The pseudopotentials are generated with the Perdew-Burke-Ernzerhof exchange-correlation functional [25], and the energy cutoff for the plane-wave basis expansion is set at 70 Ry.
To study the monolayer we employ a slab of length l =20Å along the out-of-plane direction, and take the lattice constants (a 1 = 2.46Å alongx and a 2 = 4.32Å alongŷ) from Ref. [11]. For the self-consistent calculation we use a 10 × 10 × 1 k -point mesh, and switch to a 20 × 20 × 1 mesh for the non-self-consistent step from which we obtain the Bloch functions that are used as input for the Wannierization procedure [21,26].
We use the Wannier90 code package [27] to generate atom-centered WFs spanning the  [21]. Energies are measured from the Fermi level.
states around the Fermi level. To that end, we employ a one-shot projection scheme [26] combined with band disentanglement [21]. We consider two different sets of WFs: one with a p z orbital on every atom (with an average quadratic spread of 1.08Å 2 ), and another with p z orbitals on the carbon atoms only (with an average spread of 3.39Å 2 ). The corresponding Wannier-interpolated bands are plotted in the two panels of Fig. 2 along with the DFT bands.
To obtain well-converged Wannier-interpolated spectra for the dielectric function and shift photoconductivity, we employ a dense 2000 × 2000 × 1 k -point interpolation grid. We use a fixed width of 0.01 eV when broadening the delta functions in Eqs. (4) and (8), and apply a broadening of 0.04 eV to regularize the contribution of intermediate states to the shift current [22]. The occupation factors are evaluated at T = 0 K.

Electronic structure and optical spectra
The DFT band structure of monolayer BC 2 N (black lines in Fig. 2) is replotted in Fig. 3 over a narrower energy range around the Fermi level, together with the direct gap and with the density of states (DOS) projected onto atomic p z orbitals. The minimum gap of ∼ 1.6 eV occurs at the S point, and the electronic states near the band edge are composed almost entirely of p z orbitals, with only a small contribution from other orbitals -mainly p y -at higher energies (not shown). More precisely, the states at the top of the valence (bottom of the conduction) band are primarily composed of p z orbitals on the C B (C N ) atoms, with smaller but non-negligible contributions from p z orbitals on the N and B atoms. These features guided our choice of the two Wannier basis sets described in Sec. 4. In particular, the larger set with one p z orbital on every atom captures the dominant orbital character at the band edge as well as the covalent nature of the system; we will use it in most of our calculations, with the exception of Sec. 5.3.1 where we switch to the minimal basis with p z orbitals on the carbon atoms only. Next we analyze the optical response tensors evaluated by Wannier interpolation using our reference basis set. The symmetry-allowed components of the dielectric function are xx and yy , and those of the shift photoconductivity σ yxx , σ yyy , and σ xxy = σ xyx ; for the sake of clarity, we focus on σ yxx and σ yyy . We report bulklike values for aa and σ abb , by rescaling the values obtained for the slab by a factor of l/h≈ 4.0 (l = 20Å is the slab thickness, and h = a 2 1 + a 2 2 = 4.97Å is the stacking distance) [16]. The joint density of states (JDOS) is shown in Fig. 4(a). It exhibits van Hove singularities at E g ∼ 1.6 eV and E ΓX ∼ 2 eV, and a strong peak at ∼ 2.4 eV. These features carry over to the dielectric function in Fig. 4(b). Between E g and E ΓX , yy is negligibly small; this is caused by mirror-parity selection rules [28] that are exact at E g and hold to a good approximation up to E ΓX . Turning to the shift photoconductivity in Fig. 4(c), σ yxx also peaks at ∼ 2.4 eV while σ yyy remains significantly smaller over the entire range shown. In the band-edge region, σ yyy ≈ 0 due to the same selection rules that enforce yy ≈ 0; in contrast, σ yxx is sizeable and shows a step-like feature followed by a plateau. Overall, the shapes of the −σ yxx and σ yyy spectra are reminiscent of those of xx and yy , respectively.

Four-orbital tight-binding model
Since our Wannier basis contains one orbital per atom, theĤ andr hoppings in Eqs. (1) and (3) are purely interatomic. Their behavior as a function of distance between orbitals is shown in Fig. 5. As expected from the localized nature of WFs,Ĥ hoppings decay very rapidly with distance: the four NN hoppings at ∼ 1.5Å dominate over the rest by more than an order of magnitude, and hoppings beyond 3Å are negligible. The behavior ofx andŷ hoppings is more complex, with 1st NN coefficients being less than half the size of the dominant 2nd NN ones at ∼ 3Å. The largest coefficients overall areŷ hoppings between boron pairs separated along x, but hoppings as distant as ∼ 6Å remain sizeable. The longer range ofr hoppings compared toĤ hoppings, as well as their non-monotonic behavior, seem reasonable given that ther operator grows linearly with distance. (a)  Figure 5: Decay, as a function of distance between Wannier orbitals, of the hopping matrix elements ofĤ,x, andŷ, color-coded by the atomic combination forming the orbital pairs. The Wannier basis consists of one p z orbital on every atom.
Let us now analyze the impact ofr hoppings on the calculated xx and σ yxx spectra. Before coming to the full spectra, we first inspect the associated transition matrix elements between the top valence and bottom conduction bands. Their dispersions are plotted in Fig. 6;    Figure 7: (a) and (b) xx and −σ yxx spectra, for different levels of truncation of ther hopping matrix. The Wannier basis and labelling scheme are the same as in Fig. 6, with the difference that here ther hoppings are included in a cummulative way. The inset zooms in on −σ yxx at the band edge.
in panels (a,b) we compare full Wannier-interpolation results with uncorrected TB results obtained by setting allr hoppings to zero, and in panels (c,d) we show the separate corrections to both matrix element from 1st and 2nd NNr hoppings. The uncorrected TB approximation works very well for the xx matrix element: the largest relative error, which occurs around the band edge S, does not exceed 5%. That approximation is much less satisfactory for the σ yxx matrix element, especially around S where the relative error reaches ∼ 50%. In the case of the xx matrix element, the small corrections to the TB approximation come mostly from 2nd NN x hoppings; this is consistent with the behavior of those hoppings in Fig. 5(b), where 1st NN terms at ∼ 1.5Å are much smaller than 2nd NN ones at ∼ 3Å. The largest corrections to the σ yyx matrix element again come from 2nd NNr hoppings, as expected from Figs. 5(b,c) (theŷ hoppings contribute through r x;y kvc ). The full spectra xx and σ yxx , calculated both with and withoutr-hopping corrections, are displayed in Fig. 7. Consistent with the preceeding analysis, those corrections are fairly minor for xx , but they are significant for σ yxx in the band-edge region, wherer hoppings up to 2nd NN have a sizeable impact on the computed spectrum.

Two-band models for the band edge
With the aim of reproducing the electronic and optical properties at the band edge in the simplest possible way, we now turn to minimal two-band models. We consider below two such models constructed in different ways.

Tight-binding model
The first minimal model we consider is a TB model whose basis consists of p z -like WFs on the carbon atoms, forming a quasi one-dimensional chain along x (see Fig. 1). As discussed in Sec. 4, this model is expected to yield acceptable results only at the band edge, where carbon p z states are prevalent. The decay of theĤ andr hoppings in this model is shown in Fig. 8. Compared with the four-band model (Fig. 5), the decay is significantly slower. The largestĤ hoppings have similar magnitudes in both models, while the largestx andŷ hoppings are an order of magnitude larger in the minimal model; as a result,r hopping corrections are much more pronounced. This can be seen in the dispersions of the xx and σ yxx matrix elements in Fig. 9: around S the corrections are significant already for xx [panels (a,c)], while in the case of σ yxx [panels (b,d)] they reduce the matrix element by almost a factor of three, with the largest corrections coming from 1st and especially 2nd NNr hoppings.

k · p model
As an alternative approach for constructing a minimal band-edge model, we now extract a two-band k · p effective Hamiltonian from our reference four-band TB Hamiltonian. This is motivated in part by previous works [16,[29][30][31], where two-band k · p models were used to describe the band-edge photocurrent response of different types of materials.
Our k · p model is constructed as described in Appendix A, by expanding the TB Hamiltonian (11) to second order in k around the S point, and then applying Löwdin perturbation theory. The result is a transformed 2 × 2 HamiltonianH k , which we expand in terms of the identity matrix and of the Pauli matrices as Its band dispersion in Fig. 3 agrees well with the DFT dispersion near the band edge.  Near the band edge, the dielectric function and shift photoconductivity are given by the product between the transition matrix elements and the JDOS [16], The k · p expressions for the matrix elements at a band extremum read [16,32] |r a and where the coefficients f i , f i,a = ∂ a f i , and f i,ab = ∂ 2 ab f i are evaluated at the S point using Eq. (32) in Appendix A, and ε ijm is the Levi-Civita symbol. Note that these expressions only depend on the k · p Hamiltonian, which in turn is constructed from the TB Hamiltonian; thus,r hoppings are not taken into account when evaluating optical matrix elements from a TB-derived k · p model. The k · p results for the JDOS, xx , and σ yxx are shown as dashed red lines in the shaded regions of Fig. 4. Panel (a) shows a good agreement with the Wannier-interpolated JDOS around the band gap: the height of the step-like feature at E g is nicely reproduced, and although the Wannier-interpolation curve grows monotonically above E g while the k · p one stays flat, the discrepancy is small. In panel (b), the k · p curve for xx matches very well the Wannier-interpolation one over the entire band-edge region: it reproduces not only the step height at E g but also the subsequent decrease, thanks to the 1/ω 2 factor in Eq. (21). In contrast, in panel (c) the k · p curve for σ yxx deviates considerably from the Wannierinterpolation one, overshooting it by about a factor of five at E g . This is in line with our previous finding that the σ yxx matrix element at the S point gets strongly reduced whenr hopping terms are included.

Discussion
The question of how to evaluate optical matrix elements was debated in the TB literature until the early 2000s (see Ref. [4] for an overview). It eventually became clear that the "minimal TB substitution"v → (1/ )∇ k H (W) k for the velocity matrix elements leaves out important physics. In particular, it completely neglects intra-atomic dipole transitions [4], as well as corrections to interatomic transitions from off-site dipole matrix elements. Although the shortcomings of the minimal (or uncorrected) TB approach to the calculation of optical matrix elements are by now well understood, that approach continues to be widely used because it requires no additional parameters beyond the standard ones: on-site energies,Ĥ hoppings, and orbital centers.
The development of Wannier-interpolation schemes by the ab initio electronic-structure community provided an opportunity to assess the importance of those additional TB parameters (hopping matrix elements ofr) to various transport and optical responses in a wide range of materials. In the early works on Wannier interpolation [6,33], the anomalous Hall conductivity and magnetic circular dichroism spectrum of bcc Fe were found to be very well described by the uncorrected TB approach where allr hopping terms are discarded.
More recently, Wannier interpolation has been used to calculate nonlinear optical responses, and more pronouncedr-hopping corrections were found in some cases, such as the shift photoconductivities of WS 2 [34] and especially GaAs [22], and the high-harmonic generation spectrum of monolayer hexagonal BN [35]. This provided the motivation for the present work, where we carried out a systematic study of the impact of hopping matrix elements ofr on the linear (dielectric function) and quadratic (shift photoconductivity) optical responses of monolayer BC 2 N.
Our findings can be summarized as follows: (i)r hoppings decay more slowly with distance thanĤ hoppings, and 2nd NNr hoppings can be significantly larger than 1st NN ones; (ii)r hoppings are more important for the shift current than for the dielectric function, indicating that the former is more sensitive to the spatial structure of the WFs; (iii) the importance of r-hopping corrections increases as the number of Wannier basis orbitals decreases, since the Wannier orbitals tend to be more extended in minimal bases; (iv) two-band k ·p Hamiltonians constructed from TB Hamiltonians are likely to provide a poor description of the shift-current response at the band edge, due to the neglect ofr-hopping corrections.
Despite being available "for free" within the Wannier interpolation framework,r hoppings are often discarded in Wannier-based calculations of nonlinear optical responses (for some recent examples, see Refs. [31,[36][37][38]). Our results indicate that this common practice is ill-advised, since the error incurred can sometimes be significant.
One question that the present work does not address, and which could be an interesting direction for future research, is how to devise reasonable approximations for the hopping matrix elements ofr in the context of empirical TB theory, where the basis orbitals are not explicitly available.