Excitons in boron nitride single layer

# Excitons in boron nitride single layer

Thomas Galvani    Fulvio Paleari    Henrique P. C. Miranda    Alejandro Molina-Sánchez    Ludger Wirtz Physics and Materials Science Research Unit, University of Luxembourg, 162a avenue de la Faïencerie, L-1511 Luxembourg, Luxembourg, EU    Sylvain Latil CEA, IRAMIS, SPEC, GMT, 91191 Gif-sur-Yvette, France    Hakim Amara Laboratoire d’Etude des Microstructures, ONERA-CNRS, BP 72, 92322 Châtillon Cedex, France    François Ducastelle Laboratoire d’Etude des Microstructures, ONERA-CNRS, BP 72, 92322 Châtillon Cedex, France
July 12, 2019
###### Abstract

Boron nitride single layer belongs to the family of 2D materials whose optical properties are currently receiving considerable attention. Strong excitonic effects have already been observed in the bulk and still stronger effects are predicted for single layers. We present here a detailed study of these properties by combining ab initio calculations and a tight-binding-Wannier analysis in both real and reciprocal space. Due to the simplicity of the band structure with single valence () and conduction () bands the tight-binding analysis becomes quasi quantitative with only two adjustable parameters and provides tools for a detailed analysis of the exciton properties. Strong deviations from the usual hydrogenic model are evidenced. The ground state exciton is not a genuine Frenkel exciton, but a very localized “tightly-bound” one. The other ones are similar to those found in transition metal dichalcogenides and, although more localized, can be described within a Wannier-Mott scheme.

###### pacs:
71.20.Nr, 71.35.-y, 71.55.Eq, 78.20.Bh, 78.67.-n

## I Introduction

Two-dimensional materials are currently the object of many investigations concerning their electronic and optical properties. Graphene is the most known exampleCastro Neto et al. (2009) but hexagonal boron nitride,Watanabe et al. (2006); Jaffrennou et al. (2007); Watanabe and Taniguchi (2009); Museur et al. (2011); Pierret et al. (2014) and now transition metal dichalcogenides (TMD)Zhao et al. (2012); Xu et al. (2013); Zhang (2015); Molina-Sánchez et al. (2015); Wu et al. (2015) as well as new materials such as phosphorene,Ling et al. (2015); Favron et al. (2015) silicene, germanene, etc.Vogt et al. (2012); Li et al. (2014) are receiving considerable attention. Apart from graphene these materials are semiconductors which new optical properties compared to those of the familiar 3D semiconductors. Excitonic effects, in particular, are more pronounced in 2D than in 3D, with the exciton binding energies being much higher, of the order of 0.1 – 1 eV or more. The spatial extension of the excitons remains fairly large in general so that they are frequently considered as Wannier-Mott excitons. However, it has been quickly noticed that the usual hydrogenic model does not apply in 2D because of the different screening processes involved.Mak et al. (2010); Splendiani et al. (2010); Cudazzo et al. (2011); Chernikov et al. (2014); Huang et al. (2013); Molina-Sánchez et al. (2013); Qiu et al. (2013); He et al. (2014); Pulci et al. (2015) There is thus a need to understand more precisely these excitonic effects. The case of hexagonal boron nitride (hBN) is more specific still. Even in its bulk hexagonal form, very strong excitonic effects have been reported very early based on ab-initio calculations.Arnaud et al. (2006) The theoretical interpretation of a very strongly bound excitons (0.7 eV for the ground-state exciton in bulk hBN) was confirmed by various experimentWatanabe et al. (2006); Jaffrennou et al. (2007); Watanabe and Taniguchi (2009); Museur et al. (2011); Pierret et al. (2014) and was refined in theoretical calculations considering symmetry arguments.Wirtz et al. (2006, 2008); Arnaud et al. (2008) The reason for the strong binding energy is the quasi-2D nature of the hBN structureArnaud et al. (2006) consisting of stackings of hexagonal layers interacting through weak (mostly Van der Waals like) interactions. Furthermore hBN has a very large gap, 6 eV, leading to a rather weak dielectric screening, so that all ingredients conspire to enhance these excitonic effects. They have been studied recently, but the experiments are difficult because of the necessity to work in the far UV range.

The current interest in 2D materials and the development of techniques to handle few-layer materials suggest of course to study the properties of hBN as a fonction of the number of layers, as has been done in the case of graphene and TMD. In the case of TMD it has been shown that the nature of the gap, indirect or direct, depends on this number of layers. What about hBN ? Preliminary experimental studies are already available,Schue et al. (2016) but much remains to be done, and first, a precise knowledge of the single layer (SL) properties is required.

This is the main purpose of the present article. We present a detailed theoretical study of the first excitonic levels, and characterize their energies and shape by combining ab initio calculations and a simple tight-binding model. The ab initio approach is the usual one, based on a GW plus Bethe-Salpeter approach. The tight-binding approach is close to the approach put forward by Wannier long ago.Wannier (1937); Knox (1963); Toyozawa (2003); Bechstedt (2015) It turns out here that we have basically to take into account just one valence band and one conduction band. Furthermore, close to the gap, the corresponding Bloch states are concentrated on the nitrogen (N) and boron (B) atoms respectively, so that the usual orbitals can be considered as genuine Wannier functions. It is then possible to work out the Wannier equations in real space in a fairly simple but surprisingly accurate way.

The paper is organised as follows: Section II is devoted to the electronic structure of hBN-SL which is calculated using standard ab initio techniques and fitted to a simple tight-binding model. Section III contains the principal discussion of the various excitons and of their symmetry, using in particular imaging tools in real and reciprocal spaces. Finally Section IV is concerned with the calculation of optical matrix elements. This is followed by a discussion (Section V) and several appendices.

