Two-photon absorption in two-dimensional materials: The case of hexagonal boron nitride

We calculate the two-photon absorption in bulk and single layer hexagonal boron nitride (hBN) both by an ab-initio real-time Bethe-Salpeter approach and by a the real-space solution of the excitonic problem in tight-binding formalism. The two-photon absorption obeys different selection rules from those governing linear optics and therefore provides complementary information on the electronic excitations of hBN. Combining the results from the simulations with a symmetry analysis we show that two-photon absorption is able to probe the lowest energy $1s$ states in the single layer hBN and the lowest dark degenerate dark states of bulk hBN. This deviation from the"usual"selection rules based on the continuous hydrogenic model is explained within a simple model that accounts for the crystalline symmetry. The same model can be applied to other two-dimensional materials with the same point-group symmetry, such as the transition metal chalcogenides. We also discuss the selection rules related to the inversion symmetry of the bulk layer stacking.


I. INTRODUCTION
Two-photon absorption (TPA) is a nonlinear optical process in which the absorption of two photons excites a system to a higher energy electronic state. Two-photon transitions obey selection rules distinct from those governing linear excitation processes and thereby provide complementary insights into the electronic structure of excited states, as has been demonstrated in molecular systems 1 and bulk solids. 2 In particular it is frequently argued that one-photon processes are only allowed for excitons of s symmetry whereas p states can be observed in TPA. These rules can be derived within a continuous hydrogenic model for excitons where full rotational symmetry is assumed. They have been invoked to analyze the excitonic effects observed in carbon nanotubes, 3 and also recently in bulk h-BN. 4 Actually these rules are not generally valid if the genuine crystalline symmetry is taken into account, but they can at least provide interpretations in terms of high or low oscillator strength. 5 Nonlinear optical properties of two-dimensional (2D) crystals, and as such the TPA, have been recently the object of several experimental and computational studies. For example, a giant TPA has been reported 6,7 for transition metal dichalcogenides (TMDs) which has been attributed to the peculiar optical properties of 2D crystals; in another study on TMDs, TPA has been used to probe excited states which are dark in linear optics. 8 As a consequence, there is a revival of interest about the selections rules for two-photon transitions. 9 In TMDs, the symmetry of the excitonic states is related to the existence of gaps at the K points of the Brillouin zone. 10 Excitonic effects are strong and the exciton wave functions are fairly localized, so that the low threefold symmetry plays an important role. 11 Although it has been first argued that the usual selection rules based on the hydrogenic model are also valid, 12 more accurate studies have shown that this is actually not the case, for one-photon as well as for two-photon processes. [13][14][15][16] Here we analyze the case of the h-BN single layer and bulk, which have the same lattice symmetry as the TMDs and very strongly bound excitons. We combine tight-binding calculations 14 of the two-photon transition probability with sophisticated ab initio real-time Bethe-Salpeter simulations 17 of the two-photon resonance thirdorder susceptibility. This combination is a unique feature of this work: on the one hand the tight-binding calculations allow us to identify the symmetry properties of the excitons, on the other hand the ab initio real-time Bethe-Salpeter simulations provide the TPA spectra-to our knowledge the first ab initio TPA spectra at this level of theory-which can be compared quantitatively with experiment. These two very different approaches show consistently that the TPA is able to probe the lowest 1s exciton in the bulk and single layer h-BN.
In Sec. II we discuss the choice of h-BN as a case study and describe the tight-binding modeling of its electronic and optical properties. In Sec. III we detail how the twophoton transition probabilities and the two-photon resonance third-order susceptibility are obtained respectively within the tight-binding and the ab initio framework. We then show and compare the results of the ab initio realtime simulations (Sec. IV A) and of the tight-binding calculations (Sec. IV B) for the single layer h-BN. We also contrast the case of the single layer with the bulk, highlighting the role of the inversion symmetry. Finally, we discuss the selection rules for one-and two-photon processes and on the basis of our results we clarify few recent experimental works on nonlinear optical properties of 2D crystals and bulk h-BN(Sec. IV C).

