Direct and indirect excitons in boron nitride polymorphs: a story of atomic configuration and electronic correlation

We compute and discuss the electronic band structure and excitonic dispersion of hexagonal boron nitride (hBN) in the single layer configuration and in three bulk polymorphs (usual AA' stacking, Bernal AB, and rhombohedral ABC). We focus on the changes in the electronic band structure and the exciton dispersion induced by the atomic configuration and the electron-hole interaction. Calculations are carried out on the level of \textit{ab initio} many-body perturbation theory (GW and Bethe Salpeter equation) and by means of an appropriate tight-binding model. We confirm the change from direct to indirect electronic gap when going from single layer to bulk systems and we give a detailed account of its origin by comparing the effect of different stacking sequences. We emphasize that the inclusion of the electron-hole interaction is crucial for the correct description of the momentum-dependent dispersion of the excitations. It flattens the exciton dispersion with respect to the one obtained from the dispersion of excitations in the independent-particle picture. In the AB stacking this effect is particularly important as the lowest-lying exciton is predicted to be direct despite the indirect electronic band gap.

All first principles calculations predict that the electronic gap is direct for the single layer [28,29], but indirect in the bulk, and actually even for a bilayer [30,31]. The most advanced band structure calculations have been obtained within the GW approximation which include electron-electron correlations explicitly. However, when discussing two-particle response functions (optical spectroscopy, energy loss spectroscopy, or x-ray scattering), this scheme is often insufficient because it lacks the important electron-hole interactions. This is particularly true in systems where the screening is weak, like in thin films and large gap bulk semiconductors, which is the case for hBN. Indeed, a key element to gain quantitative insight in hBN properties has been the analysis of the excitonic dispersion at finite Q along the M and the K lines [32].
In this paper we push the analysis further, discussing in detail the dispersion of the excitons especially in the K direction, relevant for optical spectroscopy, with the intent of tracking down the structural elements affecting the excitonic properties and the excitonic dispersion in particular. To this aim, the case of hBN monolayer is investigated and compared with three different polymorphs of bulk hBN, reported in Fig. 1. The usual structure is the so-called AA stacking where B and N atoms alternate along the stacking axis. Another stable structure, although less common, is the ABC rhombohedral one [33], and finally the AB Bernal stacking has been reported for few layers [34]. Our analysis combines ab initio Bethe-Salpeter calculations [35] and a tight-binding Wannier model which has already been shown to be fairly accurate for hBN [29]. Our work integrates and completes the discussion of excitonic effects in multilayer hBN at Q = 0 published elsewhere [31], as well as the investigation of optical spectra in different stacking sequences [11,26].
The paper is organized as follows. In Sec. II we detail the computational parameters used for the ab initio calculations and we devise the tight-binding (TB) model used. In Sec. III we focus on the monolayer calculation. By comparing the ab initio and the TB results, we are able to validate the approach and to appreciate the reliability of the TB predictions in this system. Moreover, the monolayer hosts the fundamental inplane physics and constitutes the building block of the three bulk polymorphs. In Sec. IV, the comparison of the different bulk phases is carried out with a highlight on the way the stacking sequence affects the reference in-plane electronic FIG. 1. The atomic structure of hBN monolayer, the AA , the AB, and the ABC stacking. In the latter, the unitary lattice vectors of the rhombohedral cell are also reported. A red line joining the centers of three B 3 N 3 hexagons highlights the difference between the three bulk structures. and excitonic properties. Finally, in Sec. V we draw our conclusions.

II. THEORETICAL METHODS
In this section we report the parameters used in the ab initio calculations and introduce the tight-binding model developed for the calculation of one-particle and two-particle properties of the four systems considered.

Free-standing single layer
For the boron nitride single layer we use an in-plane lattice parameter of a = 2.50 Å and a distance between the periodic replica of 13 Å. The Kohn-Sham states and energies, the GW corrections, and the excitonic properties have been computed with the GPAW code [36]. Projector-augmented wave (PAW) methods have been used for both atomic species. DFT energies and wave functions have been obtained within the PBE exchange-correlation potential, using a plane-wave cutoff energy of 40 Ha. We did not find relevant differences in the exciton dispersion with respect to the local density approximation (LDA) results. For the density calculation we used a 12 × 12 × 1 -centered grid. The Bethe-Salpeter equation (BSE) has been solved using a truncated Coulomb potential on a 36 × 36 × 1 -centered k-point grid, with a cutoff energy of 100 eV and including 60 bands in the calculation of the dielectric constant. Three valence and three conduction bands have been included in the excitonic Hamiltonian and quasiparticle energies have been approximated with a scissor operator of 2.75 eV adjusted on a recently published GW result [29].

Bulk structures
For the three bulk structures we used the same in-plane parameter a as for the monolayer and an interlayer distance c = 3.25 Å. The value of c is in agreement with measures of the AA stacking published in a previous work of ours [37]. The Kohn-Sham equations and the GW corrections have been calculated with the plane-wave simulation package ABINIT [38]. Norm-conserving Troullier-Martins pseudopotentials have been used for both atomic species. DFT energies and wave functions have been obtained within the local density approximation (LDA) to the exchange-correlation potential, using a plane-wave cutoff energy of 30 Ha for the three stackings. The hexagonal Brillouin zone of the AA and the AB stackings have been sampled with a 8 × 8 × 4 k-point grid, while the Brillouin zone of the rhombohedral cell of the ABC stacking has been sampled with a 9 × 9 × 9 grid. All k-point grids are centered.
Quasiparticle corrections have been obtained within the perturbative G 0 W 0 approach. They have been computed on all points of a 6 × 6 × 4 -centered grid for AA and AB stacking and on a 9 × 9 × 9 -centered grid for ABC. A cutoff energy of 30 Ha defines the matrix dimension and the basis of wave functions for the calculation of the exchange part of the self-energy. The correlation part has been computed including 600 and 150 bands in hexagonal and rhombohedral structures, respectively, and using the same wave function basis as for the exchange part. To model the dielectric function, the contour deformation method has been used for AA , computing the dielectric function up to 60 eV, but this showed to give negligible improvements with respect to the Godby-Needs plasmon pole approximation, so the latter has been used in the other structures. To obtain the GW energies along high-symmetry lines and on finer grids, the GW corrections have been interpolated. Note that the Brillouin zone of the ABC rhombohedral stacking is larger than the hexagonal cell of AA and AB systems, so the band structure has been folded to the hexagonal cell for a consistent comparison. For a given point on the hexagonal cell, three points in the rhombohedral cell have been taken: k and k ±k withk = (1/3, 1/3, 1/3) expressed in the reciprocal coordinates of the rhombohedral cell.
The macroscopic dielectric function M (q, ω) has been calculated on the GW-BSE level using the EXC code [39]. For the hexagonal structures AA and AB, we included six valence bands and three conduction bands, fixing a cutoff energy of 360 eV for both the matrix dimension and the wave function basis. The static dielectric matrix entering the BSE kernel has been computed within the random phase approximation with local fields, including 350 bands and with cutoff energies of 120 eV and 200 eV for the matrix dimension and the wave function basis, respectively. In the case of the ABC structure, the BSE has been solved including four valence and four conduction bands and keeping the same cutoff as for the hexagonal structures. The static screening used has the same parameters as for the underlying GW calculation. With these parameters, the energies of the excitons are converged within 0.05 eV in all structures. The dispersion of the exciton as a function of q in the ABC cell required the same folding procedure explained above.
The two triangular sublattices B (green) and N (gray) forming the hBN honeycomb lattice. Red arrows mark the 1nn vectors τ α and the positions R, R of the hopping electron when the hole is placed on an N site at the origin (blue circle). (b) Scheme of the four hopping terms t ⊥ , t 2⊥ , t , and t 2 in the AA stacking.

Free-standing single layer
We first consider the TB model of the free-standing monolayer, already introduced by some of us [29]. Let us recall its main features. The real space honeycomb lattice can be divided into two triangular sublattices B and N , either connecting all sites of B or N atoms, respectively. Vectors τ {1,2,3} connect neighboring sites of the two sublattices. A scheme of the structure is reported in Fig. 2(a). From the energetic point of view, we note ± the on-site energies on B (+) and N (−). First-and second-nearest neighbor (from now on contracted in 1nn and 2nn) hopping integrals are t ⊥ and t 2⊥ , respectively. The latter is assumed to be equal for B − B and N − N hopping. With these ingredients and by passing to the Bloch representation as in Ref. [29], the TB energies in the single layer can be approximated for k ⊥ along the KM direction as: where γ (k ⊥ ) = α=1,2,3 e ik ⊥ τ α and k ⊥ is strictly in-plane. In the above expression the (+) sign is for conduction e states localized on B sites, and (−) is for valence h states localized on N . Note that without the 2nn contribution valence and conduction bands would be symmetric. The description of the exciton relies on a Wannier TB model [29], but in this work we extend it beyond the optical limit and add improvements on the electron-hole interaction. Let us first recall the basics of the model in the optical limit Q = 0. The fact that we can describe the concerned oneelectron π bands with atomiclike Wannier functions allows us to work directly with excitonic Wannier equations in real space [35,[41][42][43][44]. Assuming the hole to be fixed on a N site, and using relative coordinates R for the electron-hole pair, we have reduced the problem to the one of the electron hopping on sites of B in the presence of an attractive impurity located at the origin, as represented schematically in Fig. 2(a). The corresponding Bethe-Salpeter-like Hamiltonian H eh = H 0 eh + + V contains a kinetic part H 0 eh , a screened Coulomb term that can be taken as a fitting parameter and an exchange term V ∝ Q that actually vanishes in the Q = 0 limit.
We extend now the model to Q = 0 hence allowing the exciton to move across the layer. Through the definition of appropriate Bloch states |R, Q defined in Appendix B, the momentum Q is a good quantum number related to the propagation of the center of mass of the electron-hole pair. We can then adopt a mixed representation, where the motion of the pair is treated in reciprocal space Q and the relative electron-hole distance in real space R. For the monolayer, it can be demonstrated that the simple model above can be extended to a Hamiltonian H eh = H 0 eh + + V where the kinetic part reads: τ being |τ α | for any α = 1, 2, 3. The Coulomb part reads with R an appropriate attractive potential. In 2D, an adjustable Keldysh potential seems to be the most pertinent choice [29]. Since the effective hopping integrals now depend on the direction of the hopping step, the symmetry of the problem is much lower, but calculations can easily be done with the same techniques as for Q = 0. See Galvani et al. [29] and Appendix B for the details of the derivation. At variance with the optical limit, at finite Q the exchange term V shall not be neglected. In fact it is responsible for dipolelike coupling between different sites which induces a singularity in the exciton dispersion at Q → 0 [43]. In our TB model it produces effective interactions R , Q|V |R, Q when the electron and the hole are sufficiently close (within the 1nn shell). In the mixed representation these interactions are therefore local with respect to the electron-hole relative distance R, but have short-and long-range components with respect to the propagation of the center of mass Q. The short-range components produce analytic terms that shift upward the dispersion curve with respect to states where V = 0, such as triplet states. Instead, the long range contribution is a dipoledipole term whose Fourier transform is linear for small Q and singular at Q = 0, as expected. As shown in Appendix C, an approximate expression valid for the single layer reads: where J will be considered here as a parameter to be fitted to ab initio calculations, and where Q = |Q| and τ = |τ α | for any α.

Bulk structures
When passing to bulk systems, some aspects complicate the model. The first is that the corresponding Bloch states are characterized by k = k ⊥ +ẑ k . Moreover in the AA and AB stacking, the basis is formed of four Bloch states because of the four atoms in the unitary cell. This problem is not encountered in the ABC stacking as long as one works in the elementary rhombohedral cell with only two atoms. The second complication is connected to the anisotropy of the layered structures that requires us to distinguish between in-plane and interplane screening. Finally, interlayer coupling is accounted for by a hopping term t linking two neighboring sites that are vertically aligned and a parallel 2nn hopping t 2 between second-nearest neighbors of different planes, again assumed equal for B − B and N − N hopping. A sketch of the hopping terms in the AA case is given in Fig. 2 Once properly generalized and retaining only 1nn intraand interplane couplings, the energy of the π * conduction states reads: where we recall that c is the interlayer distance. For valence states E σ kh = −E σ ke with σ = AA , AB, or ABC. In the case of AA and AB stacking the total number of bands is four (two e and two h bands), while the ABC stacking has only two bands consistently with the number of atoms in the respective unitary cells. However, in order to report the ABC dispersion in the smaller hexagonal cell, bands have to be folded by substituting k c in Eq. (6) with k c + 2πλ/3 with λ = 0, ±1 which results in a total of six bands. The corresponding 2nn expressions are unnecessarily complicated to be reported here, but exact expressions including all terms are reported in Appendix A. In the successive calculations, t 2 and t 2⊥ have been included. While t 2⊥ breaks the E h = −E e symmetry in the same way as in the monolayer, the t 2 term breaks the valence-conduction symmetry of the splitting between bands of the same character (valence or conduction). These aspects will be discussed more in detail in Sec. IV.
The generalization of the Bethe-Salpeter-Wannier equation passes through the following steps.
(i) Since we will not consider vertical dispersion, the momentum Q is still in-plane as in the single layer.
(ii) The most delicate part is the indexing of the relative distances between the electron and the hole since now they can localize independently on different planes. In the AA stacking the B and N sublattices are swapped from one layer to the other and are shifted in the AB and ABC cases. As a consequence the relative distance between the two particles belongs to a lattice that depends on the localization of the two particles. We keep track of this information following the approach introduced recently in few-layer hBN [31] consisting in (i) fixing the origin on a N site of one plane, (ii) introducing shifted triangular sublattices + u with u a proper translation vector [40] (including in-plane and vertical shiftsẑc) and (iii) generalizing R to R α,β where α and β label the sublattice occupied, respectively, by the electron and the hole [45]. The example of the AA stacking is given in Fig. 3.
(iii) The matrix elements of the kinetic term (2) are expressed in terms of and in-plane and out-of-plane hopping terms according to the exciton state |R α,β , Q .
(iv) In the Coulomb matrix elements (3) we use a standard 3D R = 1/( R) potential instead of the Keldysh potential. We took into account the anisotropy of the material treating and ⊥ as distinct fitting parameters.
(v) In the Q = 0 limit in 3D, the Q = 0 singularity of the Coulomb potential is stronger than in 2D since at low Q it varies as 1/Q 2 . This induces discontinuities at Q = 0 leading to an upward shift of longitudinal modes that in bulk hBN is about 1 eV [43,46]. Since the neglect of this effect does not change the conclusions regarding the nature of the gapwidth, we neglected the exchange term (4).

Single-particle band structure
In Fig. 4(a) we report both the DFT and the quasiparticle results for the hBN monolayer band structure. The quasiparticle gap, equal to 7.25 eV [29], is direct between the π and π * bands at point K in the Brillouin zone, while the bands are very flat along the KM lines. These valence and conduction bands are not completely symmetric, the dispersion being larger for the valence band. As discussed elsewhere, these results are in agreement with several previous calculations [21,[28][29][30][47][48][49][50].
The TB model introduced previously gives the best fit to the KM region of the ab initio band structure with = 3.625 eV, t ⊥ = −2.33 eV, and t 2⊥ = −0.4 eV. The fit to ab initio data and the results presented here have not been obtained from the approximate expression (1), but from the exact diagonalization of the full TB Hamiltonian reported in Appendix A.
The inclusion of t 2⊥ breaks the electron-hole symmetry, as clearly shown in equation (1), by reducing the effective hopping integral of the conduction band and increasing that of the valence band. As a consequence the conduction band is flatter than the valence band, in agreement with the ab initio results. The sensitivity of the band structure to the 2nn contribution is exemplified in Fig. 4(b) where different values of t 2⊥ have been used. Note that the (direct) gap at K is not modified because the symmetry of the crystal leads to γ (k ⊥ = K ) = 0. To conclude, this simple TB model is able to reproduce the π states in the regions where excitons are relevant for the optical properties.

Exciton dispersion
The dispersion curves E exc ( Q) obtained from ab initio calculations are shown in Fig. 5(a). When Q = we recognize the excitons already characterized in previous works, with the doubly degenerate ground state exciton of symmetry E (or 1s in the atomiclike notation) [29,46]. At higher energy we have the dispersion of the six 2s and 2p states. Note that in either direction one can recognize some additional parabolic bands at high energy close to the zone boundary. We will come back to this characteristic later, when discussing the TB model. In both M and K directions, the second exciton is much more dispersing than the ground-state one. Actually it has been shown that its linear dispersion at Q → 0 is a peculiarity of the 2D geometry which is generated by the exchange contribution to the electron-hole interaction kernel in the Bethe-Salpeter equation [50,51]. Finally we highlight the weak dispersion of the first exciton along K and in particular the fact that the energy at Q = and Q = K basically coincide. This is expected because the most intense single-particle transitions (and thus the most important contributions in the excitonic spectrum) at Q = and Q = K come from vertical K → K and slant K → K transitions, respectively. Since the single-particle states K and K have the same energy, the resulting dispersion of the exciton along K is expected to attain the same value at the extrema of the path. Our results along M are very similar to those published in literature [46,50].
The exciton dispersion has been also computed within the TB approach and is reported in Fig. 5(b). The J term of the exchange contribution (4) has been treated as a fitting parameter, fixed here at 5 eV. The curves have been obtained by diagonalizing a matrix involving 860 sites with a Keldysh potential ranging up to the ninth shell. At moderate values of Q the TB model agrees very well with ab initio results. In particular it reproduces the double degeneracy of the ground state at Q = 0. The two corresponding states | + and | − can be taken as two "circular" states whose components on the three τ 1,2,3 sites τ | + and τ | − are proportional to the cubic roots of unity (1, ω, ω 2 ) and (1, ω 2 , ω), respectively.
As soon as one moves away from Q = 0 the degeneracy is lifted according to their different dipolar orientation. The lower-energy exciton has a transverse orientation at low Q, i.e., the electron-hole dipole is perpendicular to the momentum Q. This makes it optically active and its dispersion is insensitive to the exchange term. Instead, the higher energy exciton is longitudinal, its dipole being parallel to Q. This makes it optically dark and particularly sensitive to the exchange term which has a linear dependence at low Q. Indeed is mildly lifted by the Coulomb term alone (red dashed line), whereas the inclusion of the exchange term reproduces the correct linear dispersion. In Appendix D we report an analytical result predicting this behavior within a perturbative treatment of the exchange interaction [51,52].
At large Q, the agreement is less satisfying since the formula used is no longer sufficient (terms involving sums over reciprocal lattice vectors are neglected) and also because only π states are considered. This consideration allows us to point out the origin of the parabolic bands observed in the ab initio calculation close to the zone boundary. These are present also in the TB model along K, so we can ascribe these bands to π → π * excitations. Instead, along M these bands are predicted only in the ab initio solution, so we can advance the hypothesis that they are of the σ → π * or π → σ * type as the σ states are absent in our TB model.
Further improvements on the TB calculation are obtained when 2nn hopping integrals are included, as a consequence of the fact that they break the valence-conduction symmetry of the band structure. When t 2⊥ = −0.4 eV, the main effect is to decrease the dispersion of the lowest exciton along K from 0.5 eV to 0.3 eV, close to the ab initio value of about 0.27 eV, as reported in Fig. 5(b). In view of the simplicity of the TB model the overall agreement is very good.
To go beyond, we can represent the excitonic wave function as a function of R, which is the distance between the electron and the origin where the hole is fixed. Indeed, in the mixed (R, Q) representation and for fixed Q, the excitonic wave function can be expanded in the corresponding space: As discussed in Appendix B, the basis |R, Q and therefore the coefficients R Q are not uniquely defined. With our standard definition (k e , k h ) = (k, k − Q), R Q is in general a complex quantity. In the bottom part of Fig. 5 we show a map of the intensities | R Q | 2 (electronic densities) which are gauge invariant. Here we consider more particularly the first two excitons (degenerate at Q = 0) along M and K. Consider first the M direction (left panels). Since this is a mirror of the point group we expect the wave function to be either odd or even with respect to the reflection symmetry so that the dipole d is either perpendicular to Q (transverse mode, exciton 1) or parallel to it (longitudinal mode, exciton 2). In particular the intensity corresponding to the transverse mode should vanish on the symmetry axis. As shown in the plots, this is clearly the case for exciton 1 starting from the degenerate ground state exciton at Q = 0. Actually both excitons remain fairly localized, but they deform significantly as a function of Q with a tendency to become elongated in a direction normal to M. In the K case (right panels), we no longer expect definite symmetries except at and K points, but we can see in the figure that the first exciton remains fairly localized with a compact shape except at the middle of K where it tends to elongate. At low Q the system can be considered as quasi-isotropic, but still we have longitudinal and transverse modes.
Let us finally mention an interesting limit for the ground state exciton in the extremely localized case, where the excitonic wave function only extends to the 1nn of the (fixed) hole. It turns out that this model can be completely solved in real space. The resulting lowest mode does not disperse at all in the whole Brillouin zone. This type of flat band has attracted recently great interest in various fields of solid state physics [53]. The corresponding toy excitonic model is described in Appendix E.

IV. THE AA , AB, AND ABC STACKING
It is an established theoretical result that the electronic gap changes from direct at K in the monolayer, to indirect in the AA stacking [22,23,25,32,54]. This effect has important consequences in the optical properties of bulk hBN, as it has been stressed by some recent works [13,14,27,32,55].
Motivated by this, we investigate the impact of the stacking sequence on the single-particle (band structure, electronic gap) and two-particle excitations (exciton dispersion, optical gap).

A. Single-particle band structure
The GW full band structure of the three stackings is reported in Fig. 6 along high-symmetry lines of the hexagonal Brillouin (for the band folding of the rhombohedral ABC structure, see explanation in Sec. II A 2) while the relevant quasiparticle dispersion along the KMK path is plotted in the bottom panels of Fig. 6 where its principal characteristics (dispersion along M and splitting at M) are highlighted in red. Results are also summarized in Table I. The first interesting feature is the behavior of the highest valence and the lowest conduction bands at the K point. In the monolayer, the valence and conduction bands have their extrema at K, but this is not the case in the AA stacking. Here two valence bands cross each other yielding two local maxima in two points close to K. We indicate with letter T the one along K. In the conduction region, two bands also cross at K but they form only one local minimum at T . In the ABC stacking, the K point in the rhombohedral Brillouin zone has no particular symmetry [56], and even if the folded bands show extrema there, they do not correspond to global extrema of conduction and valence bands.
Qualitatively very different is the AB stacking, where the crossing is avoided in both valence and conduction bands because of symmetry reasons, leading to a pretty flat dispersion of the top valence in the vicinity of K. In the conduction band, the same avoided crossing yields a clear local minimum at K. The splitting at K between the two highest occupied states is about 0.2 eV, and it is about 0.4 eV between the two lowest empty states.
A second interesting aspect is the peculiar dispersion of the lowest conduction band along the KM direction, which is the most relevant direction for the optical properties of this material [27,32,55]. The qualitative behaviors of the AA and the ABC phases are similar: Away from K the lowest conduction band disperses almost linearly and has a minimum at M, while in the case of the AB stacking the dispersion is flatter and has a concave shape away from K. Still, beyond a local maximum between K and M it also attains its minimum at M. The dispersion of the bottom conduction is 0.75 eV in the AA , 0.51 eV in the ABC and only 0.10 eV in the AB. It is also worth reporting the energy splitting between the two lowest conduction bands at M: This is 1.68 eV in the AA stacking, 1.11 eV in the ABC stacking, and 1.34 eV in the AB stacking. These data are reported also in Table I.
The nature of the gap also merits to be discussed. The smallest direct and indirect gaps extracted from the band structure are reported in Table II for the GW and the LDA band structure. One immediately sees that the three structures have similar gaps at the LDA level, and in particular there is negligible difference between the AA and the ABC structures. Instead, after the inclusion of GW corrections, the gaps (direct and indirect) of the ABC stacking are sensibly smaller than those of the other two structures, mostly because of the different quasiparticle corrections to the valence bands. This result suggests that there are significant differences in the screening properties of the ABC stacking with respect to those of the other two phases, but a detailed analysis in this respect goes beyond the scope of this paper. Finally, note that in the AA phase, the smallest direct gap is at M. However in T , where optical matrix elements are stronger, it is 4.64 eV and 6.45 eV in LDA and GW, respectively. In the AB stacking it is located at K and in the ABC stacking actually does not lie on a high-symmetry line of the hexagonal cell. In the latter case, the smallest direct gap is actually 5.75 eV, so very close anyway to the direct gap at K reported in Table II. Consider now the AA stacking treated within TB. In the simplest approximation we keep the same and t ⊥ as in the monolayer [57]. and we add only the first-neighbor interlayer hopping t . Then the two conduction π * eigenvalues can be approximated by Eq. (5) (and with opposite sign the two valence π states). From the equation we get that the splitting between the π states vanish either when γ (k ⊥ ) = 0 (H K line, not shown here) or when k c = ±π/2, i.e., on the upper and lower faces of the Brillouin zone (AH L line). In fact since the periodicity along z is 2c, then k c ∈ [− π 2 , π 2 ]. Actually the doubly degenerate state dispersing along H L has exactly the average energy of the two splitted branches on KM. The model reproduces the bottom conduction at M, where the splitting between the two conduction bands is 4t ⊥ t / (|γ (M )| = 1). The fitting procedure to the ab initio band structure yields t 0.5 eV, i.e., t /t ⊥ 0.2 which indicates that the interlayer coupling is actually fairly important.
Sticking to the first-neighbor level, we can already predict from Eqs. (6) and (7) that the conduction-band splitting at M will decrease along the series AA , AB, ABC, and similarly for the valence-band splitting. Indeed in the AB stacking the splitting 2t √ t 2 + t 2 ⊥ / is smaller than in the AA configuration because the number of heteroatomic pairs along the stacking axis is smaller on average. Even smaller it is in the ABC stacking, where the splitting 3 2 t ⊥ t / is indeed the lowest [58]. Moreover in the case of the AB stacking we can verify from Eq. (7) that the bands do not cross at K.
However, to reproduce the electron-hole asymmetry between the conduction and valence bands, and hence the presence of an indirect gap, we need at least second-neighbor interactions within the planes (t 2⊥ ) and between the planes (t 2 ), as discussed in Sec. III 1. The latter term accounts for the difference of splitting, larger in the conduction band than  in the valence band by about 0.4 eV in the AA phase. Approximated formulas are reported in Appendix A together with a summary of the parameters used. In Fig. 7 we report a study of the changes induced in the AA band structure by variations of the 2nn hopping terms. Typical orders of magnitude are t 2⊥ −0.4 eV and t 2 −0.1 eV. Using similar parameters the band structures of the other stackings are also well reproduced within the TB approximation.

B. Exciton dispersion
Let us now pass to the discussion of how the stacking sequence, and hence the changes in the band structure, affect the exciton dispersion. Preliminary results for the AA phase can be found in some previous works of ours [27,32]. In Fig. 8 we report ab initio calculations of the exciton dispersion (black curves) and the free-carrier dispersion (red curves) in the three bulk phases. Quantities related to these dispersion relations and the exciton binding energy in the three systems are also reported in Table III.

Ab initio calculations
First, let us focus on the independent-particle dispersion, or the free-carrier dispersion (GW-IP curves). In the AA case, one recognizes in the convex shape with a minimum at the middle of the K distance the dispersion of the bottom conduction along KM, with a peculiar double-dip shape reminiscent of the conduction band at T . All along the GW-IP dispersion, the energy of the transition passes from 6.28 eV at Q = , corresponding to the smallest direct gap at M, to a minimum of 5.80 eV at 1 Å −1 , corresponding to the indirect gap T M, as already reported in Table II. The resulting dispersion is 0.48 eV. Qualitatively, the same shape characterizes also the dispersion of the first exciton (GW-BSE), but the electron-hole kernel of the Bethe-Salpeter equation has the effect of enhancing the localization of the electron close to the hole and hence of quenching the exciton dispersion to only  exciton is smaller than what an analysis based uniquely on the band structure would suggest. As pointed out by some of us [27], this has an important implication in the difference between photoluminescence and absorption spectra, the former being more sensitive to the lowest (possibly indirect) excitation while the latter displays highest intensity for direct excitations.
Note that the dispersion of the first two excitons cross at the middle of the K path where they have very close energies. As the size of the dots suggests, one of them is bright (large dots) and the other dark (tiny dots).
In order to visualize the distortion of the lowest-energy exciton along its dispersion curve, we plotted the ab initio electronic part of the exciton density | Q (r e , r h )| 2 as a function of r e having fixed r h on an arbitrary N atom. The electronic densities at Q = , Q = K/2, Q = K, and also Q = M are reported in Fig. 9. The data have been obtained by solving the GW-BSE calculation with the YAMBO code [59] on appropriate supercells such that the desired Q point is folded onto . More details on this method to investigate dispersions at finite Q will be available in a future work. It is worth stressing the similarity of these plots with the corresponding TB plots of the first excitons in the monolayer (Fig. 5). This is a consequence of the fact that even in the bulk the first exciton is basically in-plane, so very similar to the first exciton of the monolayer.
In the other panels of Fig. 8, we report the corresponding free-carrier and excitonic dispersions (GW-IP and GW-BSE, respectively) in the ABC and the AB stacking. In either system, one recognizes the shape of the bottom conduction band in the GW-IP dispersion. In the ABC stacking the double dip is not observed, consistently with the absence of a local minimum of the conduction band at K. We remember also that the lowest direct gap of 5.75 eV is not in the high-symmetry lines of the hexagonal cell (so it does not coincide with the entry of Table II). The dispersion in the ABC structure passes from 0.38 eV, in the GW-IP case, to 0.14 eV in the GW-BSE case, reproducing the same flattening observed in the AA case. In this respect, the optical properties of the ABC stacking are predicted to be quite similar to those observed in the AA case except for the smaller band gap.
More interesting is the prediction of the exciton dispersion in the AB phase. In this material the dispersion of the conduction band is predicted to be very weak (see Table I and Fig. 6), so the electron-hole band flattening is expected to have a stronger impact on the exciton dispersion. This is indeed the case: The exciton dispersion is reduced from 0.12 eV in the free-carrier picture to −0.02 eV in the exciton picture, which implies an inversion of the nature of the dispersion. While the indirect nature of the band gap is reflected in the free-carrier dispersion, at the GW-BSE level the lowest exciton is at Q = 0, corresponding to a direct exciton, observable in optics. This implies that in the AB stacking the peaks of luminescence and absorption spectra are expected to coincide, at variance with the other two stackings, and one can expect a stronger luminescence in AB due to the direct transition not mediated by phonons.
Before going deeper, we refer to two appendices where the robustness of these results has been checked. In Appendix F these results have been checked against variations of the interplane distance c, while in Appendix G different approximations have been used at the DFT and the quasiparticle level.

Tight-binding model
We use now the TB model to analyze in detail the influence of the different parameters governing the interplane couplings. We begin with the AA stacking. To take into account the effective anisotropy of the interactions, we use Coulomb potentials 1/( j R) with effective dielectric constants ⊥ = 6 within the planes and = 4.5 between the planes and cutoff radii of 5 Å and 4.5 Å, respectively. These values have been adjusted so as to reproduce reasonably well the first ab initio excitons. In particular the order of magnitude of the binding energy of the lowest exciton is now about 0.5-0.6 eV instead of nearly 2 eV for the single layer [see Fig. 10(a)].
We discuss the role of the different hopping terms by introducing them gradually from the monolayer picture. In Fig. 10(a) we keep only the 1nn hopping t ⊥ term, equal to −2.33 eV. The corresponding exciton dispersion is that of the monolayer, but each curve is doubly degenerate because there are two single layers per unit cell. This also means that (without including more distant neighbor hoppings) at Q = 0 the first exciton would be fourfold degenerate. We should stress here that the dispersion reported in Fig. 10(a) differs from that of Fig. 5 because of the different potential used in the two calculations: screened 3D in the former, Keldysh in the latter. In Fig. 10(b) the interaction between planes is switched on, i.e., t = 0.5 eV, and the layers start being coupled. As a result the degeneracy of either of the two excitons is lifted along the entire K line and at Q = 0 the exciton splits into a couple of doubly degenerate excitons (Davydov splitting). The electron-hole symmetry must be broken in order to have an indirect electronic gap, so we switch on t 2⊥ = −0.4 in Fig. 10(c). As expected, the influence of this parameter is very strong; in particular the lowest exciton branch becomes nearly flat. In Fig. 10(d) a further improvement is obtained by accounting for interplane 2nn integrals t 2 = −0.1 eV which, in agreement with the ab initio results, make appear the indirect minimum of the lowest exciton by breaking the electron-hole symmetry of the TB model.
Let us now compare the different stackings, using the same parameters as for the AA and including all hopping terms. The resulting exciton dispersions are reported in Fig. 11. It can be seen that both AB and ABC stackings show a "direct" lowest exciton. This is not so surprising since the interplane coupling geometries are fairly similar in these two cases and quite different from that of AA (cfr. Fig. 1). We recover the conclusion that both the interplane coupling and the electronhole asymmetry play crucial roles to account for the directindirect excitonic transition.
The general agreement with ab initio data is still pretty good. Indeed, differences are in the range of a few tens of meV, which is remarkable considering the complexity of the BSE calculations, but still not sufficient for reliable description of optical measurements. It is worth pointing out that the quality of the agreement is much poorer with respect to the single-layer case. We believe that the reason for this lies in the different treatments of electronic screening, and in particular in its momentum dependence. In ab initio calculations this is included through a Q-dependent RPA dielectric function and similarly in the TB calculation of the monolayer the Keldysh potential is indeed Q dependent. Instead, in the TB model of the bulk we used a much rougher approximation, distinguishing parallel from perpendicular screening by means of Q-independent constants. Further studies to improve the bulk TB model are currently undertaken in this direction.

V. CONCLUSIONS
We provide a thorough study of the properties of singleparticle and two-particle excitations in hBN monolayer and in three bulk polymorphs: the AA , the ABC, and the AB stackings. We report the first ab initio calculations of the exciton dispersion in the AB and the ABC stacking. Moreover we devise a tight-binding model for the characterization and the analysis of the excitonic dispersion and wave functions in this material. Using these two theoretical approaches we highlight the impact of interlayer interactions on excitonic properties of hBN.
In the monolayer, our ab initio calculations of the band structure and the exciton dispersion along M are in agreement with previously published data [29,46,50]. We also provide the exciton dispersion along K, relevant to discuss the optical properties of the single layer. We have also derived a tight-binding model for the propagating electron-hole pair which includes kinetic, Coulomb, and exchange terms in a Bethe-Salpeter-like formalism. With this model we have highlighted the importance of second-nearest-neighbor hopping terms to describe the dispersion of valence and conduction states and we have been able to analyze the symmetry of the electron-hole pair and to image its wave function.
Concerning the bulk, beside summarizing some recent results on the AA stacking [27,32], we produce some predictions about the ABC and the AB stacking. In the ABC phase, the exciton properties are predicted to be similar to the more common AA phase, in particular in relation with the exciton dispersion. An intriguing result is the prediction of a band gap around 0.5 eV smaller than the AA one. There are indications that this has to be ascribed to specific screening properties of the ABC stacking distinct from those of the other two bulk phases, but we foster for more investigations on this subject.
We generalized our tight-binding model to the three bulk systems and we pointed out the importance of the secondnearest-neighbor in-plane and interplane hopping terms. As their introduction breaks the electron-hole symmetry, they are essential to reproduce the indirect gap and the dispersion of the lowest energy excitons. Also, we indicate a route to improve this model by including appropriate Q-dependent effective dielectric functions.
Finally we predict a peculiar exciton dispersion in the AB stacking. At the single-particle level (band structure), the material exhibits an indirect gap of ∼6.1 eV (probably underestimated of about 0.5 eV [32]), but at the two-particle level (exciton dispersion) the material is predicted to have a direct gap. This is due to the strong momentum dependence of the exciton binding energy E b (Q), which is approximately halved when passing from Q = to Q = 1 Å −1 . This reduction is observed also in the other two materials, but in the AB stacking the variation of E b (Q) is larger than the dispersion of the band structure resulting in a direct exciton transition despite the indirect band gap. This finding will have strong implication for luminescence, that we expect to be much stronger in the AB than in the other bulk phases.

APPENDIX A: TIGHT-BINDING BAND STRUCTURE FOR SINGLE LAYER AND BULK HBN SYSTEMS
In this section we report the exact and approximated formula for the systems studied, together with a summary of the parameters used. In the monolayer the exact diagonalization of TB Hamiltonian leads to expression where the signs (+) and (−) are for conduction and valence states, respectively. The approximate expression (1) valid for t 2 ⊥ |γ (k ⊥ )| 2 holds very well in the KM line. The optimization of the parameters = 3.625 eV, t ⊥ = −2.33 eV, and t 2⊥ = −0.4 eV have been done manually by comparing these exact expressions with the ab initio bands in the KM region. All the analysis on the single-particle band structure and Fig. 4(b) come from the exact expressions above. Instead, in deriving the excitonic model, it has been necessary to adopt the approximate expressions reported in the main text in order to truncate the Wannier functions to the first-nearest neighbors (1nn). Given the range of validity of the approximate formula and the agreement shown with the ab initio results, this approximation is widely justified.
In the AA stacking, the exact expressions including second-nearest-neighbor (2nn) hopping terms are given below for electron states: and hole states: The approximate equation (5) reported in the main text includes only 1nn terms. While , t ⊥ and t 2⊥ have been kept as in the monolayer, the values t = 0.5 eV and t 2 = −0.1 eV have been fitted comparing these exact expressions with ab initio results. In the same way, Fig. 7 has been obtained from these exact expressions. Table IV collects the values of the optimized parameters used in our calculations.

APPENDIX B: TIGHT-BINDING EXCITON MODEL AT Q = 0 (WITH NO EXCHANGE CONTRIBUTIONS)
A general exciton state | Q is obtained by combining electron and hole one-electron states of wave vector k e = k and k h = k − Q, respectively: where k Q is the exciton wave function in the (k, Q) representation. We can also define elementary excitonic states in real space: m,e a n,h |∅ n+R,e a n,h |∅ , with N the number of sites. In the expression above, m ∈ B are boron sites and n ∈ N are nitrogen sites. If the origin is located on an n site, the set of relative distance vectors R = m − n coincides therefore with the B = N + τ 3 lattice [40]. The state |Reh, Q is the Bloch state describing the motion of an electron-hole pair of size R. It can be checked that: For sake of simplicity, in the following we will frequently drop the eh label. The kinetic part of the excitonic Hamiltonian is obtained from the difference of the single particle Hamiltonians, H 0 eh = H 0 e − H 0 h . For the electron part the action of the Hamiltonian is given by: where h e (ρ e ) is the hopping integral associated with the vector ρ e connecting two conduction (boron) sites. In the same way the action of the hole Hamiltonian is given by: In the monolayer, the vector sets {ρ e } and {ρ h } are identical, and for the simple model used at Q = 0 [29], h h (ρ h ) = −h e (ρ e ) = t 2 ⊥ /2 = t ex /2, so that finally, and the effective hopping integral between 1nn at finite Q becomes: which reduces to t ex when Q = 0, as expected. The diagonal part of the Hamiltonian does not depend on Q and is therefore the same as before. Coming back to Eq. (B5) it is not difficult to include 2nn and even to introduce different hopping integrals between boron sites and between nitrogen sites. Finally, let RQ = R, Q| be the exciton wave function in the ( R, Q) representation, then the Bethe-Salpeter-Wannier (BSW) equation without the exchange term, becomes: where we have used the fact that R, Q| |R Q is diagonal in R [cfr. Eq. (3)]. Actually, the choice of the definition of k from the pair (k h , k e ) is not unique. Instead of the pair (k − Q, k), we could have chosen the pair (k − Q/2, k + Q/2). This should of course not modify the eigenvalues of the BSW equation, but modifies the phase of the "real space" state |Reh, Q defined in Eq. (B2). In particular, in the latter gauge, the effective hopping integral R, Q|H 0 eh |R , Q becomes real, which may be convenient in some cases. Finally, let us calculate the full wave function in real space: Remembering that we fixed the hole position on the valence site at the origin, r h = n = 0, so that: Actually r h should be taken slightly above n since the π wave functions φ(r ) vanish at the origin. Since the φ(r e − R) are localized on sites R (and also a little bit on their neighbors if genuine Wannier functions are used), we see that the full wave function can be (partially) represented as the superposition of localized contributions weighted by the amplitudes R Q . Notice also that the chosen gauge ensures that the wave function is invariant when Q is replaced by Q + G, where G is a reciprocal lattice vector. The previous formalism can be extended to bulk stackings by introducing interlayer hopping integrals and extending Eq. (B5) accordingly.

APPENDIX C: EXCHANGE CONTRIBUTIONS
Using the tight-binding Bloch states, the exchange kernel k e , k h |V |k e , k h involved in the Bethe-Salpeter equation can be written: Translation invariance implies Q = Q . Then the most important integrals are those where all involved lattice sites are as close to each other as possible. In usual treatments it is assumed that they are all identical, but here we know that the n and m sites belong to distinct triangular sublattices, so that the best we can do is to assume that n − m = τ and n − m = τ , where τ and τ are the shortest vectors joining the conduction sites to the fixed hole site. At least this is true to lowest order in t ⊥ /2 , in which case the Wannier functions are localized atomic orbitals centered on the lattice sites (see Ref. [29] for a discussion). Coming back to a real space discussion it is then found that V couple states |τ , Q and |τ , Q : with: (C3) Then we use the Fourier development of 1/r = q (4π/ q 2 )e iq·r , where is the volume, and then: having introduced the matrix elements M (τ , q) = d r ϕ e (r − τ )ϕ h (r )e iq·r and remembering that ϕ are real functions.
Let us now calculate the integrals M. To be consistent they should be neglected in the simplest tight-binding model, but actually the ϕ(r ) are here Wannier functions with components on the neighboring sites, where the ϕ at j (r ) are the genuine localized atomic orbitals (see Ref. [29] for a discussion). As a consequence in the macroscopic limit q → 0 we can write: which clearly corresponds to dipolar integrals, i.e., overlap integrals weighted by e ±iq·r . Finally, using (C6) and (C4) into (C2) one gets where at is the volume of the unit cell, and G are reciprocal lattice vectors of a system with periodic boundary conditions. This derivation is actually similar to the one by Qiu et al. [51] carried out within a k · p formalism and also to the result obtained by Cudazzo et al. [50] with a simplified model for the electronic structure not based on TB formalism. As derived, the sum over G is done on three dimensions (3D), but in the two-dimensional (2D) limit, the sum over the z components of G transforms into an integral which can be performed analytically. In the remaining sum over in-plane components G ⊥ , only the G ⊥ = 0 term is singular in the Q → 0 limit, and we finally obtain: where the area of the 2D unit cell is equal to 3 √ 3τ 2 , τ = |τ | = a/ √ 3 being the 1nn distance and a the lattice parameter of the triangular lattice.

APPENDIX D: PERTURBATIVE TREATMENT AT SMALL Q OF EXCHANGE TERMS IN THE MONOLAYER
Let us evaluate the effect of the exchange V on the first two (degenerate) excitons when one goes slightly away from Q = 0. At Q = 0 the two degenerate states are labeled + and − as in the main text. Let θ be the polar angle of Q, then the matrix of V in this ( + , − ) space reads where d = τ τ τ | is the dipolar matrix element of the | state. The eigenvalues of the matrix above are equal to 0 and 2. The lower branch (eigenvalue 0) is a "transverse" eigenstate for which Q · d = 0 at small Q, so it can be optically active. Moreover its dispersion is not affected by the exchange term. Instead the dispersion of the upper branch (eigenvalue 2) is affected by the singular component of the exchange, which is linear in Q [51,52], so the degeneracy is indeed lifted by the effect of the exchange interaction. Moreover it is a longitudinal eigenstate, hence dark at normal incidence but possibly active for oblique incidence. Notice that an opposite terminology is sometimes used, for example in electron energy loss spectroscopy the active modes are longitudinal so they are the bright ones.

APPENDIX E: A TOY MODEL FOR LOCALIZED EXCITONS: EXCITONS ON A KAGOME LATTICE AND FLAT BANDS
We have seen that at low Q the degeneracy of the ground state exciton of symmetry E is lifted into transversal and longitudinal states. The upper branch rises linearly with Q due to exchange effects. It is frequently argued that this is the main reason for the splitting of the states. Actually within a simple Wannier-Mott continuous model where the degenerate states are associated to separate valleys around point K or K the dispersion of the branches are expected to be similar since the effective masses are identical in both valleys. This is not necessarily the case when intervalley interactions are significant. We show below that the effect can be huge in the limit of strong localization of the exciton.
Let us assume therefore that the considered exciton wave function is confined to B-N 1nn pairs. Fixing the position of the hole, we have then three possible orientations of the excitonic pair. In the formalism used up to now, they are labeled by the corresponding vectors τ . In the exciton Hamiltonian the "kinetic energy" part related to the free motion of the hole and of the electron, the jump of the exciton from point τ to a neighboring site is precisely accounted for within the TB model by the hopping integral t ex . But we can also move the hole which will jump from one site to a neighboring one on its triangular sublattice. Usually it is very difficult to represent both motions of holes and electrons. Here this is possible because of the constraint that they should remain 1nn (on the honeycomb lattice). The method is to mark each pair by the position of its center. All these positions lie in fact on a so-called Kagome lattice where each site has four 1nn. It is easy to realize then that the motion of an excitonic pair on this lattice corresponds to first neighbor jumps on this Kagome lattice [ Fig. 12(a)]. The problem of describing the dispersion of the exciton states has been reduced to a single particle TB band problem on the Kagome lattice with 1nn interactions. The on-site matrix elements are all the same and equal to the Coulomb energy of the pair, taken here as the origin of energies.
The solution of this problem is known and shows very interesting features. The dispersion curves are shown in Fig. 12(b). The Brillouin zone is still the hexagonal one, and there are three branches since the Kagome lattice has three sites per unit cell. In particular there is a completely flat band. In general flat bands indicate the presence of localized states. This is of course what is obtained if the interatomic jumps are forbidden. What is surprising here is that jumps are allowed, but a basis of localized states should exist. They do exist and are actually localized on the hexagons of the underlying honeycomb lattice. This has been discussed in many places. Flat bands may produce surprising effects as easy self-localization of extended states in the presence of small perturbations. For a review see for instance Ref. [53]; see also Ref. [60]. It is tempting to apply this model to the behavior of our ground state excitons by looking at the two lowest states of the Kagome lattice. The third and highest level is unphysical in our context since anyway there are other, more extended, excitons at high energy and the continuum of single particle excitations. By this model we demonstrate that if excitons are very localized a flat band may appear, which undoubtedly will induce peculiar effects. In the case of hBN, although the ground state exciton is fairly localized, we are clearly not in this extreme limit and the difference in dispersion of the two branches at low Q is principally due to exchange contributions. This does not mean that the Kagome model is useless, since for example it provides explicit solutions for the energy and the wave functions.

APPENDIX F: VARYING THE INTERPLANE DISTANCE
We have decided not to perform structural optimization in the bulk structures. For the comparison to be meaningful, we have used the same cell parameters in the three structures. In particular, the interlayer distance c = 3.25 Å has been used in the three cases. However, we have explored in the case of the AB stacking the effect on the band structure and the exciton dispersion of a variation of the interlayer spacing. In Fig. 13, we show that a variation from c = 3.25 Å to c = 3.35 Å, which corresponds to an increase of 3%, does not change the conclusions of the main text.
It can be seen that the changes in the exciton dispersion induced by the variation of c are negligible. In particular, the prediction that the nature of the gap changes from indirect to direct when the electron-hole interaction is taken into account is not compromised by this change. In fact, the value c = 3.35 Å is probably closer to the actual interplane distance in this material, so we expect this effect to be larger than what was predicted in the main text on the basis of c = 3.25 Å. We note, by the way, that the difference between IP-transition data of the two sets (distance between dashed curves in the figure) is larger than the corresponding difference between excitonic data (solid curves). This can be seen for instance in the or in the K points. It indicates the expected trend of an increase of both the IP-transition energies (i.e., an increase of the gap) and of the exciton binding energy at higher c.

APPENDIX G: ROBUSTNESS OF THE PREDICTION
Here we compare two calculations of the exciton dispersion in the AA stacking done with different approximations for the quasiparticle correction and the DFT exchange-correlation potentials. The aim is clearly not to analyze the differences of the two approaches neither to make a thorough comparison of the two results. The objective is to assess the robustness of the results discussed in the main text.
In Fig. 14 we report four dispersion curves of the first exciton in the AA stacking for Q K. All excitons have been aligned at to better visualize the variations on the energy dispersion. Two calculations are based on LDA Kohn-Sham structure where the quasiparticle corrections have been approximated with a perturbative GW correction or a scissor operator. In the other two, the quasiparticle corrections have been computed within the same approximations but from PBE Kohn-Sham results. The dispersions computed within the scissor operator are quite different, but still below 0.2 eV. However the two dispersions basically coincide when quasiparticle corrections are modelled within the GW approximation.
Furthermore it is remarkable that the two sets of calculations have been done with different codes. The LDA set of simulations have been obtained using the plane-wave codes EXC [39] and ABINIT [38], as described in II A 2. The PBE set of simulations have been obtained using the GPAW [36] package. The very good agreement of the two GW results indicate that all parameters have been carefully converged in all calculations.
This comparison demonstrates the reliability and the robustness of our results, in particular regarding the claim of direct exciton formation in the AB bulk phase.