## Ii Electronic structure of hBN single layer

### ii.1 Band structure

We first specify a few notations. The structure of the hBN single layer is shown in Fig. 1.

The band structure of the hBN single layer is shown in Fig. 2. The Density Functional Theory (DFT) calculations have been made using the Quantum ESPRESSO code with the local density approximation (LDA) for the exchange-correlation functional. Giannozzi et al. (2009) The GW corrections were computed in the GW approximation, using the YAMBO codeMarini et al. (2009) with the plasmon-pole approximation for the frequency dependence of the dielectric function. These corrections have been applied to the last two valence bands and to the first four conduction ones. The lattice parameter has been fixed at the experimental lattice constant a.u. (2.50 Å). The computational details are the same used for the subsequent Bethe-Salpeter calculation and are found in Section III. The gap, equal to 7.25 eV, is direct between the and band at point in the Brillouin zone, while the bands are very flat along the lines. This is in agreement with several previous calculations.Blase et al. (1995); Wirtz and Rubio (2009); Ribeiro and Peres (2011); Berseneva et al. (2013); Hüser et al. (2013); Cudazzo et al. (2016) Notice however that the GW approximation is known to underestimate large band gaps. In addition, after the GW corrections are made, the bottom of the conduction band at the point is lower than at . This is due to the nearly-free electron statesBlase et al. (1995) that are forming around isolated layers of hBN: with increasing inter-layer distance a larger number of states appears. The transition matrix elements from localized valence band states into these states are very low, so for all practical purposes, the isolated sheet of hBN can be considered to be a direct band gap material. Regardless of the nature of the quasiparticle gap (direct or indirect), the optical gap is direct at point . The ab initio results confirm that the contributions to all exciton states of interest in this work come from transitions near and away from .

As in the case of graphene, these two bands can be reproduced fairly well using a simple tight-binding model. Let us denote the atomic state at site . The corresponding atomic orbital is . Then, we define the two Bloch functions on the and sublattices:

 |kA(B)⟩=1√N∑n∈A(B)eiK.n|n⟩,

where is the number of unit cells, i.e. half the number of atoms. As usual, in most cases we just keep first neighbour hopping integrals , and the nitrogen and boron atoms are distinguished by their on-site matrix elements, equal to on the sites for the N atoms, and to on the sites for the B atoms. The matrix elements of the hamiltonian in the Bloch basis are therefore written as:

 ⟨kA|H|kA⟩=−Δ⟨kB|H|kB⟩=+Δ⟨kA|H|kB⟩=⟨kB|H|kA⟩∗=−tγ(k)γ(k)=∑α=1,2,3eik.τα. (1)

The energy eigenvalues are then given by:

 E=sEk;Ek=√Δ2+t2|γ(k)|2;s=sgn(E),

and the eigenstates are:

 |ks⟩=CAs|kA⟩+CBs|kB⟩.

Finally, up to a phase factor the coefficients et are given by:

 CAs=⟨kA|ks⟩=−sγ(k)|γ(k)|√Ek−sΔ2EkCBs=⟨kB|ks⟩=√Ek+sΔ2Ek. (2)

Thus the electronic structure of hBN-SL can be characterized by only two parameters and . Their order of magnitude is eV but more precise values can be obtained by fitting the valence and conduction bands to those provided by ab initio calculations. is fixed so that the gap is equal to the ab initio one, eV, and is then obtained using standard fitting procedures. Different values are obtained depending on the energy range where the fit is optimized. A global fit, disregarding the nearly-free electron states, leads to eV, but here we are more interested to have a good fit along the MK line, in which case we obtain eV (similar to the values for recent fitting of band structures in Ref. Ribeiro and Peres, 2011).

As usual, and as shown in Fig. 3 the fit is better for the valence band than for the conduction band. The fit can easily be improved by adding further neighbour interactions. Neighbours on the same sublattice contribute to diagonal matrix elements whereas neighbours on different sublattices contribute to off-diagonal elements.

Before considering second neighbour interactions explicitly, let us see a simple way to deduce the band structure of hBN-SL from that of graphene. Let be the hamiltonian of graphene with only hopping integrals. The full hamiltonian of hBN-SL is given by where is the “atomic” diagonal hamiltonian with matrix elements equal to on the sublattice and to on the sublattice. It is clear that so that where is a simple constant (multiplied by the unit matrix). on the other hand is an hamiltonian connecting sites entirely on sublattice or on sublattice . On these triangular lattices connects the first neighbours with a hopping integral equal to but it has also diagonal on-site matrix elements equal to . In the Bloch basis is therefore diagonal:

 (H0)2=∑k|kA⟩|tγ(k)|2⟨kA|+|kB⟩|tγ(k)|2⟨kB|,

and the eigenvalues are indeed equal to .

Adding second neighbour interactions in the hBN-SL structure is then equivalent to introducing first neighbour interactions on the triangular sublattices, and the eigenvalues are therefore given by:

 Ek=−t2(|γ(k)|2−3)±√Δ2+t2|γ(k)|2.

Since both and are positive, second neighbour interactions induce an asymmetry between the valence and the conduction band: The conduction band becomes flatter than the valence band, in agreement with ab initio calculations (Fig. 3). The best local fit along is provided by eV ; eV. Actually, under the approximation that the valence and conduction bands are pure and states, as shown below, only the energy difference between these two bands enters the tight binding excitonic hamilitonian derived in section III, and the second nearest neighbours hopping term does not contribute to this difference. For this reason we limit our tight-binding model for the excitons to first nearest neighbour hopping and keep the simplest previous fit with eV, eV.