A. Choice of the system
We have chosen h-BN monolayer as a case study for several reasons. One reason is the abundance of experimental studies (luminescence, 4,18-25 X-rays, 26,27 or electron energy loss spectroscopy. [28][29][30] ) on the electronic and optical properties of both the bulk and the single layer. These studies, supported by theoretical investigations, [31][32][33][34][35][36] have shown that h-BN is an indirect band gap insulator whose optical absorption spectrum is dominated by strongly bound excitons. Then, the electronic structures of the h-BN monolayer and of the bulk structure are fairly well known-which is convenient for the tight-binding model. 31,37 A h-BN monolayer is simply the honeycomb lattice where B and N atoms alternate on the hexagons. The bands close to the gap are built from the π states. In the case of the monolayer, there are just a single valence π band and a conduction π * band which are nearly parallel along the KM direction of the Brillouin zone. The gap is direct at K point and about 7 eV. 14 It is of interest as well to compare the TPA in the single layer and in the layered bulk system. The stablest bulk structure corresponds to a so-called AA stacking where B and N atoms alternate along lines parallel to the stacking axis. The periodicity along this axis is twice the inter-planar distance and the lattice parameters used in our simulation are a = 2.5Å and c/a = 2.6. 38 In the bulk, there are then two π and two π * bands and the gap becomes indirect between a point close to K (valence band) and M (conduction band). 10 Another reason for choosing h-BN is the strong bound excitons in the absorption spectrum of both the monolayer and the bulk 14,33,34 since as mentioned before, we expect the effect of the low threefold symmetry to be more visible for very localized excitons. More practically, spectral features corresponding to strongly localized exciton converge easily with the numerical parameters within the ab initio framework allowing for accurate and not too cumbersome calculations. Furthermore, while in the absorption spectrum of the single layer there is only a pair of degenerate excitons, 14 in the absorption spectrum of bulk h-BN there are two degenerate pairs of excitons: one dark pair, which is the lowest in energy, and one bright pair (Davydov splitting). 31,32,34 Probing the lowest excitons in bulk h-BN remains a challenge. We expect that since the system has an inversion symmetry, the dark states in the absorption become bright in the TPA and vice versa, so that it may be possible to probe the lowest excitons in the bulk h-BN in TPA. In Sec.IV we show that this is indeed the case.
Finally, other systems of interest, such as the TMDs, have the same lattice symmetry, so that results depending on the symmetry properties can be generalised or can be used as a starting point when studying those systems.

B. Tight-binding model
In tight-binding formalism, we used a simple model that describes top valence and bottom conduction bands close to the gap, and which is stable to catch the most relevant physics of optical excitation in h-BN. It turns out that the valence states close to the gap are almost completely centred on the nitrogen sites whereas the conduction states are centred on the boron sites. For the monolayer this has allowed us to derive a very simple tight-binding model characterized by two parameters: an atomic parameter ∆ related to the atomic levels, +∆ for boron, and −∆ for nitrogen, and a transfer integral between first neighbours −t, t > 0. In this model the direct gap is exactly equal to 2∆. Close to the gap the electronic structure can be simplified and the electrons in the conduction band (valence band) can be considered as moving on the B (N) sublattice with an effective transfer integral equal to t 2 /2∆. The Wannier states associated with the π and π * bands are then mainly localized on N and B sublattices with small components ±(t/2∆) on the other ones. The best fit to ab initio data are obtained with t = 2.30 eV, ∆ = 3.625 eV. 14 Notice than when ∆ = 0 the model is that of graphene where the gap vanishes at K points. This discrete tightbinding model has a continuous counterpart when expanding the equations close to K in reciprocal space. The band structure is then described within a Dirac model for massless 2D electrons. Conversely, in the case of the h-BN monolayer we can use a continuous massive Dirac model. The simplest extension of the 2D TB model to bulk h-BN consists in introducing transfer integrals between nearest neighbours along the stacking axis. 10,36

