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

Lorenzo Sponza, Hakim Amara, Claudio Attaccalite, Sylvain Latil, Thomas Galvani, Fulvio Paleari, Ludger Wirtz, François Ducastelle

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Lorenzo Sponza, Hakim Amara, Claudio Attaccalite, Sylvain Latil, Thomas Galvani, et al.. Direct and indirect excitons in boron nitride polymorphs: a story of atomic configuration and electronic correlation. Physical Review B: Condensed Matter and Materials Physics (1998-2015), 2018, 98 (12), pp.125206. 10.1103/PhysRevB.98.125206. hal-01820825

HAL Id: hal-01820825

https://amu.hal.science/hal-01820825

Submitted on 31 Jan 2019

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Direct and indirect excitons in boron nitride polymorphs: A story of atomic configuration and electronic correlation

Lorenzo Sponza,¹ Hakim Amara,¹ Claudio Attaccalite,²,³ Sylvain Latil,⁴ Thomas Galvani,⁵

Fulvio Paleari,⁵ Ludger Wirtz,⁵ and François Ducastelle¹

¹Laboratoire d'Etude des Microstructures, ONERA-CNRS, UMR104, Université Paris-Saclay, BP 72, 92322 Châtillon Cedex, France

²Aix Marseille University and CNRS, CINAM, UMR7325, Campus de Luminy, case 913, 13288 Marseille, France

³University of Rome Tor Vergata, Rome, Italy

⁴CEA, IRAMIS, SPEC, GMT, 91191 Gif-sur-Yvette, France

⁵Physics and Materials Science Research Unit, University of Luxembourg, 162a avenue de la Faïencerie, L-1511 Luxembourg, Luxembourg

<figure><img src="image_3.png" /></figure>

(Received 5 July 2018; published 28 September 2018)

We present a detailed discussion of 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 at the level of *ab initio* many-body perturbation theory (GW and Bethe Salpeter equation) and of a purposely developed 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. As a result the electron-hole dispersion is flatter than the one obtained from the band structure. 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.

DOI: 10.1103/PhysRevB.98.125206

## I. INTRODUCTION

Hexagonal boron nitride (hBN) is a wide band gap semiconductor ($E_g > 6$ eV [1]) whose exceptional optical properties are attracting a growing interest [2–5]. Many recent experimental [6–20] and theoretical studies [11,21–26] focused on the understanding of the high intensity of the luminescence signal in the bulk phase. This property is indeed in apparent contradiction with the indirect band gap and the observation of the so-called Stokes shift between absorption and luminescence [27].

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 $\Gamma M$ and the $\Gamma K$ lines [32].

In this paper we push the analysis further, discussing in detail the dispersion of the excitons especially in the $\Gamma 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 in-plane 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
<figure><img src="image_4.png" /><figcaption>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₃N₃ hexagons highlights the difference between the three bulk structures.</figcaption></figure>

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.

### A. Ab initio simulations

#### 1. 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].

#### 2. 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_0W_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** with **k** = (1/3, 1/3, 1/3) expressed in the reciprocal coordinates of the rhombohedral cell.

