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

Claudio Attaccalite, Myrta Grüning, Hakim Amara, Sylvain Latil, François Ducastelle

► To cite this version:

Claudio Attaccalite, Myrta Grüning, Hakim Amara, Sylvain Latil, François Ducastelle. Two-photon absorption in two-dimensional materials: The case of hexagonal boron nitride. Physical Review B: Condensed Matter and Materials Physics (1998-2015), 2018, 98 (16), pp.165126. 10.1103/Phys-RevB.98.165126. hal-01755879

HAL Id: hal-01755879

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

Submitted on 31 Jan 2019

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Two-photon absorption in two-dimensional materials: The case of hexagonal boron nitride

Claudio Attaccalite,¹,² Myrta Grüning,³,⁴ Hakim Amara,⁵ Sylvain Latil,⁶ and François Ducastelle⁵

¹CNRS/Aix-Marseille Université, Centre Interdisciplinaire de Nanoscience de Marseille UMR 7325 Campus de Luminy, 13288 Marseille Cedex 9, France

²University of Rome Tor Vergata, Rome, Italy

³School of Mathematics and Physics, Queen's University Belfast, Belfast BT7 1NN, Northern Ireland, United Kingdom

⁴European Theoretical Spectroscopy Facilities (ETSF)

⁵LEM, UMR 104 CNRS-ONERA, Université Paris Saclay F-92322 Châtillon, France

⁶Service de Physique de l'Etat Condensé (SPEC), CEA, CNRS UMR 3680, Université Paris Saclay, Orme des Merisiers, CEA Saclay, 91191 Gif-sur-Yvette Cedex, France

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

(Received 29 March 2018; revised manuscript received 19 September 2018; published 16 October 2018)

We calculate the two-photon absorption in bulk and single-layer hexagonal boron nitride (*h*-BN) both by an *ab initio* real-time Bethe-Salpeter approach and by a 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 *h*-BN. Combining the results from the simulations with an analysis of the crystal symmetries, we show that two-photon absorption is able to probe the lowest-energy 1s state in the single-layer *h*-BN and the lowest degenerate exciton of bulk *h*-BN. This result indicates that in *h*-BN multilayer stackings with inversion symmetry one can measure the Davydov splitting by means of a combination of one and two-photon excitations. The same analysis can be applied to other two-dimensional materials with the same point-group symmetry—such as the transition metal chalcogenides.

DOI: 10.1103/PhysRevB.98.165126

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. Nonlinear optical properties of two-dimensional (2D) crystals, and as such the TPA, have been recently the object of several experiments. For example, a giant TPA has been reported [1,2] for transition metal dichalcogenides (TMDs), which has been attributed to the peculiar optical properties of 2D crystals; further, TPA has been used to image single quantum emitters embedded in a *h*-BN multilayer flake [3]; moreover, in another study on TMDs, TPA has been used to probe excited states, which are dark in linear optics [4]. In fact, 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 [5,6]. 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 [7] and also recently in bulk *h*-BN [8]. Though it is recognized that these rules are not generally valid if the genuine crystalline symmetry is taken into account, it is claimed they can at least guide interpretations in terms of high or low oscillator strength [9]. For 2D materials [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–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 (TB) calculations [14] of the two-photon transition probability with sophisticated *ab initio* real-time Bethe-Salpeter simulations [17] of the two-photon resonance third-order susceptibility. This combination is a unique feature of this work; on the one hand, the TB calculations allow us to identify the symmetry properties of the excitons, on the other hand the *ab initio* real-time Bethe-Salpeter (RT-BSE) simulations provide the TPA spectra—one of the first *ab initio* TPA spectra at this level of theory [18]—which can be compared quantitatively with experiments. From these two very different approaches (TB and RT-BSE), we consistently show 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 to obtain the TPA within both TB and RT-BSE: within the TB we expand the real-space formalism to the second order in the external perturbation and within the RT-BSE we introduce a post-processing technique based on Richardson extrapolation to extract the two-photon resonant third-order susceptibility, and thus the TPA, from the real-time polarization. We then show and compare the results of the *ab initio* real-time simulations (Sec. IV A) and of the TB 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).

## II. ELECTRONIC, OPTICAL PROPERTIES AND TIGHT BINDING MODEL OF h-BN

### 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 [8,19–25], x rays [26,27], or electron energy loss spectroscopy [28,29]) on the electronic and optical properties of both the bulk and the single layer. These studies, supported by theoretical investigations [30–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 [30,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) [39].

Another reason for choosing h-BN is the strong bound excitons in the absorption spectrum of both the monolayer and the bulk [14,31,32] 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) [30,32]. 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 generalized or can be used as a starting point when studying those systems.

### B. Tight-binding model

Many TB studies are dedicated to the electronic and optical properties of h-BN (see, e.g., Refs. [14,40]). Here, 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 centered on the nitrogen sites whereas the conduction states are centered on the boron sites. For the monolayer, this has allowed us to derive a very simple TB model characterised by two parameters: an atomic parameter Δ related to the atomic levels, +Δ for boron, and −Δ for nitrogen, and a transfer integral between first neighbors −t, t > 0. For these parameters, we used the same values of Ref. [14]. 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Δ. 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. This discrete TB 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 for graphene or continuous massive Dirac model for h-BN. The simplest extension of the 2D TB model to bulk h-BN consists in introducing transfer integrals between nearest neighbors along the stacking axis [34,39].

### C. One-photon absorption for independent single particles

To first order (one-photon processes), the coupling with the electromagnetic field is described using the Hamiltonian $H_{I1} = -e \mathbf{p} \cdot \mathbf{A}/m$, where A is the vector potential varying as $e^{-i\omega t}$. To this order, $\mathbf{p}/m = \mathbf{v} = \frac{1}{i\hbar}[\mathbf{r}, H]$, where v is the velocity operator and H is the Hamiltonian of the unperturbed system. Using the tight-binding scheme, where $|\mathbf{n}\rangle$ denotes a π state at site n, we have $\mathbf{v}_{nm} = \langle \mathbf{n}|\mathbf{v}|\mathbf{m}\rangle = (\mathbf{n} - \mathbf{m})t_{nm}/i\hbar$, 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 neighbor integrals −t and the matrix element couples valence states $|\mathbf{m}_v\rangle$ to conduction states $|\mathbf{n}_c\rangle$. Actually, as mentioned above, the atomic states should be replaced by Wannier states:

$$|\mathbf{m}_v\rangle_w \approx |\mathbf{m}_N\rangle + \frac{t}{2\Delta} \sum_{\tau} |\mathbf{m}_N + \tau\rangle, \quad (1)$$

$$|\mathbf{n}_c\rangle_w \approx |\mathbf{n}_B\rangle - \frac{t}{2\Delta} \sum_{\tau} |\mathbf{n}_B - \tau\rangle, \quad (2)$$

where $|\mathbf{m}_N\rangle$ and $|\mathbf{n}_B\rangle$ denote the genuine atomic states centered on the N and B atoms, and where the three τ vectors are the first neighbor vectors $\mathbf{n}_B - \mathbf{m}_N$. The corresponding Bloch functions are given by

$$|\mathbf{k}_v\rangle \approx |\mathbf{k}_N\rangle + \frac{t}{2\Delta} f^*(\mathbf{k})|\mathbf{k}_B\rangle, \quad (3)$$