### ii.2 Wave functions, densities of states

Many electronic properties of hBN-SL only depend on the electronic states close to the gap, i.e. in energy ranges where is small compared to the gap . This means that in a first approximation, the coefficients are equal to one or zero. In other words close to the gap, the valence states are concentrated on the N sites whereas the conduction states are concentrated on the B sites, and the eigenvalues can be approximated by:

 Ek≃±(Δ+t22Δ|γ(k)|2).

To examine the validity of this approximation, we have calculated the local densities of states on both N and B sites. They are shown in Fig. ï¿½4. In our simple tight-binding model, one can easily see that . Furthermore shows a step-function-like onset at , whereas has its onset at . As a result the states are indeed quasi-pure B states in a fairly broad energy range above . And of course they are quasi-pure N states below .

To summarize we can assume that the Wannier functions associated with the valence band and the conduction band can be identified, to lowest order, with the atomic functions centred on the corresponding sites of their triangular lattices. The effective hamiltonians for the and for the are then given by:

 Hv =−∑n|nA⟩(Δ+3tv)⟨nA|−∑n,m′|nA⟩tv⟨mA| (3) Hc =∑n|nB⟩(Δ+3tc)⟨nB|+∑n,m′|nB⟩tc⟨mB|, (4)

where the primes indicate sums over nearest neighbours on the triangular lattices, and .

## Iii Excitons in hBN single layer

Ab initio excitonic calculations are based on the Bethe-Salpeter formalism,Rohlfing and Louie (2000); Onida et al. (2002); Bechstedt (2015) which in practice leads to an effective Schrödinger or Wannier equation for electron-hole pairs: 111Standard treatments of excitons can be found for example in [Knox, 1963],[Toyozawa, 2003], and [Bechstedt, 2015].

 (Ekc−Ekv)Φkvc+∑k′v′c′⟨kvc|Keh|k′v′c′⟩Φk′v′c′=E Φkvc,

where and are the conduction and valence band energy, respectively, is the electron-hole interaction kernel and is the electron-hole wave function in space. In this paper, we only consider vertical excitations where the electron and the hole have the same wave vector , i.e. we consider excitons with vanishing wave vector of their centre of mass. The are the coefficients in the expansion of the excitonic state in terms of electron-hole excitations:

 |Φ⟩=∑kΦkvca+ckavk|∅⟩,

where the vacuum state is the state where, at zero temperature, all valence states are full and all conduction states are empty. Only singlet states are considered here so that spin indices are omitted.

The Bethe-Salpeter equation has been solved using the YAMBO code. Marini et al. (2009) A Coulomb cutoff of the screened potential in the vertical direction has been used in order to avoid long-range interaction between repeated copies of the monolayer. Rozzi et al. (2006) In this way, we find that both the GW corrections and the first excitonic peaks are already converged (with about eV accuracy) with an inter-layer separation of atomic units. The same level of convergence was achieved by sampling the two-dimensional Brillouin zone with a -point grid. The internal YAMBO parameters for many-body perturbation theory (MBPT) calculations were carefully converged as well. We also carried out calculations with a -point grid in order to show the higher-level excitonic wavefunctions in real space without any overlap between repeated copies on the same monolayer. We verified furthermore that choosing an inter-layer separation of a.u. does not modify the results.

### iii.1 Ground state exciton

As an introduction, we present ab initio results concerning the ground state exciton level. Its binding energy measured with respect to the bottom of the conduction band is huge: 1.9 eV. In Fig. 5 we show an image of the excitonic wave function , where the hole (at ) is localized just above a nitrogen atom. The plot represents the total probability , i.e. the probability to find the electron at position if the hole is located at . Since this exciton is doubly degenerate, we sum the total probability over the two degenerate states in order to preserve the trigonal symmetry of the crystal lattice.Wirtz et al. (2008) As expected, the electron density is centred on the boron atoms, with a high probability — about 30% — on the first nearest neighbours. Although not in the genuine Frenkel limit, the exciton is well localized in real space. The shape of this exciton is actually quite similar to that found for the 3D hBN crystal, which is not surprising since in the latter case the exciton was already found to be well confined in a single layer although with a lower binding energy.Arnaud et al. (2006); Wirtz et al. (2006, 2008); Wirtz and Rubio (2009) We also show the wave function in reciprocal space. Here we plot the (summed) weight of the electron-hole pairs of wave vector that constitute the bound exciton. The distribution is peaked around the high-symmetry point but extends toward the boundaries of the Brillouin zone, i.e. along the lines.

### iii.2 Wannier tight-binding model

We use now the result derived in Section II that the valence and conduction states are pure (N atoms) and (B atoms) states, respectively. This approximation is fully justified when looking at the ab initio results of the previous section. Then we can write:

 a+kcakv|∅⟩≃a+kBakA|∅⟩≃1N∑n,ma+mBanAeik.(m−n)|∅⟩.

The sum over the electron and hole positions can be decomposed into a sum over the hole position and over the electron-hole distance which is then a vector joining a site on the (hole) sublattice to a (electron) sublattice site. For simplicity we use “bra” and “ket” notations:

 |kvc⟩ =a+kcakv|∅⟩=1√N∑Reik.R|Rvc⟩ (5) |Rvc⟩ =1√N∑na+nA+RanA|∅⟩. (6)