The macroscopic dielectric function $\epsilon_M(q, \omega)$ 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.
<figure><img src="image_5.png" /><figcaption>FIG. 2. (a) The two triangular sublattices $\Lambda_B$ (green) and $\Lambda_N$ (gray) forming the hBN honeycomb lattice. Red arrows mark the 1nn vectors $\tau_\alpha$ and the positions $\mathbf{R}, \mathbf{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_\perp, t_{2\perp}, t_\parallel$, and $t_{2\parallel}$ in the AA' stacking.</figcaption></figure>

## B. Tight-binding model

### 1. 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 $\Lambda_B$ and $\Lambda_N$, either connecting all sites of B or N atoms, respectively. Vectors $\tau_{\{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 $\pm\Delta$ 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_\perp$ and $t_{2\perp}$, 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 $\mathbf{k}_\perp$ along the $KM$ direction as:

$$E_{\mathbf{k}_{\perp} e(h)} \approx \pm \underbrace{\left( \Delta + \frac{|t_{\perp} \gamma(\mathbf{k}_{\perp})|^2}{2\Delta} \right)}_{\text{first neighbors}} + \underbrace{t_{2\perp} (|\gamma(\mathbf{k}_{\perp})|^2 - 3)}_{\text{second neighbors}}, \quad (1)$$

where $\gamma(\mathbf{k}_\perp) = \sum_{\alpha=1,2,3} e^{i\mathbf{k}_\perp \cdot \tau_\alpha}$ and $\mathbf{k}_\perp$ is strictly in-plane. In the above expression the (+) sign is for conduction $e$ states localized on $\Lambda_B$ sites, and (-) is for valence $h$ states localized on $\Lambda_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 one-electron $\pi$ bands with atomiclike Wannier functions allows us to work directly with excitonic Wannier equations in real space [35,41–44]. Assuming the hole to be fixed on a N site, and using relative coordinates $\mathbf{R}$ for the electron-hole pair, we have reduced the problem to the one of the electron hopping on sites of $\Lambda_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_{eh}^0 + \Xi + V$ contains a kinetic part $H_{eh}^0$, a screened Coulomb term $\Xi$ that can be taken as a fitting parameter and an exchange term $V \propto Q$ that actually vanishes in the $Q=0$ limit.

We extend now the model to $Q \neq 0$ hence allowing the exciton to move across the layer. Through the definition of appropriate Bloch states $|\mathbf{R}, \mathbf{Q}\rangle$ defined in Appendix B, the momentum $\mathbf{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 $\mathbf{Q}$ and the relative electron-hole distance in real space $\mathbf{R}$. For the monolayer, it can be demonstrated that the simple model above can be extended to a Hamiltonian $H_{eh} = H_{eh}^0 + \Xi + V$ where the kinetic part reads:

$$\begin{align}
\langle \mathbf{R}', \mathbf{Q} | H_{eh}^0 | \mathbf{R}, \mathbf{Q} \rangle &= \begin{cases} 2\Delta + 3t_{\perp}^2/\Delta & \text{if } \mathbf{R} = \mathbf{R}' \\ \frac{t_{\perp}^2}{2\Delta}(1 + e^{i\mathbf{Q} \cdot (\mathbf{R} - \mathbf{R}')}) & \text{if } |\mathbf{R} - \mathbf{R}'| = \tau, \\ 0 & \text{otherwise} \end{cases} \tag{2} \\
\tau &\text{ being } |\tau_\alpha| \text{ for any } \alpha = 1, 2, 3. \text{ The Coulomb part reads}
\end{align}$$

$$\langle \mathbf{R}', \mathbf{Q} | \Xi | \mathbf{R}, \mathbf{Q} \rangle = \begin{cases} \Xi_R & \text{for } \mathbf{R} = \mathbf{R}' \\ 0 & \text{otherwise,} \end{cases} \quad (3)$$

with $\Xi_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 \to 0$ [43]. In our TB model it produces effective interactions $\langle \mathbf{R}', \mathbf{Q} | V | \mathbf{R}, \mathbf{Q} \rangle$ 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 $\mathbf{R}$, but have short- and long-range components with respect to the propagation of the center of mass $\mathbf{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 dipole-dipole 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:

$$\begin{align}
\langle \mathbf{R}', \mathbf{Q} | V | \mathbf{R}, \mathbf{Q} \rangle & \approx \begin{cases} 0 & \text{for } \mathbf{R} \text{ or } \mathbf{R}' \notin \{\tau_1, \tau_2, \tau_3\} \\ \frac{J}{4Q\tau} (e^{-i\mathbf{Q}\cdot\mathbf{R}} - 1)(e^{i\mathbf{Q}\cdot\mathbf{R}'} - 1) & \text{otherwise,} \end{cases} \tag{4}
\end{align}$$

where $J$ will be considered here as a parameter to be fitted to ab initio calculations, and where $Q = |\mathbf{Q}|$ and $\tau = |\tau_\alpha|$ for any $\alpha$.

### 2. Bulk structures

When passing to bulk systems, some aspects complicate the model. The first is that the corresponding Bloch states are characterized by $\mathbf{k} = \mathbf{k}_\perp + \hat{z}k_\parallel$. Moreover in the AA' and AB stacking, the basis is formed of four Bloch states
<figure><img src="image_6.png" /><figcaption>FIG. 3. The four configurations $R_{α,β}$ of electron-hole pairs in AA' stacking. The origin is marked with a blue spot on the $Λ_N$ sublattice of layer 1. Electrons and holes localize on different $Λ + u$ sublattices with $u$ appropriate translation vector ($Λ$ is defined as $Λ_N$ on layer 1, and $Λ_B = Λ + τ_3$ on layer 1) [40]. The relative distance $R_{α,β}$ belongs also to a shifted lattice $Λ + u$.</figcaption></figure>

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(b).

Once properly generalized and retaining only $1nn$ intra- and interplane couplings, the energy of the $π^*$ conduction states reads:

$$E_{ke}^{AA'} \approx \Delta + \frac{(t_{\perp}|\gamma(\mathbf{k}_{\perp})| \pm 2t_{\parallel}\cos(k_{\parallel}c))^2}{2\Delta}, \quad (5)$$

$$E_{ke}^{ABC} \approx \Delta + \frac{|t_{\perp}\gamma(\mathbf{k}_{\perp}) + t_{\parallel}e^{ik_{\parallel}c}|^2}{2\Delta}, \quad (6)$$

$$E_{ke}^{AB} \approx \Delta + \frac{(t_{\parallel}\cos(k_{\parallel}c) \pm \sqrt{(t_{\parallel}\cos(k_{\parallel}c))^2 + t_{\perp}^2|\gamma(\mathbf{k}_{\perp})|^2})^2}{2\Delta}, \quad (7)$$

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 $\hat{z}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\rangle$.

(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).

### III. THE FREE-STANDING SINGLE LAYER

#### 1. 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
<figure><img src="image_7.png" /><figcaption>FIG. 4. The single-particle band structure of free-standing h-BN single layer (a) *ab initio* DFT (dotted red), quasiparticle (solid black) calculation, and optimal TB fitting (dashed blue). (b) Dependence of the TB band structure on the parameter $t_{2\perp}$.</figcaption></figure>

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–30,47–50].

The TB model introduced previously gives the best fit to the KM region of the *ab initio* band structure with $\Delta = 3.625$ eV, $t_{\perp} = -2.33$ eV, and $t_{2\perp} = -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\perp}$ 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\perp}$ have been used. Note that the (direct) gap at K is not modified because the symmetry of the crystal leads to $\gamma(\mathbf{k}_{\perp} = \mathbf{K}) = 0$. To conclude, this simple TB model is able to reproduce the $\pi$ states in the regions where excitons are relevant for the optical properties.

## 2. Exciton dispersion

The dispersion curves $E_{exc}(Q)$ obtained from *ab initio* calculations are shown in Fig. 5(a). When $Q = \Gamma$ 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 $\Gamma M$ and $\Gamma 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 \to 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 $\Gamma K$ and in particular the fact that the energy at $Q = \Gamma$ 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 = \Gamma$ and $Q = K$ come from vertical $K \to K$ and slant $K \to K'$ transitions, respectively. Since the single-particle states $K$ and $K'$ have the same energy, the resulting dispersion of the exciton along $\Gamma K$ is expected to attain the same value at the extrema of the path. Our results along $\Gamma 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 $|\Phi^+\rangle$ and $|\Phi^-\rangle$ can be taken as two “circular” states whose components on the three $\tau_{1,2,3}$ sites $\langle\tau|\Phi^+\rangle$ and $\langle\tau|\Phi^-\rangle$ are proportional to the cubic roots of unity $(1, \omega, \omega^2)$ and $(1, \omega^2, \omega)$, 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 Fig. 5(b) shows that the degeneracy between the two excitons
<figure><img src="image_8.png" /><figcaption>FIG. 5. The exciton dispersion for $q \parallel \Gamma M$ and $q \parallel \Gamma K$ computed with (a) ab initio and (b) tight-binding methods. In the (b) panel, results without 2nn hopping and including the exchange term are reported in black, calculations without the exchange are reported in dashed-red ($t_{2\perp} = 0$ eV) and solid blue ($t_{2\perp} = -0.4$ eV). Bottom panels: Real-space representation of the tight-binding electronic density $|\Phi_{RQ}|^2$ of the first two excitons with the hole located at the center of the axis (blue dot). At $\Gamma$, we report in dashed red the direction parallel to the Q vector of the series and a red hexagon marks the real-space hexagonal cell.</figcaption></figure>

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 $\pi$ 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 $\Gamma K$, so we can ascribe these bands to $\pi \to \pi^*$ excitations. Instead, along $\Gamma M$ these bands are predicted only in the *ab initio* solution, so we
can advance the hypothesis that they are of the $\sigma \to \pi^*$ or $\pi \to \sigma^*$ type as the $\sigma$ 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 $\tau_{2\perp} = -0.4$ eV, the main effect is to decrease the dispersion of the lowest exciton along $\Gamma 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:

$$|\Phi\rangle = \sum_R \Phi_{RQ} |R, Q\rangle. \qquad (8)$$

As discussed in Appendix B, the basis $|R, Q\rangle$ and therefore the coefficients $\Phi_{RQ}$ are not uniquely defined. With our standard definition ($k_e, k_h$) = ($k, k - Q$), $\Phi_{RQ}$ is in general a complex quantity. In the bottom part of Fig. 5 we show a map of the intensities $|\Phi_{RQ}|^2$ (electronic densities) which are gauge invariant. Here we consider more particularly the first two excitons (degenerate at $Q = 0$) along $\Gamma M$ and $\Gamma K$.

Consider first the $\Gamma 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_\Phi$ 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 $\Gamma M$. In the $\Gamma K$ case (right panels), we no longer expect definite symmetries except at $\Gamma$ 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 $\Gamma 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 $\Gamma 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 $\Gamma 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.
<figure><img src="image_9.png" /><figcaption>FIG. 6. The GW band structure of the three bulk phases drawn in the hexagonal Brillouin zone.</figcaption></figure>

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

<table><caption>TABLE I. Dispersion along *KM* and splitting at *M* of the GW bottom conduction (see bottom panels of Fig. 6). All energies are in eV.</caption><thead><tr><th>Stacking</th><th>Dispersion</th><th>Splitting</th></tr></thead><tbody><tr><td>AA'</td><td>0.75</td><td>1.68</td></tr><tr><td>AB</td><td>0.10</td><td>1.34</td></tr><tr><td>ABC</td><td>0.51</td><td>1.11</td></tr></tbody></table>

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 $\Delta$ and $t_{\perp}$ as in the monolayer [57], and we add only the first-neighbor interlayer hopping $t_{\parallel}$. Then the two conduction $\pi^*$ eigenvalues can be approximated by Eq. (5) (and with opposite sign the two valence $\pi$ states). From the equation we get that the splitting between the $\pi$ states vanish either when $\gamma(k_{\perp}) = 0$ (HK line, not shown here) or when $k_{\parallel}c = \pm\pi/2$, i.e., on the upper and lower faces of the Brillouin zone (AHL line). In fact since the periodicity along $z$ is 2c, then $k_{\parallel}c \in [-\pi/2, \pi/2]$. Actually the doubly degenerate state dispersing along *HL* has exactly the average energy of the two split branches on $\Gamma KM$. The model reproduces the bottom conduction at *M*, where the splitting between the two conduction bands is $4t_{\perp}t_{\parallel}/\Delta$ ($|\gamma(M)| = 1$). The fitting procedure to the *ab initio* band structure yields $t_{\parallel} \simeq 0.5$ eV, i.e., $t_{\parallel}/t_{\perp} \simeq 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_{\parallel}\sqrt{t_{\parallel}^2 + t_{\perp}^2}/\Delta$ 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 $3t_{\perp}t_{\parallel}/\Delta$ 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\perp}$) and between the planes ($t_{2\parallel}$), as discussed in Sec. III 1. The latter term accounts for the difference of splitting, larger in the conduction band than

<table><caption>TABLE II. The smallest direct and indirect gaps (GW and LDA) extracted from the bandplot of the three structures. Energies are in eV.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Smallest direct</th><th colspan="3">Indirect</th></tr><tr><th>GW</th><th>LDA</th><th>point</th><th>GW</th><th>LDA</th><th>points</th></tr></thead><tbody><tr><td>AA'</td><td>6.28</td><td>4.46</td><td><i>M</i></td><td>5.80</td><td>4.02</td><td><i>T</i><i>M</i></td></tr><tr><td>AB</td><td>6.13</td><td>4.31</td><td><i>K</i></td><td>6.01</td><td>4.20</td><td><i>T</i><i>M</i></td></tr><tr><td>ABC</td><td>5.82</td><td>4.51</td><td><i>K</i></td><td>5.27</td><td>3.95</td><td><i>K</i><i>M</i></td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>FIG. 7. Electronic band structure of the AA' phase. Influence of the TB 2nn interactions: (a) in-plane $t_{2\perp}$ and (b) between planes $t_{2\|}$.</figcaption></figure>

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\perp} \simeq -0.4$ eV and $t_{2\|} \simeq -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.

### 1. *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 $\Gamma 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 = \Gamma$, corresponding to the smallest direct gap at $M$, to a minimum of 5.80 eV at $1/\AA^{-1}$, corresponding to the indirect gap $TM$, 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

<figure><img src="image_11.png" /><figcaption>FIG. 8. From top to bottom: Dispersion of the first excitons in the AA', ABC, and AB stacking for $Q \parallel \Gamma K$ (black bullets: GW-BSE). The size of the spots is proportional to the logarithm of the intensity in arbitrary units. The dispersion of the lowest independent-particle transition (free-carrier dispersion) is reported with red empty circles (GW-IP).</figcaption></figure>

0.11 eV. In fact the exciton binding energy $E_b(Q)$ exhibits a strong dependence on the exchanged momentum, varying by almost 0.4 eV throughout the $\Gamma K$ path. The principal consequence is that the difference between direct and indirect
<table><caption>TABLE III. Dispersion of the first free-carrier transition (GW-IP) and first exciton (GW-BSE). Minimum and maximum of the binding energy $E_b(Q)$. All energies are in eV.</caption><thead><tr><th></th><th>AA'</th><th>ABC</th><th>AB</th></tr></thead><tbody><tr><td>GW-IP</td><td>0.48</td><td>0.38</td><td>0.12</td></tr><tr><td>GW-BSE</td><td>0.11</td><td>0.14</td><td>-0.02</td></tr><tr><td>Min. of E<sub>b</sub></td><td>0.30</td><td>0.35</td><td>0.35</td></tr><tr><td>Max. of E<sub>b</sub></td><td>0.67</td><td>0.63</td><td>0.50</td></tr></tbody></table>

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 $\Gamma 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 $|\Phi_Q(\mathbf{r}_e, \mathbf{r}_h)|^2$ as a function of $\mathbf{r}_e$ having fixed $\mathbf{r}_h$ on an arbitrary N atom. The electronic densities at $Q = \Gamma$, $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 $\Gamma$. 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.

## 2. *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/(\epsilon_j R)$ with effective dielectric constants $\epsilon_\perp = 6$ within the planes and $\epsilon_\parallel = 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_\perp$ term, equal to -2.33 eV. The corresponding exciton dispersion is that of the monolayer, but each curve is doubly degenerate because there

<figure><img src="image_12.png" /><figcaption>FIG. 9. *Ab initio* electronic density ($|\Phi_Q(\mathbf{r}_e, \mathbf{r}_h)|^2$) with fixed hole on a N atom at different exchanged momenta in the AA' bulk phase. At finite $Q$, a red arrow indicates a direction parallel to $Q$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 10. Tight-binding dispersion along $\Gamma K$ in the AA' stacking. The 1$n$ and 2$n$ hopping terms $t_{\perp}$, $t_{\parallel}$, $t_{2\perp}$, and $t_{2\parallel}$ are switched on sequentially from panel (a) to panel (d).</figcaption></figure>

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_{\parallel} = 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 $\Gamma 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\perp} = -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 2$n$ integrals $t_{2\parallel} = -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 electron-hole asymmetry play crucial roles to account for the direct-indirect 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 single-particle 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 $\Gamma M$ are in agreement with previously published data [29,46,50]. We also provide the exciton dispersion along $\Gamma 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.

<figure><img src="image_14.png" /><figcaption>FIG. 11. Excitonic dispersion along $\Gamma K$ for the three different stackings, AA', ABC, and AB, from left to right.</figcaption></figure>
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 second-nearest-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 = \Gamma$ 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.

## ACKNOWLEDGMENTS

The authors want to acknowledge Dr. J. Barjon of UVSQ-CNRS, Université Paris-Saclay. The French National Agency for Research (ANR) is acknowledged for funding this work under the project GoBN (Graphene on Boron Nitride Technology), Grant No. ANR-14-CE08-0018. F.P., T.G. and L.W. acknowledge support from the National Research Fund, Luxembourg (Projects No. EXCPHON/11280304 and No. INTER/RCN/16/11402387/2Defect). The research leading to these results has received funding from the European H2020 Framework Program under Grant Agreements No. 696656 Graphene Core1 and No. 785219 Graphene Core2.

## 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

$$E_{\mathbf{k}_{\perp} e(h)} = \pm \sqrt{\Delta^2 + t_{\perp}^2 |\gamma(\mathbf{k}_{\perp})|^2 + t_{2\perp} (|\gamma(\mathbf{k}_{\perp})|^2 - 3)}, \quad (A1)$$

where the signs (+) and (−) are for conduction and valence states, respectively. The approximate expression (1) valid for $\Delta \gg t_{\perp}^2 |\gamma(\mathbf{k}_{\perp})|^2$ holds very well in the KM line. The

<table><caption>TABLE IV. Optimal TB parameters. Energies in eV</caption><thead><tr><th>Δ</th><th>t<sub>⊥</sub></th><th>t<sub>2⊥</sub></th><th>t<sub>‖</sub></th><th>t<sub>2‖</sub></th></tr></thead><tbody><tr><td>3.625</td><td>-2.33</td><td>-0.4</td><td>0.5</td><td>-0.1</td></tr></tbody></table>

optimization of the parameters $\Delta = 3.625$ eV, $t_{\perp} = -2.33$ eV, and $t_{2\perp} = -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:

$$E_{\mathbf{k}e}^{AA'} = \sqrt{\Delta^2 + [t_{\perp}|\gamma(\mathbf{k}_{\perp})| \pm 2t_{\parallel}\cos(k_{\parallel}c)]^2 + t_{2\perp}(|\gamma(\mathbf{k}_{\perp})|^2 - 3) \pm 2t_{2\parallel}|\gamma(\mathbf{k}_{\perp})|\cos(k_{\parallel}c)}, \quad (A2)$$

and hole states:

$$E_{\mathbf{k}h}^{AA'} = -\sqrt{\Delta^2 + [t_{\perp}|\gamma(\mathbf{k}_{\perp})| \pm 2t_{\parallel}\cos(k_{\parallel}c)]^2 + t_{2\perp}(|\gamma(\mathbf{k}_{\perp})|^2 - 3) \pm 2t_{2\parallel}|\gamma(\mathbf{k}_{\perp})|\cos(k_{\parallel}c)}. \quad (A3)$$

The approximate equation (5) reported in the main text includes only 1nn terms. While $\Delta$, $t_{\perp}$ and $t_{2\perp}$ have been kept as in the monolayer, the values $t_{\parallel} = 0.5$ eV and $t_{2\parallel} = -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 \neq 0$ (WITH NO EXCHANGE CONTRIBUTIONS)

A general exciton state $|\Phi_Q\rangle$ is obtained by combining electron and hole one-electron states of wave vector $\mathbf{k}_e = \mathbf{k}$ and $\mathbf{k}_h = \mathbf{k} - Q$, respectively:

$$|\Phi_Q\rangle = \sum_k \Phi_{kQ} |keh, Q\rangle \\ |\mathbf{k}_e, \mathbf{k}_h\rangle = |keh, Q\rangle = a_{\mathbf{k}_e}^{\dagger} a_{\mathbf{k}_e - Q} |\emptyset\rangle, \quad (B1)$$

where $\Phi_{kQ}$ is the exciton wave function in the $(\mathbf{k}, Q)$ representation. We can also define elementary excitonic states in real space:

$$|\mathbf{Reh}, \mathbf{Q}\rangle = \frac{1}{\sqrt{N}} \sum_{\substack{n,m \\ m-n=R}} e^{-i\mathbf{Q} \cdot \mathbf{n}} a_{\mathbf{m},e}^{\dagger} a_{\mathbf{n},h} |\emptyset\rangle \\ = \frac{1}{\sqrt{N}} \sum_{n} e^{-i\mathbf{Q} \cdot \mathbf{n}} a_{\mathbf{n}+\mathbf{R},e}^{\dagger} a_{\mathbf{n},h} |\emptyset\rangle,$$
with $N$ the number of sites. In the expression above, $m \in \Lambda_B$ are boron sites and $n \in \Lambda_N$ are nitrogen sites. If the origin is located on an $n$ site, the set of relative distance vectors $\mathbf{R} = \mathbf{m} - \mathbf{n}$ coincides therefore with the $\Lambda_B = \Lambda_N + \tau_3$ lattice [40]. The state $|\mathbf{Reh}, \mathbf{Q}\rangle$ is the Bloch state describing the motion of an electron-hole pair of size $\mathbf{R}$. It can be checked that:

$$|\mathbf{keh}, \mathbf{Q}\rangle = \frac{1}{\sqrt{N}} \sum_{\mathbf{R}} e^{i\mathbf{k}\cdot\mathbf{R}} |\mathbf{Reh}, \mathbf{Q}\rangle . \quad (B2)$$

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_{eh}^0 = H_e^0 - H_h^0$. For the electron part the action of the Hamiltonian is given by:

$$H_e^0 |\mathbf{R}, \mathbf{Q}\rangle = \sum_{\rho_e} h_e(\rho_e) |\mathbf{R} + \rho_e, \mathbf{Q}\rangle, \quad (B3)$$

where $h_e(\rho_e)$ is the hopping integral associated with the vector $\rho_e$ connecting two conduction (boron) sites. In the same way the action of the hole Hamiltonian is given by:

$$H_h^0 |\mathbf{R}, \mathbf{Q}\rangle = \sum_{\rho_h} h_h(\rho_h) e^{-i\mathbf{Q}\cdot\rho_h} |\mathbf{R} - \rho_h, \mathbf{Q}\rangle. \quad (B4)$$

In the monolayer, the vector sets $\{\rho_e\}$ and $\{\rho_h\}$ are identical, and for the simple model used at $\mathbf{Q}=0$ [29], $h_h(\rho_h) = -h_e(\rho_e) = t_{\perp}/2\Delta = t_{ex}/2$, so that finally,

$$H_{eh}^0 |\mathbf{R}, \mathbf{Q}\rangle = \sum_{\rho_e} (h_e(\mathbf{\rho}) - h_h(\mathbf{\rho})e^{i\mathbf{Q}\cdot\mathbf{\rho}}) |\mathbf{R} + \mathbf{\rho}, \mathbf{Q}\rangle, \quad (B5)$$

and the effective hopping integral between $1nn$ at finite $\mathbf{Q}$ becomes:

$$\langle \mathbf{R}', \mathbf{Q} | H_{eh}^0 | \mathbf{R}', \mathbf{Q} \rangle = (h_e(\mathbf{R}' - \mathbf{R}) - h_h(\mathbf{R}' - \mathbf{R})e^{i\mathbf{Q} \cdot (\mathbf{R}' - \mathbf{R})}) \\ = \frac{t_{ex}}{2} (1 + e^{i\mathbf{Q} \cdot (\mathbf{R}' - \mathbf{R})}),$$

which reduces to $t_{ex}$ when $\mathbf{Q}=0$, as expected. The diagonal part of the Hamiltonian does not depend on $\mathbf{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 $\Phi_{RQ} = \langle \mathbf{R}, \mathbf{Q} | \Phi \rangle$ be the exciton wave function in the $(\mathbf{R}, \mathbf{Q})$ representation, then the Bethe-Salpeter-Wannier (BSW) equation without the exchange term, becomes:

$$E \Phi_{RQ} = \sum_{\mathbf{R}'} \langle \mathbf{R}, \mathbf{Q} | H_{eh}^0 | \mathbf{R}' , \mathbf{Q} \rangle \Phi_{R'Q} + \Xi_R \Phi_{RQ}, \quad (B6)$$

where we have used the fact that $\langle \mathbf{R}, \mathbf{Q} | \Xi | \mathbf{R}' \mathbf{Q} \rangle$ is diagonal in $\mathbf{R}$ [cfr. Eq. (3)]. Actually, the choice of the definition of $\mathbf{k}$ from the pair $(\mathbf{k}_h, \mathbf{k}_e)$ is not unique. Instead of the pair $(\mathbf{k}-\mathbf{Q}, \mathbf{k})$, we could have chosen the pair $(\mathbf{k}-\mathbf{Q}/2, \mathbf{k}+\mathbf{Q}/2)$. This should of course not modify the eigenvalues of the BSW equation, but modifies the phase of the “real space” state $|\mathbf{Reh}, \mathbf{Q}\rangle$ defined in Eq. (B2). In particular, in the latter gauge, the effective hopping integral $\langle \mathbf{R}, \mathbf{Q} | H_{eh}^0 | \mathbf{R}', \mathbf{Q} \rangle$ becomes real,

which may be convenient in some cases. Finally, let us calculate the full wave function in real space:

$$\begin{aligned} \Phi_Q(\mathbf{r}_h, \mathbf{r}_e) &= \sum_{\mathbf{R}} \Phi_{RQ} \langle \mathbf{r}_h, \mathbf{r}_e | \mathbf{R}, \mathbf{Q} \rangle \\ &= \frac{1}{\sqrt{N}} \sum_{\mathbf{R}} \Phi_{RQ} \sum_{\mathbf{n}} e^{-i\mathbf{Q}\cdot\mathbf{n}} \phi_e(\mathbf{r}_e - \mathbf{n} - \mathbf{R}) \\ & \quad \times \phi_h(\mathbf{r}_h - \mathbf{n}). \end{aligned} \quad (B7)$$

Remembering that we fixed the hole position on the valence site at the origin, $\mathbf{r}_h = \mathbf{n} = 0$, so that:

$$\Phi_Q(\mathbf{r}_h = 0, \mathbf{r}_e) \propto \sum_{\mathbf{R}} \Phi_{RQ} \phi_e(\mathbf{r}_e - \mathbf{R}). \quad (B8)$$

Actually $\mathbf{r}_h$ should be taken slightly above $\mathbf{n}$ since the $\pi$ wave functions $\phi(\mathbf{r})$ vanish at the origin. Since the $\phi(\mathbf{r}_e - \mathbf{R})$ are localized on sites $\mathbf{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 $\Phi_{RQ}$. Notice also that the chosen gauge ensures that the wave function is invariant when $\mathbf{Q}$ is replaced by $\mathbf{Q} + \mathbf{G}$, where $\mathbf{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 $\langle \mathbf{k}'_e, \mathbf{k}'_h | \mathbf{V} | \mathbf{k}_e, \mathbf{k}_h \rangle$ involved in the Bethe-Salpeter equation can be written:

$$\begin{aligned} \langle \mathbf{k}'_e, \mathbf{k}'_h | \mathbf{V} | \mathbf{k}_e, \mathbf{k}_h \rangle &= \int d\mathbf{r} d\mathbf{r}' \phi_{k'_e e}^*(\mathbf{r}') \phi_{k'_h e}(\mathbf{r}') \frac{2e^2}{|\mathbf{r} - \mathbf{r}'|} \phi_{ke}(\mathbf{r}) \phi_{kh}^*(\mathbf{r}) \\ &= \frac{1}{N^2} \sum_{\mathbf{n},m,n',m'} \exp\{i[\mathbf{k} \cdot \mathbf{n} - (\mathbf{k} - \mathbf{Q}) \cdot \mathbf{m} - \mathbf{k}'n'] \\ & \quad + (\mathbf{k} - \mathbf{Q}') \cdot \mathbf{m}']\} \int d\mathbf{r} d\mathbf{r}' \phi_e(\mathbf{r}' - \mathbf{n}') \phi_h^*(\mathbf{r}' - \mathbf{m}') \\ & \quad \times \frac{2e^2}{|\mathbf{r} - \mathbf{r}'|} \phi_e(\mathbf{r} - \mathbf{n}) \phi_h^*(\mathbf{r} - \mathbf{m}). \end{aligned} \quad (C1)$$

Translation invariance implies $\mathbf{Q} = \mathbf{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 $\mathbf{n}$ and $\mathbf{m}$ sites belong to distinct triangular sublattices, so that the best we can do is to assume that $\mathbf{n} - \mathbf{m} = \mathbf{\tau}$ and $\mathbf{n}' - \mathbf{m}' = \mathbf{\tau}'$, where $\mathbf{\tau}$ and $\mathbf{\tau}'$ are the shortest vectors joining the conduction sites to the fixed hole site. At least this is true to lowest order in $t_{\perp}/2\Delta$, 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 $\mathbf{V}$ couple states $|\mathbf{\tau}, \mathbf{Q}\rangle$ and $|\mathbf{\tau}', \mathbf{Q}\rangle$:

$$\langle \mathbf{\tau}', \mathbf{Q} | \mathbf{V} | \mathbf{\tau}, \mathbf{Q} \rangle = \sum_{\rho} J_{\mathbf{\tau}'\mathbf{\tau}}(\rho) e^{-i\mathbf{Q}\cdot\mathbf{\rho}}, \quad (C2)$$
with:

$$J_{\tau'\tau}(\rho) = \int dr dr' \frac{\varphi_e(\mathbf{r}-\tau)\varphi_h(\mathbf{r})2e^2\varphi_e(\mathbf{r}'-\tau')\varphi_h(\mathbf{r}')}{|\mathbf{r}'-\mathbf{r}+\rho|} \quad (C3)$$

Then we use the Fourier development of $1/r = \sum_q (4\pi/\Omega q^2)e^{iq \cdot r}$, where $\Omega$ is the volume, and then:

$$J_{\tau'\tau}(\rho) = \sum_{\mathbf{q}} \frac{8\pi e^2}{\Omega q^2} e^{i\mathbf{q} \cdot \mathbf{r}} M(\tau, \mathbf{q})^* M(\tau', \mathbf{q}), \quad (C4)$$

having introduced the matrix elements $M(\tau, \mathbf{q}) = \int dr \varphi_e(\mathbf{r}-\tau)\varphi_h(\mathbf{r})e^{iq \cdot r}$ and remembering that $\varphi$ 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 $\varphi(\mathbf{r})$ are here Wannier functions with components on the neighboring sites,

$$\begin{aligned} \varphi_e(\mathbf{r}) &\approx \varphi_e^{\text{at}}(\mathbf{r}) - \frac{t_{\perp}}{2\Delta} \sum_{\tau} \varphi_h^{\text{at}}(\mathbf{r} + \tau) \text{ and} \\ \varphi_h(\mathbf{r}) &\approx \varphi_h^{\text{at}}(\mathbf{r}) + \frac{t_{\perp}}{2\Delta} \sum_{\tau} \varphi_e^{\text{at}}(\mathbf{r} - \tau), \end{aligned} \quad (C5)$$

where the $\varphi_j^{\text{at}}(\mathbf{r})$ are the genuine localized atomic orbitals (see Ref. [29] for a discussion). As a consequence in the macroscopic limit $\mathbf{q} \to 0$ we can write:

$$M(\tau, \mathbf{q}) \approx \frac{t_{\perp}}{2\Delta} (e^{i\mathbf{q} \cdot \tau} - 1), \quad (C6)$$

which clearly corresponds to dipolar integrals, i.e., overlap integrals weighted by $e^{\pm i\mathbf{q} \cdot \tau}$.

Finally, using (C6) and (C4) into (C2) one gets

$$\begin{aligned} \langle \tau', \mathbf{Q} | V | \tau, \mathbf{Q} \rangle = & \sum_G \frac{8\pi e^2}{\Omega_{\text{at}} |\mathbf{Q} + \mathbf{G}|^2} \left(\frac{t_{\perp}}{2\Delta}\right)^2 \\ & \times (e^{-i(\mathbf{Q} + \mathbf{G}) \cdot \tau} - 1)(e^{i(\mathbf{Q} + \mathbf{G}) \cdot \tau'} - 1), \end{aligned} \quad (C7)$$

where $\Omega_{\text{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 $\mathbf{k} \cdot \mathbf{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_{\perp}$, only the $G_{\perp} = 0$ term is singular in the $Q \to 0$ limit, and we finally obtain:

$$\langle \tau', \mathbf{Q} | V | \tau, \mathbf{Q} \rangle \approx J \frac{1}{4Q\tau} (e^{-i\mathbf{Q} \cdot \tau} - 1)(e^{i\mathbf{Q} \cdot \tau'} - 1) \quad (C8)$$

with $J = 4\pi e^2 t^2 / (3\sqrt{3}\tau\Delta^2)$ where the area of the 2D unit cell is equal to $3\sqrt{3}\tau^2$, $\tau = |\tau| = a/\sqrt{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 $\Phi^+$ and $\Phi^-$ as in the main text. Let $\theta$ be the polar angle of $Q$, then the matrix of $V$ in this $(\Phi^+, \Phi^-)$ space reads

$$V \propto \mathbf{Q} \cdot \mathbf{d}_{\Phi}(1 e^{2i\theta} e^{-2i\theta} 1), \quad (D1)$$

where $\mathbf{d}_{\Phi} = \sum_{\tau} \tau \langle \tau | \Phi \rangle$ is the dipolar matrix element of the $|\Phi\rangle$ state. The eigenvalues of the matrix above are equal to 0 and 2. The lower branch (eigenvalue 0) is a “transverse” eigenstate for which $\mathbf{Q} \cdot \mathbf{d}_{\Phi} = 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 $\tau$. 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 $\tau$ 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
<figure><img src="image_15.png" /><figcaption>FIG. 12. (a) The lattice points at the middle of the bonds of a honeycomb lattice form a Kagome lattice. Each point of this lattice represents an excitonic pair which can jump from one site to a neighboring one. (b) Band structure of the Kagome lattice. Energy is given in units of the hopping integral $|t|$.</figcaption></figure>

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.

<figure><img src="image_16.png" /><figcaption>FIG. 13. Exciton dispersion along $\Gamma K$ (full circles) and GW transitions (open circles) at two different interplane distances $c$ in the AB stacking. The unit $q_0 = 0.14$ Å$^{-1}$.</figcaption></figure>

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 $\Gamma$ 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$.
<figure><img src="image_17.png" /><figcaption>FIG. 14. Dispersion along $\Gamma K$ of the first exciton in AA' hBN computed with four different approximations. All curves have been aligned at $\Gamma$.</figcaption></figure>

## 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 \parallel \Gamma K$. All excitons have been aligned at $\Gamma$ 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.

[1] K. Watanabe, T. Taniguchi, T. Niiyama, K. Miya, and M. Taniguchi, Nat. Photonics **3**, 591 (2009).

[2] Z. Jacob, Nat. Mater. **13**, 1081 (2014).

[3] T. T. Tran, C. Elbadawi, D. Totonjian, C. J. Lobo, G. Grosso, H. Moon, D. R. Englund, M. J. Ford, I. Aharonovich, and M. Toth, ACS Nano **10**, 7331 (2016).

[4] R. Bourrellier, S. Meuret, A. Tararan, O. Stéphan, M. Kociak, L. H. G. Tizei, and A. Zobelli, Nano Lett. **16**, 4317 (2016).

[5] J. Wu, H. Wang, L. Jiang, J. Guo, X. Dai, Y. Xiang, and S. Wen, J. Appl. Phys. **119**, 203107 (2016).

[6] K. Watanabe, T. Taniguchi, T. Kuroda, and H. Kanda, Diam. Relat. Mater. **15**, 1891 (2006).

[7] P. Jaffrennou, J. Barjon, J.-S. Lauret, A. Loiseau, F. Ducastelle, and B. Attal-Tretout, J. Appl. Phys. **102**, 116102 (2007).

[8] K. Watanabe and T. Taniguchi, Phys. Rev. B **79**, 193104 (2009).

[9] L. Museur, G. Brasse, A. Pierret, S. Maine, B. Attal-Tretout, F. Ducastelle, A. Loiseau, J. Barjon, K. Watanabe, T. Taniguchi, and A. Kanaev, Phys. status solidi - Rapid Res. Lett. **5**, 214 (2011).

[10] A. Pierret, J. Loayza, B. Berini, A. Betz, B. Plaçais, F. Ducastelle, J. Barjon, and A. Loiseau, Phys. Rev. B **89**, 035414 (2014).

[11] R. Bourrellier, M. Amato, L. H. Galvão Tizei, C. Giorgetti, A. Gloter, M. I. Heggie, K. March, O. Stéphan, L. Reining, M. Kociak, and A. Zobelli, ACS Photonics **1**, 857 (2014).

[12] X. Z. Du, J. Li, J. Y. Lin, and H. X. Jiang, Appl. Phys. Lett. **106**, 021110 (2015).

[13] G. Cassabois, P. Valvin, and B. Gil, Nat. Photon. **10**, 262 (2016).

[14] G. Cassabois, P. Valvin, and B. Gil, Phys. Rev. B **93**, 035207 (2016).

[15] T. C. Doan, J. Li, J. Y. Lin, and H. X. Jiang, Appl. Phys. Lett. **109**, 122101 (2016).

[16] J. Li, X. K. Cao, T. B. Hoffman, J. H. Edgar, J. Y. Lin, and H. X. Jiang, Appl. Phys. Lett. **108**, 122101 (2016).

[17] H. Henck, D. Pierucci, G. Fugallo, J. Avila, G. Cassabois, Y. J. Dappe, M. G. Silly, C. Chen, B. Gil, M. Gatti, F. Sottile, F. Sirotti, M. C. Asensio, and A. Ouerghi, Phys. Rev. B **95**, 085410 (2017).

[18] L. Schué, B. Berini, A. C. Betz, B. Plaçais, F. Ducastelle, J. Barjon, and A. Loiseau, Nanoscale **8**, 6986 (2016).

[19] L. Schué, I. Stenger, F. Fossard, A. Loiseau, and J. Barjon, 2D Mater. **4**, 015028 (2017).

[20] T. Q. P. Vuong, G. Cassabois, P. Valvin, V. Jacques, R. Cuscó, L. Artús, and B. Gil, Phys. Rev. B **95**, 045207 (2017).

[21] X. Blase, A. Rubio, S. G. Louie, and M. L. Cohen, Phys. Rev. B **51**, 6868 (1995).

[22] B. Arnaud, S. Lebègue, P. Rabiller, and M. Alouani, Phys. Rev. Lett. **96**, 026402 (2006).

[23] B. Arnaud, S. Lebègue, P. Rabiller, and M. Alouani, Phys. Rev. Lett. **100**, 189702 (2008).

[24] L. Wirtz, A. Marini, and A. Rubio, Phys. Rev. Lett. **96**, 126104 (2006).

[25] L. Wirtz, A. Marini, M. Grüning, C. Attaccalite, G. Kresse, and A. Rubio, Phys. Rev. Lett. **100**, 189701 (2008).

[26] W. Aggoune, C. Cocchi, D. Nabok, K. Rezouali, M. A. Belkhir, and C. Draxl, Phys. Rev. B **97**, 241114 (2018).

[27] L. Schué, L. Sponza, A. Plaud, H. Bensalah, K. Watanabe, T. Taniguchi, F. Ducastelle, A. Loiseau, and J. Barjon, arXiv:1803.03766.
[28] R. M. Ribeiro and N. M. R. Peres, Phys. Rev. B **83**, 235312 (2011).

[29] T. Galvani, F. Paleari, H. P. C. Miranda, A. Molina-Sánchez, L. Wirtz, S. Latil, H. Amara, and F. Ducastelle, Phys. Rev. B **94**, 125303 (2016).

[30] J. Kang, L. Zhang, and S.-H. Wei, J. Phys. Chem. Lett. **7**, 597 (2016).

[31] F. Paleari, T. Galvani, H. Amara, F. Ducastelle, A. Molina-Sánchez, and L. Wirtz, 2D Mater. **5**, 045017 (2018).

[32] L. Sponza, H. Amara, F. Ducastelle, A. Loiseau, and C. Attaccalite, Phys. Rev. B **97**, 075121 (2018).

[33] M. Chubarov, H. Pedersen, H. Högberg, J. Jensen, and A. Henry, Crystal Growth & Design **12**, 3215 (2012).

[34] C.-J. Kim, L. Brown, M. W. Graham, R. Hovden, R. W. Havener, P. L. McEuen, D. A. Muller, and J. Park, Nano Lett. **13**, 5660 (2013).

[35] R. Martin, L. Reining, and D. Ceperley, Interacting Electrons (Cambridge University Press, New York, 2016).

[36] J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen, M. Dułak, L. Ferrighi, J. Gavnholt, C. Glinsvad, V. Haikola, H. A. Hansen, H. H. Kristoffersen, M. Kuisma, A. H. Larsen, L. Lehtovaara, M. Ljungberg, O. Lopez-Acevedo, P. G. Moses, J. Ojanen, T. Olsen, V. Petzold, N. A. Romero, J. Stausholm-Møller, M. Strange, G. A. Tritsaris, M. Vanin, M. Walter, B. Hammer, H. Häkkinen, G. K. H. Madsen, R. M. Nieminen, J. K. Nørskov, M. Puska, T. T. Rantala, J. Schiøtz, K. S. Thygesen, and K. W. Jacobsen, J. Phys.: Condens. Matter **22**, 253202 (2010).

[37] F. Fossard, L. Sponza, L. Schué, C. Attaccalite, F. Ducastelle, J. Barjon, and A. Loiseau, Phys. Rev. B **96**, 115304 (2017).

[38] X. Gonze, G. Rignanese, M. Verstraete, J. Betiken, Y. Pouillon, R. Caracas, F. Jollet, M. Torrent, G. Zerah, M. Mikami, P. Ghosez, M. Veithen, J.-Y. Raty, V. Olevano, F. Bruneval, L. Reining, R. Godby, G. Onida, D. Hamann, and D. Allan, Zeitschrift für Krsitallographie - Crystalline Materials **220**, 558 (2005).

[39] http://etsf.polytechnique.fr/exc/.

[40] In the single layer, the origin belongs to $\Lambda = \Lambda_N$ and $\mathbf{R} \in \Lambda_B = \Lambda + \tau_j$ with $j$ any of 1, 2, or 3. We will take $j = 3$ for simplicity.

[41] G. H. Wannier, Phys. Rev. **52**, 191 (1937).

[42] R. S. Knox, Theory of Excitons (Academic Press, New York, 1963).

[43] Y. Toyozawa, Optical Processes in Solids (Cambridge University Press, Cambridge, United Kingdom, 2003).

[44] F. Bechstedt, Many-Body Approach to Electronic Excitations (Springer, Berlin, Germany, 2015).

[45] In the cited work [31], a detailed explanation of the splitting of excitonic states is provided showing also that surface effects lead to an energetic separation between excitons localized on the inner layers and excitons localized on the outer layers.

[46] J. Koskelo, G. Fugallo, M. Hakala, M. Gatti, F. Sottile, and P. Cudazzo, Phys. Rev. B **95**, 035125 (2017).

[47] L. Wirtz and A. Rubio, Optical and vibrational properties of boron nitride nanotubes, in *B-C-N Nanotubes and Related Nanostructures* (Springer New York, New York, 2009), pp. 105–148.

[48] N. Berseneva, A. Gulans, A. V. Krasheninnikov, and R. M. Nieminen, Phys. Rev. B **87**, 035404 (2013).

[49] F. Hüser, T. Olsen, and K. S. Thygesen, Phys. Rev. B **87**, 235132 (2013).

[50] P. Cudazzo, L. Sponza, C. Giorgetti, L. Reining, F. Sottile, and M. Gatti, Phys. Rev. Lett. **116**, 066803 (2016).

[51] D. Y. Qiu, T. Cao, and S. G. Louie, Phys. Rev. Lett. **115**, 176801 (2015).

[52] F. Wu, F. Qu, and A. H. MacDonald, Phys. Rev. B **91**, 075310 (2015).

[53] L. Zheng, L. Feng, and W. Yong-Shi, Chin. Phys. B **23**, 077308 (2014).

[54] S. Galambosi, L. Wirtz, J. A. Soininen, J. Serrano, A. Marini, K. Watanabe, T. Taniguchi, S. Huotari, A. Rubio, and K. Hämäläinen, Phys. Rev. B **83**, 081413 (2011).

[55] R. Schuster, C. Habenicht, M. Ahmad, M. Knupfer, and B. Büchner, Phys. Rev. B **97**, 041201 (2018).

[56] C.-H. Ho, C.-P. Chang, and M.-F. Lin, Phys. Rev. B **93**, 075437 (2016).

[57] In principle all parameters should be adjusted; in particular the gap equal to $2\Delta$ should be reduced, but as far as orders of magnitude are concerned this will not modify our discussion.

[58] Similar analyses have been performed in the case of carbon layer stackings [61].

[59] A. Marini, C. Hogan, M. Grüning, and D. Varsano, Comput. Phys. Commun. **180**, 1392 (2009).

[60] C. Barreteau, F. Ducastelle, and T. Mallah, J. Phys.: Condens. Matter **29**, 465302 (2017).

[61] J.-C. Charlier, X. Gonze, and J.-P. Michenaud, Carbon **32**, 289 (1994).