$$|\mathbf{k}_c\rangle \approx |\mathbf{k}_B\rangle - \frac{t}{2\Delta} f(\mathbf{k})|\mathbf{k}_N\rangle, \quad (4)$$

where the $|\mathbf{k}_{N(B)}\rangle$ are the Bloch functions built from the atomic orbitals and $f(\mathbf{k}) = \sum_{\tau} e^{i\mathbf{k} \cdot \tau}$. At lowest order in
$t/\Delta$ we have that $\langle \mathbf{k}_c | \mathbf{v} | \mathbf{k}_v \rangle \simeq \langle \mathbf{k}_B | \mathbf{v} | \mathbf{k}_N \rangle$. Then following Ref. [14] in the limit $\mathbf{k} \to \mathbf{K}$, we obtain

$$ \langle \mathbf{k}_c | \mathbf{v} \cdot \mathbf{e} | \mathbf{k}_v \rangle \simeq (\pm \mathbf{K}_c | \mathbf{v} \cdot \mathbf{e} | \pm \mathbf{K}_v) \simeq v_F (e_x \pm i e_y), \quad (5) $$

where $v_F = \frac{3}{2}a_{bn}t$ and $a_{bn}$ is the distance between first neighbor B and N atoms. $\mathbf{e}$ is the light polarization vector and $\pm \mathbf{K}_c$ are the wave functions at the $\pm K$ points. Notice that the optical coupling between valence and conduction bands is not related here to the symmetry of the $\pi$ 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 $\pi$ and $\pi^*$ bands are centered on different sublattices.

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:

$$ W = \frac{2\pi}{\hbar} \frac{e^2}{\omega^2} (\varepsilon/2)^2 \sum_{\mathbf{k}} |\langle \mathbf{k}c | \mathbf{v} \cdot \mathbf{e} | \mathbf{k}v \rangle|^2 \delta(E_{\mathbf{k}c} - E_{\mathbf{k}v} - \hbar\omega), $$

where $\varepsilon$ is the electrical field amplitude. In our model, the valence and conduction bands are symmetric, $E_{\mathbf{k}_v} = -E_{\mathbf{k}_c}$, so that $E_{\mathbf{k}_c} - E_{\mathbf{k}_v} \simeq 2\Delta + \frac{t^2}{\Delta} |f(\mathbf{k})|^2$, where $|f(\mathbf{k})|^2$ is the band energy of a triangular lattice with transfer integrals equal to $t^2/\Delta$ and centered on $2\Delta + 3t^2/\Delta$ [14]. Neglecting the dependence on $\mathbf{k}$ of the matrix element and replacing $E = \hbar\omega$ by $2\Delta$ close to the gap, we recover the usual formula for the absorption, proportional to the joint density of states [14].

### D. 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. In the present formulation, we are not including the exchange term, therefore the singlet and triplet excitons are degenerate. This fact does not affect our results since triplet excitations are dark due the spin selection rules and they do not contribute to the linear and nonlinear response. Regarding the effect of the exchange term on the singlet exciton this is negligible at zero momentum (while it becomes important at finite momentum and for higher excitons) [14,41].

Let $|\Phi\rangle$ be an excitonic state. In our model for the monolayer, it can be written as

$$ |\Phi\rangle = \sum_{\mathbf{k}} \Phi_{\mathbf{k}} a_{\mathbf{k}_c}^{+} a_{\mathbf{k}_v} |\emptyset\rangle, \quad (6) $$

where $a_{\mathbf{k}_c}^+$ and $a_{\mathbf{k}_v}$ are the usual creation and annihilation operators for conduction and valence electrons and $|\emptyset\rangle$ is the unperturbed ground state.

The exciton wave function [Eq. (6)] can be found from the solution of the so-called Bethe Salpeter equation which is just an effective Schrödinger equation for electron-hole pairs:

$$ (E_{\mathbf{k}_c} - E_{\mathbf{k}_v}) \Phi_{\mathbf{k}vc} + \sum_{\mathbf{k}'v'c'} \langle \mathbf{k}vc | K_{eh} | \mathbf{k}'v'c' \rangle \Phi_{\mathbf{k}'v'c'} = E \Phi_{\mathbf{k}vc}. $$

In this equation, $E_{\mathbf{k}_c} - E_{\mathbf{k}_v}$ plays the part of a band energy. The interaction term $K_{eh}$ contains the direct Coulomb contribution.

In real space, the excitonic Hamiltonian $H_{eh}$ reads

$$ H_{eh}^{\mathbf{R},\mathbf{R}'} = \langle \mathbf{R} | H_{eh}^0 | \mathbf{R}' \rangle + U_{\mathbf{R}} \delta_{\mathbf{R},\mathbf{R}'} \quad (7) $$

where $U_{\mathbf{R}}$ is represented by a Keldysh potential, that has been parametrized as described in Ref. [14].

From the solution of [Eq. (7)] we can calculate the matrix element between the exciton state $\langle \Phi |$ and the ground state $|\emptyset\rangle$, that are the ones entering the linear response.<fnref n="1" /> Using the Wannier basis for electron and hole wave-functions it is possible to express these matrix elements as [14]

$$ W_{\Phi} \equiv \langle \Phi | \mathbf{e} \cdot \mathbf{v} | \emptyset \rangle = (it/\hbar) \sum_{\tau} \mathbf{e} \cdot \tau \langle \Phi | \tau | \rangle. \quad (8) $$

If we define a dipole matrix element $\mathbf{d}_{\Phi} = \sum_{\tau} \tau \langle \tau | \Phi \rangle$, the transition probability associated with exciton $\Phi$ is given by

$$ P_{\Phi} = \frac{2\pi e^2 (\varepsilon/2)^2 t^2}{\hbar E^2} |\mathbf{e} \cdot \mathbf{d}_{\Phi}|^2 \delta(E - E_{\Phi}), \quad (9) $$

where $E_{\Phi}$ is the exciton energy. Since $E_{\Phi} \sim \Delta$, this probability is of order $(t/\Delta)^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,” $\mathbf{d}_{\Phi} \neq 0$, if $|\Phi\rangle$ also transforms as $E$. Notice also that only the local components of the wave function are involved in $\mathbf{d}_{\Phi}$. This is the discrete equivalent of the classical Elliott theory: within a continuous model, the matrix element is proportional to $\Phi(r=0)$ [42,43]. 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 neighbors and its oscillator strength is very large. $\Phi_{\alpha}^{\pm} \propto e^{\pm i(2\alpha\pi/3)}$, these two circularly polarized components can be associated with the one-dimensional representations of the $C_3$ group at $K$ and $K'$ points. Finally, the probability transition can be calculated using the Green function $G(z) = (z - H_{eh})^{-1}$, as solution of $\langle \emptyset | \mathbf{v} G(z) \mathbf{v} | \emptyset \rangle$, for more detail see Ref. [44].

### 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 $\hat{z}$, at the lowest order beyond the linear regime, the attenuation is a nonlinear function of $I$:

$$ \frac{dI}{dz} = -\alpha I - \beta I^2, \quad (10) $$

where $\alpha$ is the linear absorption coefficient, and $\beta$ is the two-photon absorption coefficient (see, e.g., Ref. [45]).

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 treatment (Sec. III B 1).

<footnote n="1">¹The matrix element of the velocity operator does not depend on the Coulomb potential, which has been assumed to be local, and therefore commutes with the $\mathbf{r}$ operator.</footnote>
A. Nonlinear susceptibilities from ab initio real-time dynamics

The TPA coefficient in Eq. (10) is related to the imaginary part of the third-order response function $\chi^{(3)}$[45,46]:

$$ \beta(\omega) = \frac{3\omega}{2\epsilon_0 c^2 n_0^2(\omega)} \mathrm{Im}[\chi^{(3)}(-\omega; \omega, \omega, -\omega)], \quad (11) $$

where $n_0(\omega)$ is the refractive index. In two-photon measurements, the incoming laser frequency $\omega$ is set around half of the excited state energy $\omega_0$ we want to probe $2\omega \simeq \omega_0$. In this frequency region, well below the band gap, $n_0(\omega)$ is positive, monotonic and slowly varying, therefore the peaks of $\beta(\omega)$ originates only from the imaginary part of $\chi^{(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 monochromatic homogeneous field. The time evolution of the system is given by the following equation of motion for the valence-band states,

$$ i\hbar \frac{d}{dt} |v_{mk}\rangle = (H_{\mathbf{k}}^{\mathrm{MB}} + i\boldsymbol{\mathcal{E}} \cdot \tilde{\boldsymbol{\delta}}_{\mathbf{k}}) |v_{mk}\rangle, \quad (12) $$

where $|v_{mk}\rangle$ is the periodic part of the Bloch states. In the right-hand side of Eq. (12), $H_{\mathbf{k}}^{\mathrm{MB}}$ 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_{\mathbf{k}}^{\mathrm{MB}}$ is presented below. The second term in Eq. (12), $\boldsymbol{\mathcal{E}} \cdot \tilde{\boldsymbol{\delta}}_{\mathbf{k}}$, describes the coupling with the external field $\boldsymbol{\mathcal{E}}$ in the dipole approximation. As we imposed Born-von Kármán periodic boundary conditions, the coupling takes the form of a $\mathbf{k}$-derivative operator $\tilde{\boldsymbol{\delta}}_{\mathbf{k}}$. The tilde indicates that the operator is “gauge covariant” and guarantees that the solutions of Eq. (12) are invariant under unitary rotations among occupied states at $\mathbf{k}$ (see Ref. [47] 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. (12) includes both the $\mathbf{A}$ and $\mathbf{A}^2$ contributions present in the velocity gauge. The two gauges are equivalent as shown in Appendix.

From the evolution of $|v_{mk}\rangle$ in Eq. (12) we calculate the real-time polarization $\mathbf{P}_{\parallel}$ along the lattice vector $\mathbf{a}$ as

$$ \mathbf{P}_{\parallel} = -\frac{ef|\mathbf{a}|}{2\pi\Omega_c} \mathrm{Im} \ln \prod_{\mathbf{k}}^{N_{\mathbf{k}}-1} \det S(\mathbf{k}, \mathbf{k}+\mathbf{q}), \quad (13) $$

where $S(\mathbf{k}, \mathbf{k}+\mathbf{q})$ is the overlap matrix between the valence states $|v_{mk}\rangle$ and $|v_{mk+\mathbf{q}}\rangle$, $\Omega_c$ is the unit cell volume, $f$ is the spin degeneracy, $N_{\mathbf{k}}$ is the number of $\mathbf{k}$ points along the polarization direction, and $\mathbf{q} = 2\pi/(N_{\mathbf{k}}\mathbf{a})$.

The polarization can be expanded in a power series of the incoming field $\boldsymbol{\mathcal{E}}_j$ as

$$ \mathbf{P}_i = \chi_{ij}^{(1)} \boldsymbol{\mathcal{E}}_j + \chi_{ijk}^{(2)} \boldsymbol{\mathcal{E}}_j \boldsymbol{\mathcal{E}}_k + \chi_{ijkl}^{(3)} \boldsymbol{\mathcal{E}}_j \boldsymbol{\mathcal{E}}_k \boldsymbol{\mathcal{E}}_l + O(\boldsymbol{\mathcal{E}}^4), \quad (14) $$

where the coefficients $\chi^{(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 susceptibility $\chi_{ijkl}^{(3)}(-\omega; \omega, \omega, -\omega)$: i.e., the outgoing polarization has the same frequency $\omega$ of the incoming laser field,

as a result of the absorption of two and the emission of one virtual photons.

In order to extract the $\chi^{(3)}$ coefficients, we resort to a technique similar to Richardson extrapolation [48]. In practice, we perform three different simulations with the incoming electric field at frequency $\omega$ and intensities corresponding to the amplitudes $\boldsymbol{\mathcal{E}}, \boldsymbol{\mathcal{E}}/2$ and $\boldsymbol{\mathcal{E}}/4$. The polarization resulting from each simulation can be expanded in the field as

$$ \mathbf{P}(\boldsymbol{\mathcal{E}}) = \chi^{(1)}\boldsymbol{\mathcal{E}} + \chi^{(2)}\boldsymbol{\mathcal{E}}^2 + \chi^{(3)}\boldsymbol{\mathcal{E}}^3 + O(\boldsymbol{\mathcal{E}}^4), \quad (15) $$

$$ \mathbf{P}\left(\frac{\boldsymbol{\mathcal{E}}}{2}\right) = \chi^{(1)}\frac{\boldsymbol{\mathcal{E}}}{2} + \chi^{(2)}\frac{\boldsymbol{\mathcal{E}}^2}{4} + \chi^{(3)}\frac{\boldsymbol{\mathcal{E}}^3}{8} + O(\boldsymbol{\mathcal{E}}^4), \quad (16) $$

$$ \mathbf{P}\left(\frac{\boldsymbol{\mathcal{E}}}{4}\right) = \chi^{(1)}\frac{\boldsymbol{\mathcal{E}}}{4} + \chi^{(2)}\frac{\boldsymbol{\mathcal{E}}^2}{16} + \chi^{(3)}\frac{\boldsymbol{\mathcal{E}}^3}{64} + O(\boldsymbol{\mathcal{E}}^4). \quad (17) $$

Then we combine the three polarizations so to cancel out the linear and quadratic contributions and we obtain

$$ \chi^{(3)} = \frac{8 \mathbf{P}(\boldsymbol{\mathcal{E}}) - 6\mathbf{P}\left(\frac{\boldsymbol{\mathcal{E}}}{2}\right) + 8\mathbf{P}\left(\frac{\boldsymbol{\mathcal{E}}}{4}\right)}{3 \boldsymbol{\mathcal{E}}^3}. \quad (18) $$

The calculation is repeated for all $\omega$ in the desired range of frequencies.

The level of approximation of the so-calculated susceptibilities depends of the effective Hamiltonian that appears in the right-hand side of Eq. (12). Here we work in the so-called real-time Bethe-Salpeter framework that was introduced in Ref. [49]. In this framework, the Hamiltonian $H_{\mathbf{k}}^{\mathrm{MB}}$ reads

$$ H_{\mathbf{k}}^{\mathrm{MB}} = H_{\mathbf{k}}^{\mathrm{KS}} + \Delta H_{\mathbf{k}} + V_h(\mathbf{r})[\Delta\rho] + \Sigma_{\mathrm{SEX}}[\Delta\gamma], \quad (19) $$

where $H_{\mathbf{k}}^{\mathrm{KS}}$ is the Hamiltonian of the unperturbed (zero-field) Kohn-Sham system [50], $\Delta H_{\mathbf{k}}$ is the scissor operator that has been applied to the Kohn-Sham eigenvalues, the term $V_h(\mathbf{r})[\Delta\rho]$ is the real-time Hartree potential [17] and is responsible for the local-field effects [51] originating from system inhomogeneities. The term $\Sigma_{\mathrm{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 $\mathbf{W}$ and $\Delta\gamma$. In the same equation,

$$ \Delta\rho \equiv \rho(\mathbf{r}; t) - \rho(\mathbf{r}; t=0) $$

is the variation of the electronic density and

$$ \Delta\gamma \equiv \gamma(\mathbf{r}, \mathbf{r}'; t) - \gamma(\mathbf{r}, \mathbf{r}'; t=0) $$

is the variation of the density matrix induced by the external field $\boldsymbol{\mathcal{E}}$. Equations (19) and (12) are a generalization of time-dependent Hartree-Fock ones, where the exchange is screened by the static dielectric constant, calculated within random phase approximation [51] and the single-particle levels have been shifted to reproduce quasiparticle levels obtained from $G_0W_0$ calculations. In the limit of small perturbation, these equations reproduce the optical absorption calculated with the standard $G_0W_0$ + BSE approach [51], as shown both analytically and numerically in Ref. [49].

B. Tight-binding model

In tight-binding model, the second order (two-photon processes) appears in the development of $(\mathbf{p} - e\mathbf{A})^2/2m$. In the second-order perturbation theory with respect to $\mathbf{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<fnref n="2" /> term treated to first order. The latter one is frequently considered as nonrelevant in the optical regime, which is not the case here as explained below. Let us begin with the first contribution.

1. Second-order perturbation theory

To second order, the transition probability $P_{i \to f}$ from an initial state $|i\rangle$ towards a final state $|f\rangle$ is given by [52]

$$P_{i \to f} = \frac{2\pi}{\hbar} \left| \frac{1}{4} \sum_{j \neq i, f} \frac{W_{fj} W_{ji}}{E_i - E_j + \hbar\omega} \right|^2 \delta(E_f - E_i - 2\hbar\omega), \quad (20)$$

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 $\langle f|W^2|i\rangle$. Since we are looking at transition close to the gap, $\hbar\omega \sim \Delta$, and the denominator is also of order $\Delta$. The relevant matrix element is equal to $\langle \Phi|(\mathbf{v} \cdot \mathbf{e})(\mathbf{v} \cdot \mathbf{e})|\emptyset\rangle$. 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 $-\frac{i}{\hbar} \frac{t^2}{\Delta} \mathbf{a}$, where $\mathbf{a}$ is a first neighbor distance on the triangular lattice. The final complete matrix element $\langle \Phi|(\mathbf{v} \cdot \mathbf{e})(\mathbf{v} \cdot \mathbf{e})|\emptyset\rangle/\Delta$ is therefore of order $t^3/\Delta^2$. There is no reason that it vanishes identically for the ground state exciton. Actually the (tensorial) product of velocity operators transform as $E \times E$, which also contains $E$. In the continuous limit, however, it can be shown [12] that $\langle \Phi|\mathbf{v} \otimes \mathbf{v}|\emptyset\rangle \sim \nabla\Phi(\mathbf{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 [9].

2. A<fnref n="2" /> term

In principle, the A<fnref n="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 nonlocal coupling at all orders in A. More precisely in our case, we make the substitution

$$t_{nm} \rightarrow t_{nm} e^{ie(n-m)\cdot A/\hbar}. \quad (21)$$

To first order, we can check that this generates the coupling $H_{I1}$ used above, i.e.,

$$\langle \mathbf{n} | H_{I1} | \mathbf{m} \rangle = i \left( \frac{e}{\hbar} \right) (\mathbf{m} - \mathbf{n}) \cdot \mathbf{A} \ t_{nm}, \quad (22)$$

and therefore,

$$\langle \emptyset | H_{I1} | \Phi^{\pm} \rangle = i \frac{eA}{\hbar} \sum_{\tau} (\tau \cdot \mathbf{e}) \Phi_{\tau}^{\pm} \quad (23)$$

$$\sum_{\tau} (\tau \cdot \mathbf{e}) \Phi_{\tau}^{\pm} = \mathbf{e} \cdot \mathbf{d}_{\Phi^{\pm}} = -\frac{3}{2} a_{bn} (e_x \pm i e_y) |\Phi_{\tau}^{\pm}|, \quad (24)$$

where $a_{bn} = |\tau|$ and $\Phi_{\tau}^{\pm}$ is the amplitude $\langle \tau | \Phi^{\pm} \rangle$ of the circularly polarized state defined previously. To second order, we obtain

$$\langle \mathbf{n} | H_{I2} | \mathbf{m} \rangle = -\frac{1}{2} \left(\frac{e}{\hbar}\right)^2 [(\mathbf{m} - \mathbf{n}) \cdot \mathbf{A}]^2 t_{nm},$$

from which we deduce that

$$\langle \emptyset | H_{I2} | \Phi^{\pm} \rangle = -\frac{1}{2} \left(\frac{eA}{\hbar}\right)^2 t \sum_{\tau} (\tau \cdot \mathbf{e})^2 \Phi_{\tau}^{\pm} \quad (25)$$

$$\sum_{\tau} (\tau \cdot \mathbf{e})^2 \langle \tau | \Phi_{+} \rangle = \mathbf{e} \cdot \bar{\mathbf{q}} \cdot \mathbf{e} = \frac{3}{8} a_{bn}^2 (e_x - i e_y)^2 |\Phi_{\tau}^{\pm}|, \quad (26)$$

where $\bar{\mathbf{q}}_{\Phi^{\pm}} = \sum_{\tau} (\tau \otimes \tau) \Phi_{\tau}^{\pm}$ is a quadrupolar matrix element.

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<fnref n="2" /> is stronger than the contribution discussed above coming from the second order perturbation theory. The latter one is of order $t^2/\Delta^2$ less than the former one. It is interesting also to analyze what happens in reciprocal space. Then the Peierls substitution amounts to replace $\mathbf{k}$ by $\mathbf{k} - e\mathbf{A}/\hbar$ and the nonvanishing A<fnref n="2" /> contributions appear only if $f(\mathbf{k})$ is expanded up to $\mathbf{k}^2$ terms. In the context of studies on graphene, this corresponds to including so-called warping effects. In the context of $\mathbf{k} \cdot \mathbf{p}$ methods, such effects would also appear if coupling with other conduction and valence bands are included.<fnref n="2" />

3. 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

<footnote n="2">²M. Glazov (private communication); this is also discussed in Ref. [15] for the case of TMD.</footnote>
described elsewhere [34], but if we are interested in the ground state excitons, the discussion becomes simpler. Actually, the 1s exciton of the monolayer remains confined within a single plane in bulk h-BN if the hole is fixed in this plane [30]. 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 $|\Phi_A\rangle$ and $|\Phi_{A'}\rangle$. 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 $|\Phi_A\rangle \pm |\Phi_{A'}\rangle$. The 1s state of the monolayer becomes two split states (Davydov splitting) separated by a small energy (see also Refs. [34,55]). *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_\Phi$ dipole for the AA' stacking is therefore equal to $d_{\Phi_A} \pm d_{\Phi_{A'}}$. On the other, due to the inversion symmetry $d_{\Phi_A} = -d_{\Phi_{A'}}$ and finally the dipole vanishes for the bonding state whereas the upper antibonding state is bright. In the case of the two-photon process, we can use a similar argument, but since now $\bar{q}_{\Phi_{A'}} = \bar{q}_{\Phi_A}$, the situation is reversed: the bright exciton is the lower bonding one.

# IV. RESULTS AND DISCUSSION

## A. Ab initio calculations

All operators in Eqs. (12) and (19) 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 [56] where the wave functions 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. (19) 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. [31].

We expanded $|v_{mk}\rangle$ 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. (12). The dielectric constant that enters in the calculation of $\Sigma_{SEX}$ was calculated using 40 bands and 3 Ha cutoff both for the bulk and the monolayer. We used a 15 × 15 × 5 (12 × 12 × 1 in the single layer) k-point Γ-centred sampling in the real-time simulations, which guarantee the convergence of the first peak in the spectra.<fnref n="3" /> Equation (12) is solved numerically [57] for a time interval of 120 fs using the numerical approach

<figure><img src="image_3.png" /><figcaption>FIG. 1. Two-photon absorption and imaginary part of the dielectric constant in single layer (a) and bulk h-BN (b), obtained from real-time *ab initio* dynamics. The two curves have been rescaled in such a way to have the same intensity at the maximum position. Vertical lines in (b) indicate the position of the maximum.</figcaption></figure>

described in Ref. [47] (originally taken from Ref. [58]) with a time step of $\Delta 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 introduced in order to simulate the experimental broadening [17].

In Fig. 1(a), the two-photon resonant third-order susceptibility at $\chi^{(3)}_{yyyy}(-\omega; \omega, \omega, -\omega)$ proportional to the TPA is compared with the imaginary part of the dielectric constant $\epsilon_2(\omega)$ 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 $\epsilon_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 $\epsilon_2(\omega)$ by solving the standard GW+Bethe-Salpeter equation [59] and diagonalized the excitonic two-particle Hamiltonian. We found, in agreement with previous studies [32], that the lowest exciton in the linear optical response of bulk h-BN is indeed dark. The

<footnote>could require a finer k-point grid and larger distance between the layers. These effects are compensated by the choice of the scissor operator.</footnote>

<footnote n="3">³Notice that the convergence parameters used in this work are sufficient to converge the spectra shape, while its absolute position</footnote>
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- 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 antibonding combination of excitons in each layer and thus obey different selection rules.

Experimentally, the TPA has been investigated for bulk h-BN in Ref. [8] by two-photon photoluminescence excitation (PLE) spectroscopy. In agreement with our results, the lowest peaks in one- and two-photon PLE spectra do not coincide, indicating they correspond to different excitons. However, the peak in the one-photon PLE is below the peak in the two-photon PLE spectrum, apparently contradicting our findings. On the other hand, the experimental results indicate as well that phonons are playing a crucial role and the peak in the two-photon PLE spectrum is interpreted as a *phonon-assisted* two-photon absorption, which could explain the difference with our simulations that do not include phonon scattering. Considering the possible phonon modes that can contribute to phonon assisted two-photon excitations, we estimate a shift of 150 meV that could indeed account for the observed difference. In Sec. IVC, we further discuss these experimental results and their interpretation in connection with the 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 with the parameters of the simulation.<fnref n="4" />

## B. Tight-binding calculations

### 1. Monolayer

In the case of the monolayer, we have used the TB Hamiltonian (7) with the parameters determined in Ref. [14] and have calculated the one-photon absorption spectrum from the Green function $\langle I|G(z)|I\rangle$, $|I\rangle = \sum_{\tau}(\tau \cdot e)|\tau\rangle$, which is actually independent of the choice for e. This is conveniently done using the recursion method. A cluster of $4 \times 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. 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 $\langle \tau | G | \tau \rangle$) so that all excitons have a finite weight,

<figure><img src="image_4.png" /><figcaption>FIG. 2. (a) Excitonic joint density of states of the 2D BN monolayer; (b) zoom in the exciton region; top: excitonic joint density; middle: one-photon optical spectrum; bottom: two-photon optical spectrum. For convenience, a broadening of $10^{-2}$ eV has been applied, via the imaginary part of $z$ in the Green functions. The $A_1$ exciton is only dark in the one-photon spectrum. The origin of the energy scale is taken at the position of the lowest exciton, which can also be labeled as a $1s$ exciton whereas the other ones derive from $2s$ and $2p$ states. Full *ab initio* calculations including exchange effects predict the $A_1$ level to be at a higher energy [14].</figcaption></figure>

we can check the selection rules: the nondegenerate 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_{12}$, 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'\rangle = \sum_{\tau}(\tau \cdot e)^2|\tau\rangle$. The ground state exciton is then found to be bright as expected. This is in agreement with the *ab initio* spectrum in Fig. 1(a) where the lowest peak in the TPA and optical absorption coincides.

### 2. Bulk AA'

We have generalised 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 \times 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\rangle$, which are bonding and antibonding combinations of excitons in $A$ and $A'$ planes. The results are shown in Fig. 3.

<footnote n="4">⁴Excitons at higher energy are more difficult to characterize because they require a denser k-point sampling. As specified above, the k-point samplings are chosen to guarantee the convergence of the first peak in the spectra that are the focus of this study.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 3. One- (green) and two-photon (blue) absorption spectrum of AA' h-BN. The results are in agreement with the *ab initio* results in Fig. 1(b).</figcaption></figure>

We find that the lower bonding state is dark for one-photon ab-
sorption and bright for two-photon absorption, and conversely
for the upper antibonding state. Both the excitons position
(with respect to the lowest exciton) and Davydov splitting
agree with the *ab initio* results [Fig. 1(b)]. In fact, in both
calculations, the lowest exciton of the TPA spectrum is about
0.1 eV below the lowest exciton of the absorption spectrum.
Further, in the TB calculations a second peak is visible in
the TPA at about 0.3 eV above the lowest peak. Similarly,
in the *ab initio* TPA spectrum, a strong feature is visible at
about 0.5 eV above the lowest peak, which is absent in the
absorption spectrum. This 0.2 eV difference can be attributed
to neglecting of the exchange term in the TB model. In fact,
an underestimation of the energy of higher excitations has also
been also observed for this model in Ref. [14]. As well, as we
remarked before, the position of higher peaks in the *ab initio*
spectra may not be fully converged with the parameters of
the simulation, though improving the convergence parameters
would likely increase the difference between the TB and *ab
initio* results rather than reduce it.<fnref n="5" />

To summarize, the two approaches show consistently that
the ground-state exciton, which is usually labeled as 1s exci-
ton, is bright in both the monolayer and bulk TPA spectra.
For the bulk, the comparison between the one-photon and
two-photon spectra provides a measure of the Davydov split-
ting between the two 1s excitons, which we estimate to be
~0.1 eV.

C. Selection rules

Let us summarize the selection rules which apply to h-BN
as well as to TMD. We will briefly discuss the hydrogenic
model, that although inadequate [10], is the first approxima-
tion used to interpret optical spectra, then we consider the full
symmetry of h-BN.

1. One-photon selection rules

We consider first the case of the monolayer and Wannier-
Mott excitons within the continuous hydrogenic model. The

exciton wave function is written in the form Φ(,**r**<sub>h</sub>, **r**<sub>e</sub>) =
φ<sub>k<sub>0c</sub></sub>(**r**<sub>e</sub>)φ<sub>k<sub>0v</sub></sub>(**r**<sub>h</sub>)g(**r**<sub>e</sub> − **r**<sub>h</sub>), where the φ<sub>k<sub>0</sub></sub> are the single par-
ticle Bloch functions at point **k**<sub>0</sub> corresponding to the con-
sidered direct gap, and g(**r**), the envelope function, is the
solution of the hydrogeniclike excitonic equation. 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<sub>F</sub>(e<sub>x</sub> ± ie<sub>y</sub>) [see Eq. (5)],
corresponding to polarizations σ<sub>±</sub>, depending on the valley **K**
or **K'**. These rules are weakened if the dependence on **k** of the
matrix elements is taken into account [16].

The excitonic states characterized by the wave function
Φ<sub>R</sub> = 〈**R**|Φ<sup>±</sup>〉 for the monolayer can then be classified
according to the representations of the C<sub>3v</sub> point group.
Among the three representations A<sub>1</sub>, A<sub>2</sub>, and E, only the
two-dimensional representation E is optically active, as
discussed in Sec. IID. Let us precise the correspondence
between the continuous description in terms of *s*, *p*, . . . states
characterising 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 φ<sub>k<sub>0c</sub></sub>(**r**<sub>e</sub>)φ<sub>k<sub>0v</sub></sub>(**r**<sub>h</sub>) at **K** point. The relevant group of vector
**K** is the C<sub>3</sub> point group so that in our case this product is
multiplied by e<sup>±2iπ/3</sup> under a rotation of ±2π/3. For a level
of symmetry characterized by an angular momentum *m*, the
envelope function of the corresponding exciton states varies
as e<sup>±imφ</sup>.

When *m* ≠ 0, the level is twofold degenerate and the full
wave function varies as e<sup>i(1±m)φ</sup>. The same is true at *K*'
point provided *φ* is changed into −*φ*. So finally the level
shows a fourfold degeneracy. Since the crystal symmetry is
lower than the continuous symmetry 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*<sub>1</sub> + *A*<sub>2</sub> representations. These rules are fairly
well-known [11,13–16] but have been recently rediscovered
and discussed in detail [60,61].

Then, if intervalley coupling is accounted for, the ±1 states
combine to produce a global *E* symmetry, whereas the *m*<sub>tot</sub> =
0 states combine to form *A*<sub>1</sub> and *A*<sub>2</sub> states. Now, only *E* states
are one-photon bright, since the velocity also transforms as *E*.
The ground-state 1s exciton is obviously bright and forms an
*E* state with a very strong oscillator strength, but we see that
the 2*p* states give rise also to a bright exciton.

Consider now the bulk AA' stacking. According to the dis-
cussion given in Sec. IIIB 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. 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 while the second depends on the
constructive or destructive arrangement of the dipoles in the
stacking.

2. Two-photon selection rules

Within the continuous model, the general statement con-
cerning selection rules is that only *p* states are visible,
corresponding to ∇*g*(**r**) ≠ 0 [12]. What is changing here is

<footnote n="5">⁵In fact, as higher energy excitons are more delocalized, their accurate description may require a denser *k*-point sampling. If the *k*-point sampling is too small, as it may be the case here, the excitons are artificially localized and their binding energy overestimated.</footnote>
<figure><img src="image_6.png" /><figcaption>FIG. 4. Allowed transitions towards the lowest excitonic states with circularly polarized light. (Left) Monolayer; the $Φ^+(Φ^-)$ 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.</figcaption></figure>

the symmetry of the coupling with light, which has now a tensorial character and varies therefore as $e^{\pm i\varphi}e^{\pm i\varphi}$ 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 \otimes 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 center, 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 2*p* states. Since they consider only circular polarizations, they do not discuss either the possibility of exciting *m* = 0 (or *A*₁, *A*₂) states, but the corresponding oscillator strengths should be weak since they imply intervalley interactions.

### 3. 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) [62–64]. 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,65]. In the case of *h*-BN, only two-photon PLE spectra are available for the bulk phase [8]. As already commented in

Sec. IV A, 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 [33]. The difference between one-photon and two-photon spectra has been interpreted as the signature of *s* and *p* states, respectively. The present analysis shows that this is probably not true because of the nonconventional selection rules. Another argument is that the 2*p* and 2*s* 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.<fnref n="6" /> But the precise interpretation of the observed splitting remains to be correlated with the calculated Davydov splitting.

To summarize, we have performed TB and ab initio calculations of two-photon absorption in monolayer and bulk boron nitride and found that at low energy, the spectra are dominated by excitonic effects. The ground-state 1s exciton is predominant—as for one-photon absorption—indicating strong deviations from the selection rules based on the hydrogenic model, which are frequently employed to interpret the experimental spectra. For both the bulk and 2D case, we have explained the selection rule within a simple TB model that takes into account the crystalline symmetry. Finally, the result obtained for the bulk indicates that one can measure the Davydov splitting in multilayer stackings with inversion symmetry of *h*-BN and 2D crystals with the same point-group 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 H2020 Framework Program under Grant Agreements No. 696656 Graphene Core1 and No. 785219 Graphene Core2. C.A. acknowledges PRACE for computational resources on Marconi at CINECA (Grant No. Pra16_4181L). L. Sponza and J. Barjon are gratefully acknowledged for fruitful discussions.

### APPENDIX: ON THE EQUIVALENCE OF LENGTH AND VELOCITY GAUGE IN THE TIGHT-BINDING MODEL

Optical properties are usually treated using the so-called velocity gauge: **p** is replaced by **p** − *e***A**. 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 perturbation term with respect to **A**. In the length gauge where the interaction term in the Hamiltonian is equal to *e* **r** ⋅ **ε**, there is no quadratic term. Since approximations are made, it is useful to check

<footnote n="6">⁶The 2*p* state seems inaccessible to experiments because it is at too high energy and we are not aware of experiments that measure this state in *h*-BN.</footnote>
gauge invariance. If we refer to the discussion made in the velocity gauge in Sec. III B 1, we have therefore just to calculate the second order perturbation term proportional to $\Delta\langle\Phi|(\mathbf{r}\cdot\mathbf{e})(\mathbf{r}\cdot\mathbf{e})|\emptyset\rangle$. We have used the fact that $E = i\omega\mathbf{A}$ and that $\omega \simeq \Delta$. The first $\mathbf{r}$ generates electron-hole pairs $\mathbf{r}_{nm}a_n^+ a_m^-|\emptyset\rangle$. To lowest order in $t/\Delta$, $\mathbf{r}_{nm}$ vanishes, and we must use the improved Wannier basis defined in Eqs. (1) and (2), and then $\mathbf{r}_{nm} \simeq -(t/2\Delta)(\mathbf{n}-\mathbf{m})$, where $(\mathbf{n}-\mathbf{m})$ is a first neighbor vector $\tau$, so that

$$ \mathbf{r}|\emptyset\rangle \simeq -\frac{t}{2\Delta} \sum_{\tau} \tau|\tau\rangle. $$

The second $\mathbf{r}$ operator connects intraband conduction and valence states, and to lowest order in $t/\Delta$ it can be checked that $\mathbf{r}|\tau\rangle \simeq \tau|\tau\rangle$, so that, finally,

$$ \Delta\langle\Phi|(\mathbf{r}\cdot\mathbf{e})(\mathbf{r}\cdot\mathbf{e})|\emptyset\rangle \simeq -\frac{t}{2} \sum_{\tau} (\tau\cdot\mathbf{e})^2 \Phi_{\tau}^{\pm}, $$

which is exactly the result derived in Eq. (25). Thus to the lowest order, the second-order perturbation theory in the length gauge reproduces the first-order term as a function of $(\mathbf{A}\cdot\mathbf{e})^2$ in the velocity gauge. The general equivalence of both gauges for nonlinear responses of higher order has also been discussed recently in Ref. [66].

[1] Y. Li, N. Dong, S. Zhang, X. Zhang, Y. Feng, K. Wang, L. Zhang, and J. Wang, Giant two-photon absorption in monolayer mos2, Laser Photonics Rev. 9, 427 (2015).

[2] S. Zhang, N. Dong, N. McEvoy, M. O'Brien, S. Winters, N. C. Berner, C. Yim, Y. Li, X. Zhang, Z. Chen et al., Direct observation of degenerate two-photon absorption and its saturation in ws2 and mos2 monolayer and few-layer films, ACS nano 9, 7142 (2015).

[3] A. W. Schell, T. T. Tran, H. Takashima, S. Takeuchi, and I. Aharonovich, Non-linear excitation of quantum emitters in hexagonal boron nitride multiplayers, APL Photonics 1, 091302 (2016).

[4] K. He, N. Kumar, L. Zhao, Z. Wang, K. F. Mak, H. Zhao, and J. Shan, Tightly Bound Excitons in Monolayer WSe$_2$, Phys. Rev. Lett. 113, 026803 (2014).

[5] P. R. Callis, Two-photon-induced fluorescence, Annu. Rev. Phys. Chem. 48, 271 (1997).

[6] Y.-R. Shen, The Principles of Nonlinear Optics (Wiley Inter- science, New York, 1984), p. 575.

[7] F. Wang, G. Dukovic, L. E. Brus, and T. F. Heinz, The optical resonances in carbon nanotubes arise from excitons, Science 308, 838 (2005).

[8] G. Cassabois, P. Valvin, and B. Gil, Hexagonal boron nitride is an indirect bandgap semiconductor, Nat. Photonics 10, 262 (2016).

[9] E. B. Barros, R. B. Capaz, A. Jorio, G. G. Samsonidze, A. G. Souza Filho, S. Ismail-Beigi, C. D. Spataru, S. G. Louie, G. Dresselhaus, and M. S. Dresselhaus, Selection rules for one-and two-photon absorption by excitons in carbon nanotubes, Phys. Rev. B 73, 241406 (2006).

[10] G. Wang, A. Chernikov, M. M. Glazov, T. F. Heinz, X. Marie, T. Amand, and B. Urbaszek, Colloquium: Excitons in atomically thin transition metal dichalcogenides, Rev. Mod. Phys. 90, 021001 (2018).

[11] F. Wu, F. Qu, and A. H. MacDonald, Exciton band structure of monolayer mos$_2$, Phys. Rev. B 91, 075310 (2015).

[12] T. C. Berkelbach, M. S. Hybertsen, and D. R. Reichman, Bright and dark singlet excitons via linear and two-photon spectroscopy in monolayer transition-metal dichalcogenides, Phys. Rev. B 92, 085413 (2015).

[13] J. Xiao, Z. Ye, Y. Wang, H. Zhu, Y. Wang, and X. Zhang, Nonlinear optical selection rule based on valley-exciton locking in monolayer ws2, Light: Sci. Appl. 4, e366 (2015).

[14] T. Galvani, F. Paleari, H. P. C. Miranda, A. Molina-Sánchez, L. Wirtz, S. Latil, H. Amara, and F. Ducastelle, Excitons in boron nitride single layer, Phys. Rev. B 94, 125303 (2016).

[15] M. M. Glazov, L. E. Golub, G. Wang, X. Marie, T. Amand, and B. Urbaszek, Intrinsic exciton-state mixing and nonlinear optical properties in transition metal dichalcogenide monolayers, Phys. Rev. B 95, 035311 (2017).

[16] P. Gong, H. Yu, Y. Wang, and W. Yao, Optical selection rules for excitonic rydberg series in the massive dirac cones of hexagonal two-dimensional materials, Phys. Rev. B 95, 125420 (2017).

[17] C. Attaccalite and M. Grüning, Nonlinear optics from an ab initio approach by means of the dynamical berry phase: Application to second-and third-harmonic generation in semiconductors, Phys. Rev. B 88, 235113 (2013).

[18] J. Maultzsch, R. Pomraenke, S. Reich, E. Chang, D. Prezzi, A. Ruini, E. Molinari, M. S. Strano, C. Thomsen, and C. Lienau, Exciton binding energies in carbon nanotubes from two-photon photoluminescence, Phys. Rev. B 72, 241402 (2005).

[19] K. Watanabe and T. Taniguchi, Jahn-Teller effect on exciton states in hexagonal boron nitride single crystal, Phys. Rev. B 79, 193104 (2009).

[20] P. Jaffrennou, J. Barjon, J.-S. Lauret, A. Loiseau, F. Ducastelle, and B. Attal-Tretout, Origin of the excitonic recombinations in hexagonal boron nitride by spatially resolved cathodoluminescence spectroscopy, J. Appl. Phys. 102, 116102 (2007).

[21] L. Museur, G. Brasse, A. Pierret, S. Maine, B. Attal-Tretout, F. Ducastelle, A. Loiseau, J. Barjon, K. Watanabe, T. Taniguchi, and A. Kanaev, Exciton optical transitions in a hexagonal boron nitride single crystal, Phys. Status Solidi RRL 5, 214 (2011).

[22] A. Pierret, J. Loayza, B. Berini, A. Betz, B. Plaçais, F. Ducastelle, J. Barjon, and A. Loiseau, Excitonic recombinations in hBN: From bulk to exfoliated layers, Phys. Rev. B 89, 035414 (2014).

[23] L. Schue, B. Berini, A. C. Betz, B. Plaçais, F. Ducastelle, J. Barjon, and A. Loiseau, Dimensionality effects on the luminescence properties of hbn, Nanoscale 8, 6986 (2016).

[24] J. Li, X. K. Cao, T. B. Hoffman, J. H. Edgar, J. Y. Lin, and H. X. Jiang, Nature of exciton transitions in hexagonal boron nitride, Appl. Phys. Lett. 108, 122101 (2016).

[25] T. Q. P. Vuong, G. Cassabois, P. Valvin, V. Jacques, R. Cuscó, L. Artús, and B. Gil, Overtones of interlayer shear modes in the phonon-assisted emission spectrum of hexagonal boron nitride, Phys. Rev. B 95, 045207 (2017).
[26] S. Galambosi, L. Wirtz, J. A. Soininen, J. Serrano, A. Marini, K. Watanabe, T. Taniguchi, S. Huotari, A. Rubio, and K. Hämäläinen, Anisotropic excitonic effects in the energy loss function of hexagonal boron nitride, Phys. Rev. B **83**, 081413 (2011).

[27] G. Fugallo, M. Aramini, J. Koskelo, K. Watanabe, T. Taniguchi, M. Hakala, S. Huotari, M. Gatti, and F. Sottile, Exciton energy-momentum map of hexagonal boron nitride, Phys. Rev. B **92**, 165122 (2015).

[28] R. Schuster, C. Habenicht, M. Ahmad, M. Knupfer, and B. Büchner, Direct observation of the lowest indirect exciton state in the bulk of hexagonal boron nitride, Phys. Rev. B **97**, 041201 (2018).

[29] L. Sponza, H. Amara, F. Ducastelle, A. Loiseau, and C. Attaccalite, Exciton interference in hexagonal boron nitride, Phys. Rev. B **97**, 075121 (2018).

[30] B. Arnaud, S. Lebègue, P. Rabiller, and M. Alouani, Huge Excitonic Effects in Layered Hexagonal Boron Nitride, Phys. Rev. Lett. **96**, 026402 (2006).

[31] L. Wirtz, A. Marini, and A. Rubio, Excitons in Boron Nitride Nanotubes: Dimensionality Effects, Phys. Rev. Lett. **96**, 126104 (2006).

[32] L. Wirtz, A. Marini, M. Grüning, C. Attaccalite, G. Kresse, and A. Rubio, Comment on “Huge Excitonic Effects in Layered Hexagonal Boron Nitride”, Phys. Rev. Lett. **100**, 189701 (2008).

[33] L. Schue, L. Sponza, A. Plaud, H. Bensalah, K. Watanabe, T. Taniguchi, François Ducastelle, A. Loiseau, and J. Barjon, Direct and indirect excitons with high binding energies in hbn, arXiv:1803.03766 [cond-mat.mtrl-sci].

[34] F. Paleari, T. Galvani, H. Amara, F. Ducastelle, A. Molina-Sánchez, and L. Wirtz, Excitons in few-layer hexagonal boron nitride: Davydov splitting and surface localization, 2D Mater. **5**, 045017 (2018).

[35] B. Huang, X. K. Cao, H. X. Jiang, J. Y. Lin, and S.-H. Wei, Origin of the significantly enhanced optical transitions in layered boron nitride, Phys. Rev. B **86**, 155202 (2012).

[36] G. Cappellini, G. Satta, M. Palummo, and G. Onida, Optical properties of bn in cubic and layered hexagonal phases, Phys. Rev. B **64**, 035104 (2001).

[37] X. Blase, A. Rubio, S. G. Louie, and M. L. Cohen, Quasiparticle band structure of bulk hexagonal boron nitride and related systems, Phys. Rev. B **51**, 6868 (1995).

[38] V. L. Solozhenko, G. Will, and F. Elf, Isothermal compression of hexagonal graphite-like boron nitride up to 12 gpa, Solid State Commun. **96**, 1 (1995).

[39] J. Kang, L. Zhang, and S.-H. Wei, A unified understanding of the thickness-dependent bandgap transition in hexagonal two-dimensional semiconductors, J. Phys. Chem. Lett. **7**, 597 (2016).

[40] E. Doni and G. P. Parravicini, Energy bands and optical properties of hexagonal boron nitride and graphite, Il Nuovo Cimento B (1965-1970) **64**, 117 (1969).

[41] P. Cudazzo, L. Sponza, C. Giorgetti, L. Reining, F. Sottile, and M. Gatti, Exciton Band Structure in Two-Dimensional Materials, Phys. Rev. Lett. **116**, 066803 (2016).

[42] P. Y. Yu and M. Cardona, Fundamentals of Semiconductors (Springer, Berlin, Germany, 2010).

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

[44] G. Grosso and G. Parravicini, Solid State Physics, 2nd ed. (Academic Press, San Diego, 2014).

[45] R. W. Boyd, *Nonlinear Optics*, 3rd ed. (Academic Press, San Diego, 2008).

[46] Q. Lin, O. J. Painter, and G. P. Agrawal, Nonlinear optical phenomena in silicon waveguides: Modeling and applications, Opt. Express **15**, 16604 (2007).

[47] I. Souza, J. Iñiguez, and D. Vanderbilt, Dynamics of beryll phase polarization in time-dependent electric fields, Phys. Rev. B **69**, 085106 (2004).

[48] L. F. Richardson and J. A. Gaunt, VIII. The deferred approach to the limit, Philos. Trans. R. Soc., A **226**, 299 (1927).

[49] C. Attaccalite, M. Grüning, and A. Marini, Real-time approach to the optical properties of solids and nanostructures: Time-dependent bethe-salpeter equation, Phys. Rev. B **84**, 245110 (2011).

[50] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. **140**, A1133 (1965).

[51] G. Strinati, Application of the Green's functions method to the study of the optical properties of semiconductors, Nuovo Cimento Rivista Serie **11**, 1 (1988).

[52] G. Grynberg, A. Aspect, and C. Fabre, Introduction to Quantum Optics (Cambridge University Press, 1990).

[53] G. D. Mahan, Theory of two-photon spectroscopy in solids, Phys. Rev. **170**, 825 (1968).

[54] A. Shimizu, Two-photon absorption in quantum-well structures near half the direct band gap, Phys. Rev. B **40**, 1403 (1989).

[55] J. Koskelo, G. Fugallo, M. Hakala, M. Gatti, F. Sottile, and P. Cudazzo, Excitons in van der waals materials: From monolayer to bulk hexagonal boron nitride, Phys. Rev. B **95**, 035125 (2017).

[56] P. Giannozzi et al., QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials, J. Phys.: Condens. Matter **21**, 395502 (2009).

[57] The approach is implemented in the LUMEN extension of the YAMBO many-body code [59]. The source code is available from http://www.attaccalite.com/lumen/.

[58] S. E. Koonin and D. C. Meredith, Computational Physics (Addison-Wesley, Reading, MA, 1990).

[59] A. Marini, C. Hogan, M. Grüning, and D. Varsano, yambo: An ab initio tool for excited state calculations, Comput. Phys. Commun. **180**, 1392 (2009).

[60] T. Cao, M. Wu, and S. G. Louie, Unifying Optical Selection Rules for Excitons in Two Dimensions: Band Topology and Winding Numbers, Phys. Rev. Lett. **120**, 087402 (2018).

[61] X. Zhang, W.-Y. Shan, and D. Xiao, Optical Selection Rule of Excitons in Gapped Chiral Fermion Systems, Phys. Rev. Lett. **120**, 077401 (2018).

[62] C.-J. Kim, L. Brown, M. W. Graham, R. Hovden, R. W. Havener, P. L. McEuen, D. A. Muller, and J. Park, Stacking order dependent second harmonic generation and topological defects in h-bn bilayers, Nano Lett. **13**, 5660 (2013).

[63] Y. Li, Y. Rao, K. F. Mak, Y. You, S. Wang, C. R. Dean, and T. F. Heinz, Probing symmetry properties of few-layer mos2 and h-bn by optical second-harmonic generation, Nano Lett. **13**, 3329 (2013).
[64] M. Grüning and C. Attaccalite, Second harmonic generation in *h*-bn and *mos*₂ monolayers: Role of electron-hole interaction, *Phys. Rev. B* **89**, 081102 (2014).

[65] G. Wang, X. Marie, I. Gerber, T. Amand, D. Lagarde, L. Bouet M. Vidal, A. Balocchi, and B. Urbaszek, Giant Enhancement of the Optical Second-Harmonic Emission of *wse*₂ Monolayers by

Laser Excitation at Exciton Resonances, *Phys. Rev. Lett.* **114**, 097403 (2015).

[66] G. B. Ventura, D. J. Passos, J. M. B. Lopes dos Santos, J. M. V. P. Lopes, and N. M. R. Peres, Gauge covariances and nonlinear optical responses, *Phys. Rev. B* **96**, 035431 (2017).