We see that is the linear superposition of exciton amplitudes for pairs separated by . This is nothing but the Bloch wave function for excitons with a wave vector . We can then rewrite Bethe-Salpeter-Wannier equation in real space using the basis and the wave function coefficients . The “kinetic energy” term , diagonal in space, becomes a tight-binding-like term in space equal to with . To be more precise let be the hamiltonian acting in the excitonic space. The kinetic energy (or free single-particle) term is the difference of the two hamiltonians (3) and (4). We rewrite it here to indicate that each term in the r.h.s. acts either on the electronic or on the hole component of the electron-hole states. Dropping now the indices within the “bras” and “kets”, the matrix elements of are equal to and are therefore equal to if , and to if and are first neighbours on the triangular lattice.

The electron-hole interaction kernel on the other hand contains a direct term and an exchange contribution. Let us first consider the direct term which is the most important one, as will be checked later (see Table 1). The corresponding integral can be expanded in real space and involve integrals of type:

 −∫drdr′φv(r−Rn)φv(r−Rp)e2|r−r′|φc(r′−Rm)φc(r′−Rq),

where the are (real) valence or conduction orbitals.

Usually the largest integral is the one where all indices are identical, but here this on-site integral is forbidden since the conduction and valence orbitals belong to different sublattices (actually, this is not completely true as will be discussed below in Section III.3). The next most important integrals are those where and , and finally the (direct) Coulomb term becomes where:

 UR=⟨R|Kdeh|R⟩=−∫drdr′φ2c(r)e2|r−r′|φ2v(r′−R),

which means that acts as a local potential on “site” . The Coulomb potential should also be screened but here the will just be considered as parameters to be fitted to ab initio data.

We have therefore reduced our problem to a very simple tight-binding problem for the relative motion of the electron and of the hole. Since the motion is relative we can fix the hole at the origin of the sublattice. The vectors lie on the sublattice: our problem becomes the problem of an electron moving on the sublattice in the presence of a hole at the origin which plays the part of an impurity, source of the attractive potential . To summarize, when exchange effects are neglected, we have to handle the standard tight-binding equations:

 EΦR=∑R′heh(R−R′)ΦR′+∑RURΦR,

which therefore depends only on and on .

### iii.3 Discussion of the Wannier model

Although standard, the Wannier equations are difficult to solve in many cases because several valence and conduction bands are involved. As a consequence the Wannier functions have no longer direct relationships with the atomic orbitals. On the other hand, in the case of strong screening and small gaps, the potential does not perturb the single particle Bloch states too much and expansions can be used, leading to the familiar hydrogenic model where the underlying lattice can finally be forgotten. This model is a posteriori justified when the extension of the excitonic states is large compared to the lattice parameter. As shown in Fig. 5 this is not the case here, but the simplicity of the electronic structure of the boron nitride single sheet will allow us to take lattice effects fully in account.

#### iii.3.1 A very simple model

As shown above, we have to solve an impurity problem in a simple tight-binding basis. In the case of localized potentials, Green function or direct diagonalization techniques are known to be very efficient. Actually we have just to transpose the methods used to study deep impurities centres in semiconductors.Lannoo and Bourgoin (1981); Yu and Cardona (2010) Although the Coulomb potential is a long range potential we expect that the energy of the lowest bound state can reasonably be obtained using truncated potentials in real space. Let us recall that the electron is moving on a triangular lattice with first neighbour hopping integrals in the presence of an impurity located at the origin taken at a lattice point of the hole sublattice, i.e. at the centre of a triangle of the B sublattice (see Fig. 6). We have therefore to diagonalize the following hamiltonian:

 Heh =H0eh+U ⟨R|H0eh|R′⟩ =2Δ+3t2/ΔifR=R′ =t2/ΔifRandR′are first neighbours =0otherwise U =∑R|R⟩UR⟨R|,

where the potential is to be fitted to ab initio data, and where exchange terms are neglected. The spectrum of is known since it is the spectrum of the triangular lattice with (positive) first neighbour hopping integrals (Fig. 6). The lowest eigenvalue (taken as the origin in Fig. 6) is at , which corresponds to the energy gap in this model. In the presence of a localized attractive potential, we expect that bound excitonic states appear when the potential is strong enough. This is indeed what happens. A brief analytical discussion is presented in Appendix A. Here, we present results obtained with a potential fitted up to the 28th neighbouring shell to ab initio data, and which is discussed below, in Sec. III.4. We have diagonalized our tight-binding hamiltonian using a box containing about sites on the triangular lattice. Many exciton states have been studied, but since our model is obviously becoming inaccurate when the energy rises, and when the extension of the exciton increases, five states are just considered in detail. The advantage of our procedure in real space is that we handle real wave functions and so we can simply image the wave function themselves. Furthermore in the case of degenerate states it is easy to show components of definite symmetry.

#### iii.3.2 The ground state exciton and the exciton symmetries

We show first in Fig. 7 the results concerning the ground state exciton which is doubly degenerate. In the tight-binding case, we show the two components which are clearly antisymmetric or symmetric with respect to the -axis. The agreement with the ab initio result is very good. At this point it is useful to comment on the symmetry of this state. Since we have fixed the position of the hole, we can use the point symmetry of the triangular lattice with respect to the origin located at a centre of a triangle. In principle the problem is not purely a 2D one since the orbitals extend in the direction and are odd with respect to a reflection. Apart from this trivial symmetry, we have just to consider the symmetry with its 3-fold rotation axis and its mirror planes . This group is known to have three different representations. Beyond the identity one, we have the familiar 2-dimensional representation , and a second representation of dimension 1, characterized by an odd character for the reflections. Our exciton has clearly here the symmetry with two (chiral) components which can be taken as varying as or .