C. Optical matrix elements
To first order (one-photon processes), the coupling with the electromagnetic field is described using the hamiltonian H I1 = −e p.A/m, where A is the vector potential varying as e −iωt . To this order, where v is the velocity operator and H is the hamiltonian of the unperturbed system. Using the tightbinding scheme, where |n denotes a π state at site n, we have v nm = n|v|m = (n − m)t nm /i , where t nm is the transfer integral between sites n and m. Using our simple model for the h-BN sheet, we keep only first neighbour integrals −t and the matrix element couples valence states |m v to conduction states |n c . Actually, as mentioned above, the atomic states should be replaced by Wannier states: where |m N and |n B denote the genuine atomic states centred on the N and B atoms, and where the three τ vectors are the first neighbour vectors n B − m N . The corresponding Bloch functions are given by: where the |k N (B) are the Bloch functions built from the atomic orbitals, |k N (B) = 1 √ N n∈N (B) e ik.n |n , and where f (k) = τ e ik.τ , and N is the number of lattice cells. To lowest order in t/∆, we have therefore: At low energy we can expand k close to the K and K (equivalent to −K) points. Up to a phase factor which depends on the choice of the K points, we obtain: where v F = 3 2 a bn t is the slope of |(f (k + q)| for small q, |(f (k + q)| v F q when q → 0, and a bn is the distance between first neighbor B and N atoms. e is the vector characterizing the light polarization, where ±K c are the wave functions at the ±K points. The matrix elements are therefore maximum for circular polarizations of opposite signs in K and K valleys. This result is now fairly well known in the case of transition metal dichalcogenides but it is worth insisting on a particular point. The optical coupling between valence and conduction bands is not related here to the symmetry of the π states of nitrogen and boron. Indeed, when considering polarization vectors within the sheet plane, they behave as s states. The coupling is due to the fact that the Wannier functions of the π and π * bands are centred on different sublattices. Since furthermore there is no symmetry centre connecting these sublattices a finite dipolar moment is allowed, of the circular type here, because of the threefold symmetry. 39 D. One-photon absorption for independent single particles Within a one-particle model and neglecting the photon wave vector, the one-photon absorption is proportional to the transition probability W given by the Fermi golden rule: where E is the electrical field amplitude. In our model the valence and conduction bands are symmetric, is the band energy of a triangular lattice with transfer integrals equal to t 2 /∆ and centred on 2∆ + 3t 2 /∆. 14 Neglecting the dependence on k of the matrix element, and replacing E = ω by 2∆ close to the gap, we recover the usual formula for the absorption, proportional to the joint density of states, which is here also proportional to the densities of states of the valence and conduction bands. 14

E. Excitons
In the presence of electron-hole interactions, excitonic effects come into play. They are usually treated using a Bethe Salpeter formalism, which, within our TB formalism, can be reduced to an effective Wannier-Schrödinger equation for electron-hole pairs.

Formalism
Let then |Φ be such an excitonic state. In our model for the monolayer, it can be written as: where a + kc and a kv are the usual creation and annihilation operators for conduction and valence electrons and |∅ is the unperturbed ground state. It is convenient to work in real space, using the Wannier basis for electrons and holes, and to write: Within this formalism the action of the velocity operator on the ground state can be written: (8) Finally we define electron-hole states in reciprocal and real spaces: One should be careful in distinguishing the variables associated with the motion of the electron-hole pair as a whole from those describing the relative motion of the electron and of the hole. The latter ones are k and R.
The position of the pair is defined here by the position of the hole, whereas its motion as a whole is characterized by Q = k c −k v which is a good quantum (conserved) number, taken here equal to 0. Actually |Rvc is nothing but the Bloch function for the motion of the corresponding pair for Q = 0. The extension of the formalism to Q = 0 will be treated elsewhere. In the following we drop the subscripts v and c. Eq.(8) can then be rewritten as:

Bethe-Salpeter Wannier equation
We recall here the formalism derived in Ref. [14]. We start from the usual simplified Bethe Salpeter equation which is just an effective Schrödinger equation for electron-hole pairs: In this equation, E kc − E kv plays the part of a band (or kinetic) energy, approximated here by 2∆ + t 2 |f (k)| 2 /∆. The interaction term K eh contains a direct Coulomb contribution and an exchange term which will be neglected here. Coming back in real space, we obtain an excitonic hamiltonian H eh : According to our definitions the vectors R have their origin on a hole (nitrogen) site and their extremities on the electron (boron) triangular sublattice. We then see that the above hamiltonian is similar to a tight-binding hamiltonian for electrons moving on a triangular lattice in the presence of an "impurity" at the origin (at the centre of a triangle). If the potential U R is known, this is a classical problem which can be solved with standard techniques, using diagonalization or Green function techniques.

Optical matrix element
In the many body excitonic language, the relevant matrix element is: 40 If we define a dipole matrix element d Φ = τ τ τ |Φ , the transition probability associated with exciton Φ is given by: where E Φ is the exciton energy. Since E Φ ∼ ∆, this probability is of order (t/∆) 2 . In the case of the BN monolayer, the point symmetry is that of the C 3v group and vectors transform as the two-dimensional E representation of this group. The exciton is therefore "bright", d Φ = 0, if |Φ also transforms as E. Notice also that only the local components of the wave function are involved in d Φ . This is the discrete equivalent of the classical Elliott theory: Within a continuous model the matrix element is proportional to Φ(r = 0). 41,42 In our case, the ground state exciton of the monolayer, discussed in detail in Ref. [14], does transform as E. If the position of the hole is fixed at the origin, the exciton wave function extends principally on the first neighbours and its oscillator strength is very large. More precisely, let Φ 1 , Φ 2 , Φ 3 be the three components τ |Φ . It is always possible to choose a basis (Φ + , Φ − ) of the E representation transforming in a "chiral" manner when rotated by 2π/3, Φ ± α ∝ e ±i(2απ/3) . In the language used for Wannier-Mott excitons, these two circularly polarized components can be associated with the one-dimensional representations of the C 3 group at K and K' points. Hence the notation of s states in each valley used in the case of TMD. Here the intervalley coupling is more important and the classification of states based on group theory is more accurate.
A more compact formulation for the probability transition W is the following: where the sum over i is over all exciton states. Finally, using the Green function G(z) = (z − H eh ) −1 , it is sufficient to calculate ∅|v G(z) v|∅ . 43

III. TWO-PHOTON ABSORPTION
The TPA is related to the nonlinearity in the attenuation of a laser beam. Considering the intensity I of a beam propagating alongẑ, at the lowest order beyond the linear regime, the attenuation is a nonlinear function of I: where α is the linear absorption coefficient, and β is the two-photon absorption coefficient (see e.g. Ref. [44]) In what follows, we detail how we evaluate the TPA from either the third-order susceptibility extracted from ab initio real-time dynamics (Sec. III A) and from the transition probability of a two-photon process within a second-order perturbation treatement (Sec. III B 1).

A. Nonlinear susceptibilities from ab initio real-time dynamics
The TPA coefficient in Eq. (20) is related to the imaginary part of the third-order response function χ (3)44, 45 : where n 0 (ω) is the refractive index. In two-photon measurements, the incoming laser frequency ω is set around half of the excited state energy ω 0 we want to probe 2ω ω 0 .
In this frequency region, well below the band gap, n 0 (ω) is positive, monotonic and slowly varying, therefore the peaks of β(ω) originates only from the imaginary part of χ (3) , that is the quantity we will extract from the real-time simulations. 17 In the real-time simulations, the electronic system is excited by a time-dependent monochromatic homogeneous field. The time-evolution of the system is given by the following equation of motion for the valence band states, where |v mk is the periodic part of the Bloch states. In the r.h.s. of Eq. (22), H MB k is the effective Hamiltonian derived from many-body theory that includes both the electron-hole interaction and local field effects. The specific form of H MB k is presented below. The second term in Eq. (22), E ·∂ k , describes the coupling with the external field E in the dipole approximation. As we imposed Born-von Kármán periodic boundary conditions, the coupling takes the form of a k-derivative operator ∂ k . The tilde indicates that the operator is "gauge covariant" and guarantees that the solutions of Eq. (22) are invariant under unitary rotations among occupied states at k (see Ref. [46] for more details).
We notice that we adopt here the length gauge, which presents several advantages for ab initio simulations. 17 Comparing with the velocity gauge approach-that is used in the tight-binding model-the second term of Eq. (22) includes both the A and A 2 contributions present in the velocity gauge. The two gauges are equivalent as shown in Appendix .
From the evolution of |v mk in Eq. (22) we calculate the time-dependent polarization P along the lattice vector a as where S(k, k + q) is the overlap matrix between the valence states |v nk and |v mk+q , Ω c is the unit cell volume, f is the spin degeneracy, N k is the number of k points along the polarization direction, and q = 2π/(N k a).
The polarization can be expanded in a power series of the incoming field E j as: (24) where the coefficients χ (i) are a function of the frequencies of the perturbing fields and of the outgoing polarization. As explained above, the two-photon absorption is proportional to the imaginary part of the two-photon resonance third-order sosceptibility χ In order to extract the χ (3) coefficients we resort to a technique similar to Richardson extrapolation. 47 In practice, we perform three different simulations with the incoming electric field at frequency ω and intensities corresponding to the amplitudes E, E/2 and E/4. The polarization resulting from each simulation can be expanded in the field as: Then we combine the three polarizations so to cancel out the linear and quadratic contributions and we obtain: . .
The calculation is repeated for all ω in the desired range of frequencies.
The level of approximation of the so-calculated susceptibilities depends of the effective Hamiltonian that appears in the r.h.s. of Eq. (22). Here we work in the so-called time-dependent Bethe-Salpeter framework that was introduced in Ref. [48]. In this framework the Hamiltonian H M B k reads: where H KS k is the Hamiltonian of the unperturbed (zerofield) Kohn-Sham system 49 , ∆H k is the scissor operator that has been applied to the Kohn-Sham eigenvalues, the term V h (r)[∆ρ] is the time-dependent Hartree potential 17 and is responsible for the local-field effects 50 originating from system inhomogeneities. The term Σ SEX is the screened-exchange self-energy that accounts for the electron-hole interaction, 51 and is given by the convolution between the screened interaction W and ∆γ. In the same equation: ∆ρ ≡ ρ(r; t) − ρ(r; t = 0) is the variation of the electronic density and: ∆γ ≡ γ(r, r ; t) − γ(r, r ; t = 0) is the variation of the density matrix induced by the external field E.
In the limit of small perturbation Eq. (29) and Eq. (22) reproduce the optical absorption calculated with the standard GW + BSE approach, 51 as shown both analytically and numerically in Ref. [48].

B. Tight binding model
In tight-binding the second order (two-photon processes) appears in the development of (p − eA) 2 /2m. In second order perturbation theory with respect to A, there are two terms. The first one is related to the linear term (p−eA) treated in a second order perturbation theory, and the second one comes from the A 2 term treated to first order. The latter one is frequently considered as non relevant in the optical regime, which is not the case here as explained below. Let us begin with the first contribution.

Second order perturbation theory
To second order, the transition probability P i→f from an initial state |i towards a final state |f is given by: 52 where W ij is the one-photon matrix element towards the intermediate state j. A convenient approximation consists here to replace E j by a mean energy between the initial energy and the exciton levels close to the bottom of the conduction band. 12,53,54 The sum in the numerator can then be freely performed, and we are back to a situation similar to the first order calculation, where now the relevant matrix element is equal to f |W 2 |i . Since we are looking at transition close to the gap, ω ∼ ∆, and the denominator is also of order ∆. The relevant matrix element is equal to Φ|(v.e)(v.e)|∅ . As seen above, the first velocity operator operating on the ground state generates electron-hole states. Therefore, the second one couples different electron-hole states. Since it is related to the kinetic energy part, it operates separately on the electron state and on the hole state, and therefore gives contributions proportional to the velocity of the one-particle states. The corresponding matrix elements in real space are equal to − i t 2 ∆ a, where a is a first neighbour distance on the triangular lattice. The final complete matrix element Φ|(v.e)(v.e)|∅ /∆ is therefore of order t 3 /∆ 2 . There is no reason that it vanishes identically for the ground state exciton. Actually the (tensorial) product of velocity operators transform as E × E which also contains E. In the continuous limit however, it can be shown 12 that Φ|v ⊗ v|∅ ∼ ∇Φ(r)| r=0 , which implies that only p states are bright. We are in a typical situation where precise selection rules based on the exact crystalline symmetry allow transitions which become forbidden if an approximate (higher) spherical symmetry is assumed. 5

A 2 term
In principle the A 2 term is local and its influence is negligible in the optical regime where the wave length is much larger than interatomic distances. At least is this true when using the full hamiltonian. In band theory we project the hamiltonian on the subspace defined by the number of bands taken into account, and the correct method to include gauge-invariant coupling with the electromagnetic field is to make the so-called Peierls substitution. This generates non-local coupling at all orders in A. More precisely in our case, we make the substitution: To first order, we can check that this generates the coupling H I1 used above, i.e. : and therefore: where a bn = |τ | and Φ ± τ is the amplitude τ |Φ ± of the circularly polarized state defined previously. To second order, we obtain: from which we deduce that: The last result is remarkable. It shows first that the matrix element is similar to that of the first order term and therefore that the corresponding two-photon process should be strongly visible. Then we see that the favored circular polarization for the two-photon process is the opposite of that of the one-photon process. Finally this direct contribution to first order in A 2 is stronger than the contribution discussed above coming from the second order perturbation theory. The latter one is of order t 2 /∆ 2 less than the former one. It is interesting also to analyze what happens in reciprocal space. Then the Peierls substitution amounts to replace k by k−eA/ and the non-vanishing A 2 contributions appear only if f (k) is expanded up to k 2 terms. In the context of studies on graphene, this corresponds to including so-called warping effects. In the context of k.p methods such effects would also appear if coupling with other conduction and valence bands are included. 55

One-and two-photon absorption in bulk h-BN
As mentioned previously and, as far as the band structure is concerned, the extension to three dimensional stackings of BN layers is fairly simple. On the other hand the excitonic formalism should also be extended to the case where there are several atoms per unit cell. In practice, if we continue to define exciton states using the separation between electrons and holes, we must add labels indicating in which type of plane they are. The general corresponding TB formalism is described elsewhere, 36 but if we are principally interested in the ground state excitons, the discussion becomes simpler. Actually the 1s exciton of the monolayer discussed above remains confined within a single plane in bulk h-BN if the hole is fixed in this plane. 31 Along the 0z stacking direction we have therefore to treat a problem similar to that of a Frenkel exciton, where the 2D exciton plays the part of an atomic excitation. We can then build two different excitonic Bloch functions along the 0z stacking direction which we call |Φ A and |Φ A . Introducing interplanar transfer integrals in the TB model couples these two states, and since there is an inversion center between any pair of A and A planes the final exciton eigenstates should be the usual bonding and antibonding states |Φ A ± |Φ A . The 1s state of the monolayer becomes two split states (Davydov splitting) separated by a small energy (see also Ref. [36 and 56]). ab initio calculations show that the bonding state is the ground state. Since the threefold symmetry is preserved, the two states are still themselves twofold degenerate.
Let us now look at the one-photon absorption process. We still consider a polarization e parallel to the planes. The total d Φ dipole for the AA stacking is therefore equal to d Φ A ± d Φ A . On the other due to the inversion symmetry d Φ A = −d Φ A and finally the dipole vanishes for the bonding state whereas the upper anti-bonding state is bright. In the case of the two-photon process, we can use a similar argument, but since nowq Φ A =q Φ A , the situation is reversed: the bright exciton is the lower bonding one.

A. ab initio calculations
All operators in Eq. (22) and Eq. (29) are expanded in the basis set of the Kohn-Sham band states which can be obtained from a standard DFT code. Specifically, we used the Quantum ESPRESSO code 57 where the wavefunctions are expanded in plane waves with a cutoff of 60 Ry and the effect of core electrons is simulated by norm-conserving pseudopotentials. A 12 × 12 × 4 (12 × 12 × 1 for the single layer) k-point shifted grid has been used to converge the electronic density. The band states are obtained from the diagonalization of the Kohn-Sham eigensystem. In order to simulate an isolated h-BN layer we used a supercell approach with a layer-layer distance of 20 a.u. in the perpendicular direction. The scissor operator entering Eq. (29) is chosen so as to reproduce the position of the first bright excitons in the absorption spectra of bulk and monolayer h-BN from Ref. [33]. We expanded |v mk in terms of Kohn-Sham eigenstates and we evolved the coefficients of the bands between 2 and 7 in the monolayer (7 and 12 in the bulk) in Eq. (22). We used a 15 × 15 × 5 (12 × 12 × 1 in the single layer) kpoints Γ-centered sampling in the real-time simulations which guarantee the convergence of the first peak in the spectra.
Equation (22) is solved numerically 58 for a time interval of 120 fs using the numerical approach described in Ref. [46](originally taken from Ref. [59]) with a time step of ∆t = 0.01 fs, which guarantees for numerically stable and sufficiently accurate simulations. A dephasing term corresponding to a finite broadening of about 0.05 eV is introduce in order to simulate the experimental broadening. 17 In Figure 1, panel (a), the two-photon resonant thirdorder susceptibility at χ (3) yyyy (−ω; ω, ω, −ω) proportional to the TPA-is compared with the imaginary part of the dielectric constant 2 (ω) for the monolayer h-BN. To facilitate the comparison there is a factor 2 between the energy scale of the spectra. The first peak of the TPA is found exactly at half the photon energy of the first peak of the imaginary part of the dielectric constant. That means that the two-and one-photon absorption are resonant with the same exciton. In panel (b) the same comparison is shown for bulk h-BN. In this case the first peak of the TPA is found 0.076 eV below half the photon energy of the first peak of 2 . That means that the two-photon absorption is resonant with an exciton at lower energy than the one of one-photon absorption. We also obtained 2 (ω) by solving the standard GW +Bethe-Salpeter equation 60 and diagonalized the excitonic two-particle Hamiltonian. We found-in agreement with previous studies 34 that the lowest exciton in the linear optical response of bulk h-BN is indeed dark. The position of this dark exciton is consistent with the splitting deduced from the TPA calculations. The ab initio results are then fully consistent with the discussion in Sec. III B 3. In the monolayer the ground state exciton is visible in both one-photon and two photon absorption. Instead, in bulk AA BN, the lowest exciton (pair) is dark in linear optics but becomes visible in two photon absorption. In fact, as explained in Sec. III B 3, the two lowest exciton pairs in the bulk are due to the bonding and anti-bonding combination of excitons in each layer and thus obey different selection rules.
Finally, broader features at higher energies are visible in the spectra. Those features are known to originate from both in-plane excitons with different symmetries 14 and from interplanar excitations. Notice that those peaks may not be fully converged. 61

Monolayer
In the case of the monolayer we have used the TB hamiltonian (12) with the parameters determined in Ref. [14] and have calculated the one-photon absorption spectrum from the Green function I|G(z)|I , |I = τ (τ .e)|τ , which is actually independent of the choice for e. This is conveniently done using the recursion method. A cluster of 4.10 4 atoms is used and 100 recursion levels are calculated. The spectrum, proportional to the imaginary part of the optical dielectric constant is also proportional to the imaginary part of the Green function. We have first calculated the spectrum in the absence of excitonic effects (U = 0 in the hamiltonian). It is actually quite similar to the density of states of the triangular lattice. 14 In the presence of excitonic effects the excitons appear as bound states below the continuum. We have checked that they appear at the same positions as those determined from a full diagonalization of the hamiltonian. 14 Furthermore by comparing the optical spectrum to that obtained from a Green function matrix element without any particular symmetry (here τ |G|τ ) so that all excitons have a finite weight, we can check the selection rules: the non degenerate exciton of symmetry A 1 is actually dark in the optical response (Fig. 2).
In the case of two-photon absorption, the calculation of the main contribution due to H I2 , the quadratic term in A, can be calculated in a similar way; it suffices to modify the matrix element of the Green function, by taking as initial vector |I = τ (τ .e) 2 |τ . The ground state exciton is then found to be bright as expected.

Bulk AA'
We have generalized the TB formalism to the 3D AA stacking. Once the appropriate hamiltonian is defined the recursion method can again be used. Here, the cluster considered contains 8.10 4 atoms and 100 recursion levels are calculated. According to our previous discussion, to study the two split ground state excitons we have now to use starting states |I which are bonding and antibonding combinations of excitons in A and A planes . The results are shown in Figs. 3. As predicted in the previous section we do find that the lower bonding state is dark for onephoton absorption and bright for two-photon absorption, and conversely for the upper antibonding state.

C. Selection rules
Let us summarize the selection rules which apply to h-BN as well as to TMD.

One-photon selection rules
We consider first the case of the monolayer, and recall the usual rules for Wannier-Mott excitons within the continuous hydrogenic model. The exciton wave function is written in the form Φ(r h , r e ) = φ k0c (r e )φ k0v (r h )g(r e − r h ), where the φ k0 are the single particle Bloch functions at point k 0 corresponding to the considered direct gap, and g(r), the envelope function, is the solution of the hydrogenic-like excitonic equation for the relative coordinate r = r e −r h . In our case k 0 is a K or K point. Then the optical matrix element is the standard one-particle matrix element at this point weighted by g(r = 0). The full excitonic symmetry is the symmetry of the above product, but notations generally use the symmetry of g(r). Optical one-photon transitions are then only allowed for s states. On the other hand the direct optical transition is allowed since the standard matrix element varies as v F (e x ±ie y ) (see Eq. (6)), corresponding to polarisations σ ± , depending on the valley K or K . Hence the rule for the exciton one-photon transitions: only s states are allowed and the optical transitions involve σ ± polarisations. These rules are weakened if the dependence on k of the matrix elements is taken into account. 16 At higher values of k, deviations from an isotropic hamiltonian appear, with so-called warping effects, and the usual orbital selection rules are modified. 62 Furthermore, when the exciton size decreases, its size in reciprocal space increases and intervalley coupling is no longer negligible. Although k.p method can be used, 15 it is then much more convenient to work in real space so as to take the triangular symmetry fully in account. 14 The excitonic states characterized by the wave function Φ R = R|Φ ± for the monolayer can then be classified according to the representations of the C 3v point group. Among the three representations A 1 , A 2 and E, only the two-dimensional representation E is optically active, as discussed in Section II E 3. Let us precise the correspondence between the continuous description in terms of s, p, . . . states characterizing the symmetry of the envelope function and the present description using the discrete crystal symmetry. For that, we have to include the symmetry of the product of φ k0c (r e )φ k0v (r h ) at K point. The relevant group of vector K is the C 3 point group so that in our case this product is multiplied by e ±2iπ/3 under a rotation of ±2π/3. In other words it transforms as e iϕ if ϕ is the azimuthal angle. For a level of symmetry characterized by an angular momentum m, the envelope function of the corresponding exciton states varies as e ±imϕ . When m = 0 the level is twofold degenerate and the full wave function varies as e i(1±m)ϕ . The same is true at K point provided ϕ is changed into −ϕ. So finally the level shows a fourfold degeneracy. This is not consistent with representation group theory which says that the degeneracy is at most equal to two in general. The crystal symmetry is lower than the continuous symmetry and the degeneracy is lifted. Consider the p states for example (m = ±1). It has been shown in Ref. [14] that the four p levels decompose according to the E + A 1 + A 2 representations. This can be analysed in a first step as an intravalley splitting into distinct 1 ± 1 states. Within the continuous model this can be attributed to Berry phase effects. 63,64 Another effect of the lower crystalline symmetry is that m tot = 1 ± m should be counted modulo three, since only rotations of order three are involved. Thus only states m tot = 0, ±1 are relevant in each valley.
Here the allowed values in valley K are 0 and 2, i.e. 0 and -1, modulo 3. In valley K the allowed values are 0 and +1. These rules are fairly well-known 11,13-16 but have been recently re-discovered and discussed in detail. 65,66 Then, if intervalley coupling is accounted for, the ±1 states combine to produce a global E symmetry, whereas the m tot = 0 states combine to form A 1 and A 2 states. Now, only E states are one-photon bright, since the velocity also transforms as E. More precisely, its compo-nents v.(x ± iŷ) transform as e ±iϕ . The ground state 1s exciton is obviously bright and form an E state with a very strong oscillator strength, but we see that the 2p states give rise also to a bright exciton. This is really a "modulo 3 " effect which transforms a forbidden m tot = 2 state into an allowed m tot = −1 state. Notice also that for a given valley the circular polarizations for the s and p states are opposite.
Consider now the bulk AA stacking. We still assume the light polarization to be within the planes. Then according to the discussion given in Section III B 3, we have just to consider a linear superposition of monolayer E states. For this stacking this gives rise to Davydov pairs of bonding and antibonding states, and only antibonding states are bright. Interlayer excitons where the hole and the electron are in different planes can be discussed as well and are relevant if the light polarization is along the stacking axis.
There are therefore two types of selection rules for such lamellar structures for a polarization in the planes. The first one ensures the existence of a dipole within the planes (exciton of E symmetry). The second depends on the constructive or destructive arrangement of the dipoles in the stacking.

Two photon selection rules
Within the continuous model, the general statement concerning selection rules is that only p states are visible, corresponding to ∇g(r) = 0. 12 What is changing here is the symmetry of the coupling with light which has now a tensorial character and varies therefore as e ±iϕ e ±iϕ for the monolayer. If the rule of m modulo 3 is added, this means that all m tot values are allowed. In the discrete case, it varies as E ⊗ E = E + A 1 + A 2 which indicates also that all excitons are in principle bright. We have seen in particular that the oscillator strength for the ground state 1s exciton is found actually to be very strong. In the case of the bulk AA stacking, the bright excitons in the Davydov pairs are the bonding states. This is a familiar selection rule: In the presence of a symmetry centre odd (even) states are one(two)-photon allowed. Then the difference with the one-photon selection rules is less important than in the usual continuous model, but at least in the case of the AA stacking combining both processes can be used to discriminate between the components of the Davydov doublets (Fig.4). The two-photon selection rules for the monolayer agree with those derived by Xiao et al. 13 who, however, do not discuss the possibility of the splitting of 2p states. Since they consider only circular polarizations, they do not discuss either the possibility of exciting m = 0 (or A 1 , A 2 ) states, but the corresponding oscillator strengths should be weak since they imply intervalley interactions.  ) is excited with a one-photon (two-photon) process, where Φ ± denote the two circularly components of the degenerate 1s exciton. Right: Bulk AA stacking; there are now two antisymmetric (AS, odd, one-photon allowed) and symmetric (S, even, two-photon allowed) degenerate states, separated by a Davydov splitting.

Experimental results
One of our main result is the prediction that the ground state 1s exciton which is bright in one-photon absorption should be bright also in two-photon processes, with opposite circular polarizations. The best experimental evidence is certainly the observation in TMD of resonant second harmonic generation (SHG). 67 In the absence of symmetry center, SHG is allowed for a dichalcogenide layer and it is actually found to be strongly resonant when the 1s level is excited in a two-photon process. Furthermore the circular polarization of the 2ω emission is actually opposite to that of the excitation. 13,68 In the case of h-BN only two-photon photoluminescence excitation (PLE) spectra are available for the bulk phase. 4 They do show a peak slightly above the one-photon main peak whereas we predict a peak below it. The situation is complicated by the fact that the gap is indirect. This is important for the interpretation of luminescence spectra but it is suspected that absorption spectra are governed by direct transitions. 35 The difference between onephoton and two-photon spectra has been interpreted as the signature of s and p states respectively. The present analysis show that this is probably not true because of the non conventional selection rules. Another argument is that the 2p and 2s are predicted to be well above the 1s states (about 1 eV for the monolayer, 0.5 eV for the bulk) and cannot therefore be involved in the observed splitting. But the precise interpretation of the observed splitting remains to be correlated with the calculated Davydov splitting.
To summarize we have performed tight-binding and ab initio calculations of two-photon absorption in monolayer and bulk boron nitride and found that at low energy it is driven by excitonic effects. The ground state 1s exciton is predominant as for one-photon absorption, indicating strong deviations from the usual selection rules which are explained within a simple TB model taking into account the crystalline symmetry.

ACKNOWLEDGMENTS
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. The research leading to these results has received funding from the European Union Seventh Framework Program under grant agreement no. 696656 GrapheneCore1. L. Sponza and J. Barjon are gratefully acknowledged for fruitful discussions.
Appendix: On the equivalence of length and velocity gauge in tight-binding Optical properties are usually treated using the socalled velocity gauge: p is replaced by p − eA. As discussed in the main text this implies that to calculate the response to second order in A we have to calculate a first order perturbation term in A 2 and a second order per-turbation term with respect to A. In the length gauge where the interaction term in the hamiltonian is equal to e r.E, there is no quadratic term. Since approximations are made, it is useful to check gauge invariance. If we refer to the discussion made in the velocity gauge in Section III B 1 we have therefore just to calculate the second order perturbation term proportional to ∆ Φ|(r.e)(r.e)|∅ . We have used the fact that E = iωA and that ω ∆. The first r generates electron-hole pairs r nm a + nc a mv |∅ . To lowest order in t/∆, r nm vanishes, and we must use the improved Wannier basis defined in Eqs (1-2), and then r nm −(t/2∆)(n − m), where (n − m) is a first neighbour vector τ , so that: The second r operator connects intraband conduction and valence states, and to lowest order in t/∆ it can be checked that r|τ τ |τ , so that, finally: ∆ Φ|(r.e)(r.e)|∅ − t 2 τ (τ .e) 2 Φ ± τ , which is exactly the result derived in Eq. (35). Thus to lowest order, the second order perturbation theory in the length gauge reproduces the first order term as a function of (A.e) 2 in the velocity gauge. The general equivalence of both gauges for non-linear responses of higher order has also been discussed recently in Ref. [69].