This exciton is very similar to the so-called or excitons met in TMD. In this case, the Wannier-Mott approach in the approximation is generally used, and the symmetry of the excitons is frequently defined as follows: The exciton wave function is written in the form , where the are the single particle Bloch functions at point corresponding to the considered direct gap, and , the envelope function, is the solution of the hydrogenic-like excitonic equation for the relative coordinate .Jorio et al. (2011) The full excitonic symmetry is the symmetry of this product, but notations generally use the symmetry of . This decomposition makes sense uniquely if the expansion around is valid, which is not necessarily the case here. In our case, the direct gap occurs at points and . Neglecting intervalley coupling (which may not be valid either), we see that and the product of Bloch functions varies as , i.e. as one component of the representation of the symmetry at point . Now, since the conduction and valence bands are non degenerate at point , the ground state envelope function is isotropic and of symmetry . This is why the exciton is denoted a exciton. With similar arguments we obtain that the exciton at has also a symmetry modulated by the Bloch function proportional to . So, in this description we obtain two degenerate excitons, but if they are considered together they form a double degenerate exciton of symmetry . Both descriptions are equivalent in the case of large excitons which can be associated separately to points and .Wu et al. (2015) In our case where the exciton is localized in real space (delocalized in reciprocal space), using directly the full point symmetry of the exciton is more accurate.

#### iii.3.3 Other excitons

##### Analysis in real space

At higher energy, a group of six states appears. All of them as well as the previous state have similar energies within 0.1 eV. Actually they do not appear in the same order in both ab initio and tight-binding calculations. Their wave functions are on the other hand very similar. We follow here the order provided by the ab initio calculations. Both ab initio and tight-binding approaches find first a similar exciton with again a twofold degeneracy (exciton #2 in Fig. 8). It has therefore also an symmetry. The agreement between both calculations is still fairly good. The two following ones are non degenerate. The TB method shows unambiguously that the first one transforms according to the representation (exciton #3, Fig. 8) and the second one according to the one (exciton #4, Fig. 8). Finally, the two upper states are degenerate and belong to the symmetry (fifth exciton in Table 1). We summarize in Table 1 the energies and symmetries of these excitons. The next excitons are found more than 0.2 eV above this group in the ab initio calculations.

Although the overall agreement between ab initio and TB calculations is fairly good, the behaviour of the exciton seems particular. This is still more obvious if we compare ab initio calculations performed with and without the exchange contribution (Table 1). Whereas the energy variation between the two calculations for the other excitons is of a few percents, the exciton is strongly perturbed, its energy moving from eV to eV when the (repulsive) exchange contribution is suppressed. This is quite unusual but can be related to the fact that the exciton is, by symmetry, the only one where the electronic intensity at the origin is non-vanishing. Actually, although we have neglected this possibility in our simplified TB model, the ab initio calculations do show such a non-vanishing intensity (Fig. 8). The point is that even if this intensity is low, it introduces a perturbation proportional to the on-site exchange term , which is very large. Actually, the local Coulomb and exchange integrals and are of the same order of magnitude since they both involve similar orbitals. They are of opposite signs however, which explains why our TB scheme which neglects on-site Coulomb interactions is not too bad even in this case.

##### 2s and 2p states: Analysis in reciprocal space

At this point it is instructive to compare our results with those obtained for dichalcogenides. In TMD the Wannier-Mott model is valid provided appropriate anisotropic potentials are used. It is then convenient, as discussed earlier, to analyse the excitons in each valley in terms of symmetries. The usual sequence is a level, and then a level nearly degenerate with a level. This level gives rise to four states because of the valley degeneracy. A careful examination of the symmetry of the and states close to the points based on the so-called massive Dirac model has shown that actually the degeneracy within each valley is lifted, but time reversal symmetry between and insures that the states are splitted into two doubly degenerate states. Furthermore the level is found to be above the levels. Wu et al. (2015); Srivastava and Imamoğlu (2015); Zhou et al. (2015); Berkelbach et al. (2015) In our case where lattice effects and therefore inter-valley effects are included a further splitting occurs. As pointed out above, the full symmetry of the exciton states is obtained from the product of the envelope function ( symmetry for states) and of the Bloch functions at points of symmetry also. Now, the decomposition of the representation gives rise precisely to the observed one: . In the language, the states are first splitted into states whose chiralities are equal or opposite to those of the Bloch functions at points and . Hence two states (one in each valley) have a vanishing global chirality. Forming bonding and antibonding states between these states lead to the and states. The two other ones remain degenerate and form a state. We expect the bonding state to be below , but, as argued before, the repulsive exchange contribution neglected in this discussion pushes the upwards. This is described in Fig. 9.

It remains to determine which exciton belongs to this family. This might be exciton #2, as shown in the figure, or exciton #5. It is not obvious to decide from the plots of the wave functions in real space, but we show now that an analysis in reciprocal space provides the answer. In the first column of Fig. 10 are shown the intensities of the five excitons considered previously. It is clear at once that excitons #2, #3, and #4 belong to the same family and are therefore the expected excitons, which means that the relevant exciton is exciton #2 as depicted in Fig. 9. A consequence is that exciton #5 is the exciton. Furthermore we check that the state which is more spread out around the origin in real space than the state is more concentrated on the points in reciprocal space.

Tight-binding calculations lead to quite similar results as can be seen in the second column of Fig. 10. The advantage of the TB method is that we can easily obtain the wave functions themselves. It is possible also to have a modulus-phase representation by placing on each point of a grid a circle with color related to the phase between and with an opacity proportional to the intensity at this k-point, which is shown in the last columns of Fig. 10. The four states have more rich structures with phases rotating within each triangular spot located at the points. They can be explained if we use the model recalled above where the symmetry of the exciton states is governed by the product of the Bloch functions and an envelope function so that the exciton wave function is proportional to where here is one particular vector among the three equivalent ones. In the discrete tight-binding model is the sum of a vector of the triangular lattice and of any first neighbour vector . Assume now that is a envelope function and only depends on the modulus of . Under a rotation of angle , is transformed into another equivalent vector modulo a vector of the reciprocal lattice. The exciton state is therefore multiplied by a phase factor equal to or , where is the cubic root of unity, depending on the initial orientation of the lattice with respect to the origin. As mentioned previously, the wave function transforms as the component of positive chirality of the representation , and the wave function is given, up to a constant by:

 Φk∝∫dre−i(K1−k).rg(r)n(r),

where is the site density, i.e. the sum of Dirac functions on the triangular lattice sites, proportional to .

For close to , , we can neglect the variation of , i.e. keep only the term, and the integral over the angle yields , where is the Bessel function of zero order. If is peaked at some average value (remember that is at least equal to the minimum hole-electron distance ), and is therefore peaked at and the intensity is maximum in a circle of radius . If is close to another vector, say , then the integral is multiplied by a factor , where , i.e by a factor or . Actually the functions in reciprocal space have he same symmetry properties than in real space.

Assume now that has a symmetry, so that , where is the angle between and the -axis. Here also, the integral over this angle can be performed explicitly, so that , where is the Bessel function of order one, and is now the angle of with the -axis. We conclude that the phase rotates within each circle centred on the points. This is clearly as shown in Fig. 10. This proves definitely that the states #2 to #4 are of symmetry . More precisely consider first the non degenerate exciton #3 and #4. Exciton does show the symmetry already evidenced in real space. Furthermore one can notice that the phases rotates in opposite directions around the and points, and that these rotations are counterbalanced by the rotations between different points of the same family ( or ), in full agreement with the arguments put forward above. The same is true for exciton #3 except that the amplitudes are odd with respect to the mirrors (phase shift of ). One can also notice the signature of the Bessel functions: The intensities vanish at the origin of the spots, disappear at larger distances than in the case of states, and are maximum in between. Actually warping effects along the lines transform the circles into triangles. The case of the degenerate state is more complex, since each component seems to mix the behaviours of and states. Mixing between different states is allowed indeed (this is also the case for the exciton). It is clear however that the main features correspond to rotating phases within circles, and one can check that here the rotations within the circles and those between the circles are in the same direction.

To summarize, although the first excitons considered here are fairly localized, with important lattice effects, they can be classified to some extent within a scheme borrowed from the 2D atomic terminology ( states) and already applied successfully to TMD. The genuine symmetry of the exciton states is however more precisely related to the representations of the triangular point group. We have also seen that exchange effects are unusually strong for fully invariant states.

### iii.4 Fit of the potential

In a first approach we have tested a screened Coulomb potential, but it was quickly apparent that it was not possible in this way to reproduce accurately more than the first exciton. Meanwhile several developments in the literature were convincingly arguing that actually it is not possible in 2D to use such a potential and that a genuine 2D electrostatic potentialKeldysh (1979) should be used instead.Cudazzo et al. (2011); Berkelbach et al. (2013); Chernikov et al. (2014); Latini et al. (2015); Rodin et al. (2014); Wu et al. (2015); Pulci et al. (2015) We have therefore used the Keldysh potential: 222More precisely, we have used the simplified form proposed in [Cudazzo et al., 2011].

 V2D(r)=πe22r0[H0(rr0)−Y0(rr0)],

where the only parameter is the “screening” length , directly related to the 2D polarisability. Finally, since we are dealing with relative binding energies, our model only depends on two parameters, the excitonic hopping integral and . Of course, the Keldysh potential is still defined within a continuous approach which has no reason to apply exactly here where lattice effects are important. In the best fit, the hopping integral is found equal to 1.50 eV, so that eV, which is completely consistent with our value eV deduced from ab initio band structures. Finally we find Å. When is much larger than the potential tends to an unscreened Coulomb potential. Below the potential is screened and becomes logarithmic. Since the first electronic shell around the hole is at a first neighbour B–N distance, about 1.45 Å, we see in Fig. 11 that we are here in the screened regime where the potential is slowly varying. In Fig. 11 we show the effective distance dependent dielectric constant defined from . Notice also that the distance between the two first excitons is completely different from that predicted by the 2D hydrogenic model. This is due in part to these screening effects, but also and more importantly to lattice effects coupled with the specific electronic structure of hBN, since such deviations have already been observed with a fixed dielectric constant (see also Ref. [Wu et al., 2015]). The main reason is probably that it is nearly forbidden for the hole and the electron to be at the same position, which penalizes principally the binding energy of the ground state exciton.

We have neglected the exchange term here. Actually its short range part contributes generally to the splitting between the spin singlet and triplet states, the latter being dark in principle in the absence of spin-orbit coupling. This is a repulsive (positive effect) absent in the triplet term whose level should therefore be below the singlet one. Ab initio calculations predict a splitting about meV for the exciton.Wirtz et al. (2005) In a first order perturbation calculation this splitting is equal to the average of in the considered excitonic state. In the case of the ground state exciton, the excitonic wave function is concentrated on the first neighbour shell, so that this splitting is equal to a fraction of . As expected then eV . This is also completely consistent with the variation shown in Table 1 of the exciton energy when suppressing the exchange term. The particular case of the exciton which is very dependent on the intra-atomic values and of the Coulomb and exchange potentials has been discussed above. A simple perturbation method improving the simplest TB model used here can be derived to discuss this effect in more detail and is described in Appendix C.

To summarize, the simple tight-binding model for the excitons in hBN-SL is remarkably successful, even at a quantitative level and the comparison with ab initio calculations shows that the effective screened Coulomb potential to be used in a continuous model is really a potential of the Keldysh type in its strongly screened regime.

## Iv Optical matrix elements

The optical absorption is related to transitions from the ground state (energy ) to final states of energy and therefore to the corresponding matrix elements of the perturbation induced by the electromagnetic field. Each transition is then characterized by an oscillator strength . Here is the (unit) vector of the light polarization, and is the velocity operator.

### iv.1 Matrix elements between single particle states

In the absence of excitonic effets we have just to calculate the matrix elements between valence and conduction Bloch states with identical vectors. It is not difficult to calculate them in the general case,Margulis et al. (2012, 2014); Berkelbach et al. (2015) but here we just detail the calculation for states close to the gap where we know that the Bloch functions “live” on separate triangular sublattices. Then:

 ⟨kv|v|kc⟩=1N∑n,meik(m−n)⟨nA|v|mB⟩.

There is no unique way of calculating the matrix elements of , depending on whether we express it using the momentum operatorJorio et al. (2011) or the relation . Both methods are equivalent in an exact treatment but not when using an incomplete basis as is the case of our tight-binding basis. The second method has the disadvantage to use the operator which is not always well defined in periodic systems, but here there is no problem and the advantage is to work directly in real space, assuming that , and therefore:

 ⟨nA|v|mB⟩=−1iℏ(m−n)t, (7)

if and are first neighbours (on the honeycomb lattice), and zero otherwise. Then:

 ⟨kv|v|kc⟩ =1N∑n,meik(m−n)itℏ(m−n)=itℏ∑αeik.τατα =itℏ∇kγ(k).

In the limit , one finds , where is here the nearest neighbour distance, i.e. the lattice parameter divided by and and are the unit vectors along the -axis and the -axis, respectively, so that finally:

 ⟨Kv|v.^e|Kc⟩|≃v|ex+iey|;ℏv=3at/2.

Notice that is nothing but the effective mass of the conduction and valence bands at point . Using eV and km/s (as the Fermi velocity of graphene precisely given by ), we obtain .

At point , is replaced by . The matrix elements are maximum for circular polarized light and, for linearly polarized light, the matrix element is constant and equal to . The oscillator strength is equal to when is larger than the gap, i.e. about 2 close to the gap (equal to ). The absorption, proportional to , is therefore proportional to the density of states of the triangular lattice divided by . Its shape is characterized by a discontinuity at the edge and a Van Hove singularity at a distance equal to above the edge, as shown in Fig. 6.

### iv.2 Matrix elements between excitonic states

In the presence of excitons we have now to calculate the matrix element , where is the exciton state. From (6) and (7), we see that:

 v|∅⟩=−tiℏ∑n,m′(m−n)a+mBanA|∅⟩=itℏ√N∑ατα|τα⟩,

so that:

 ⟨∅|v.^e|Φ⟩=−itℏ√N∑α^e.τα⟨τα|Φ⟩

Defining the dipole associated with the exciton through:

 dΦ=∑ατα⟨τα|Φ⟩,

we see that . Thus only the local components of the exciton wave fonction contribute to the optical matrix element. This is completely equivalent to the statement that, within the usual hydrogenic model, only states contribute (Elliott theory, see [Yu and Cardona, 2010; Toyozawa, 2003]). Here we have an equivalent selection rule: the dipole should not vanish; in particular the wave function should have finite components on the first neighbours .

### iv.3 Application to the five first excitons of hBN-SL

Consider first the ground state exciton. It has two components and which can be chosen as those corresponding to circular polarizations, so that the components are the cubic roots of unity, , and , and finally . Here, is the amplitude of the exciton state on the first neighbours, at most equal to . In the case of single particle transitions the oscillator strength was of the order of for a transition close to the gap; hence a total oscillator strength of the order of times this value. We see here that the oscillator strength of the exciton is of the same order of magnitude if is large, i.e. if the exciton is strongly localized, which is the case here. Actually from ab initio calculations, the weight is found about one third its maximum value 1/3. In other words 30% of the weight of the ground state is concentrated on the first triangular shell. TB calculations on the other hand find a weight about 50%. The oscillator strength of the other excitons are smaller. The second exciton as well as the last one (#5) has the same symmetry as the first one but their amplitude on the first neighbours is weak. Finally the two other ones studied above (#3-4) are dark because their symmetry are characterized by representations and different from the vectorial representation , so that and this is confirmed by the ab initio calculations. Finally the ground state exciton takes almost all the oscillator strength.

## V Discussion

The excitons of hBN-SL have been characterized in detail. The first one, of lowest energy is particularly localized. Is it a Frenkel or Wannier-Mott exciton ? This discussion is somewhat semantic. It is in some sense similar to the long standing debate between the Heitler-London (atomic) approach and the Hund-Mulliken (“molecular”) approach to single particle properties. In practice, it turns out that in solids the band Hund-Mulliken approach is more fruitful since it can deal with many situations except when correlations effects are very strong. Even then, specific approaches “à la Hubbard” can be used and compete with the methods of quantum chemistry (interaction configuration approach). In between, the tight-binding method has proven very efficient to deal with electrons sharing itinerant properties (conductivity) and localized ones (magnetism, chemical bonding). We are certainly here in a similar situation. The localized excitons of hBN-SL can be described within a TB-Wannier framework, but cannot be described accurately within a approach similar to the nearly free electron approach of electronic properties. They could simply be described as “tightly-bound excitons”. We have shown indeed that in the case of hBN-SL which is a genuine case study, the tight-binding approach can be very accurate by fitting to ab initio data a few parameters.

On the experimental side optical properties of hBN-SL are not available yet, but there are already some indications that the expected main exciton is observed. In the case of bulk hBN, stacking effects induce splittings of this exciton level which are observed. Progress in the analysis of these stacking effects are under progress. Finally dispersion effects as well as exciton-phonon coupling remain to be studied.

Note added. A recent paper presents a model for excitons in dichalcogenides whose spirit is quite similar to our tight-binding model.Gunlycke and Tseng (2016)

###### Acknowledgements.
A. Molina-Sánchez and L. Wirtz acknowledge support from the National Research Fund, Luxembourg (Projects C14/MS/773152/FAST-2DMAT and INTER/ANR/13/20/NANOTMD). F. Ducastelle and H. Amara are indebted to L. Schué, J. Barjon, and A. Loiseau for numerous fruitful discussions. The research leading to these results has received funding from the European Union Seventh Framework Programme under grant agreement no. 604391 Graphene Flagship. We acknowledge funding by the French National Research Agency through Project No. ANR-14-CE08-0018.

## Appendix A A very simplified model for the ground state exciton

The ground state exciton is so localized that its properties do not depend too much on the long range part of the potential. It is useful therefore to examine the properties of a model where the range of the potential is limited to the three first neighbours of the central hole. We have then to diagonalize the following hamiltonian:

 Heh =H0eh+U ⟨R|H0eh|R′⟩ =3texcifR=R′ =texcifRetR′are first neighbours =0otherwise, U =∑R=1,2,3|R⟩u⟨R|,

where the three sites surrounding the hole at the origin are labelled .

### a.1 Green functions

The resolvant or Green function corresponding to this hamiltonian is . is then the Green function of the triangular lattice. With the chosen origin of energies, the spectrum of starts at with a jump equal to . In the presence of the attractive potential we can have a bound state if the determinant of the operator within the space of dimension 3 generated by the three states and vanishes. Let now and be the diagonal and off-diagonal matrix elements of , respectively ; . We find a double solution and a simple solution . Using standard methods one can determine the behaviour of and close to the origin. First, assuming a constant density of states , we find that:

 F0=1Wlogzz−W,

so that when is close to 0. Similarly, is found to behave as . Then is negative and diverges logarithmically below whereas tends to a constant. As a result the first equation has always a negative solution for , solution such that . With the previous model for and , we obtain, for small values of , If is large, . The corresponding eigenstates are, as expected, the “chiral”states and , where is the cubic root of unity, . We recover our exciton of symmetry . The components of beyond the (1,2,3) cluster can be obtained from the equation , i.e. .

### a.2 Reciprocal space

We can also express in reciprocal space:

 Φ±k =⟨k|ϕ±⟩=1√N∑Re−ik.R⟨R|ϕ±⟩ ≃1√3Nγ(±K−k)),

where we have limited the sum to the first neighbours and taken into account that when provided is chosen along the -axis, as in Fig. 1. Thus, up to a normalization constant the weight is equal to . Since is maximum when , we see that and are peaked at points and , respectively. As expected the sum is maximum on the boundary of the Brillouin zone, as shown in Fig. 12.

## Appendix B Exchange contribution

The exchange contribution involves integrals of type:

 +∫drdr′φv(r−Rn)φv(r′−Rp)2|r−r′|φc(r−Rm)φc(r′−Rq).

The largest integrals correspond to cases where the overlap is minimum for and integrations. But since we have forbidden site coincidence for valence and conduction orbitals, the best we can do is to consider first neighbour overlap between and and between and . Therefore we only keep the terms and , where and are first neighbours on the honeycomb lattice. The largest terms occur when , and the integral becomes a function of , and where measures the separation between the pairs and

 J(τ,τ′ρ)=∫drdr′φv(r)φc(r−τ)2|r−r′|φv(r′−ρ)φc(r′−ρ−τ′),

which induces in the tight-binding hamiltonian an effective overlap integral which depends on the first neighbours of the origin and :

 ⟨τ|Kxeh|τ′⟩=∑ρJ(τ,τ′ρ)

To lowest order, when , this adds a local term to the direct term on the first neighbours:

 Jτ=J(τ,τ,ρ=0)=∫drdr′φv(r)φc(r−τ)2|r−r′|φv(r′)φc(r′−τ).

On the other hand when we sum all contributions corresponding to all separations of the distant and pairs we obtain a sum of dipolar contributions which are known to be singular here ( limit). In 3D this produces the so-called longitudinal-transversal splitting.Denisov and Makarov (1973); Toyozawa (2003) In 2D the singularity is weaker, with terms varying as .Wu et al. (2015); Qiu et al. (2015) This will be discussed elsewhere.

## Appendix C An improved Wannier model

The tight-binding model developed in the main text is based on several approximations. To lowest order in the Wannier functions corresponding to the valence and conduction bands are taken as the orbitals on the nitrogen and boron sites, respectively. Then the kinetic energy part of the exciton hamiltonian is approximated by its second order term in and finally the Coulomb matrix elements are calculated using atomic orbitals, i.e. Wannier function of zero order. We will see that higher order terms induce corrections of order . This is not negligible but has been implicitly taken into account via our fitting procedure for the interactions where on the triangular lattice. Problems arise because, to higher order, other interactions become allowed, in particular when . Let us then define more accurate Wannier state from the exact eigenstates defined in Eq. (2). We choose the phases such that these Wannier states reduce to the atomic states when . Then, to linear order in :

 |m+⟩w=1√N∑ke−i.k.m|k+⟩≃|mB⟩−t2Δ∑τ|mB−τ⟩|n−⟩w=1√N∑ke−i.k.n|k−⟩≃|nA⟩+t2<