# Photoluminescent properties of the carbon-dimer defect in hexagonal boron-
nitride: a many-body finite-size cluster approach

Michael Winter  Institute of Theoretical Phyics, Friedrich-Alexander
Universität Erlangen-Nürnberg, Staudtstr. 7/B2, 91058 Erlangen, Germany  Univ.
Grenoble Alpes, CNRS, Inst NEEL, F-38042 Grenoble, France  Manon H. E.
Bousquet  Université de Nantes, CNRS, CEISAM UMR 6230, F-44000 Nantes, France
Denis Jacquemin  Université de Nantes, CNRS, CEISAM UMR 6230, F-44000 Nantes,
France  Ivan Duchemin  Univ. Grenoble Alpes, CEA, IRIG-MEM-L_Sim, 38054
Grenoble, France  Xavier Blase  Univ. Grenoble Alpes, CNRS, Inst NEEL, F-38042
Grenoble, France

###### Abstract

We study the carbon dimer defect in a hexagonal boron-nitride monolayer using
the $GW$ and Bethe-Salpeter many-body perturbation theories within a finite
size cluster approach. While quasiparticle energies converge very slowly with
system size due to missing long-range polarization effects, optical
excitations converge much faster, with a $1/R^{3}$ scaling law with respect to
cluster average radius. We obtain a luminescence zero-phonon energy of 4.36
eV, including significant 0.13 eV zero-point vibrational energy and 0.15 eV
reorganization energy contributions. Inter-layer screening decreases further
the emission energy by about 0.3 eV. These results bring support to the recent
identification of the substitutional carbon dimer as the likely source of the
zero-phonon 4.1 eV luminescence line. Finally, the $GW$ quasiparticle energies
are extrapolated to the infinite h-BN monolayer limit, leading to a predicted
defect HOMO-LUMO photoemission gap of 7.6 eV. Comparison with the optical gap
yields a very large excitonic binding energy of 3 eV for the associated
localized Frenkel exciton.

Ab initio many-body theory; $GW$ formalism; Bethe-Salpeter equation ; defects
; 2D materials

##  I Introduction

There is by now a large experimental [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15, 16] and theoretical [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27,
28, 29, 30, 31, 32, 33, 34, 35, 36, 37] literature focusing on better
understanding point defects in hexagonal boron-nitride (h-BN), following
similar studies on h-BN nanotubes [38]. The very large band gap of the h-BN
host ($\simeq$ 6.1 eV experimental value for the bulk system [39]) allows in
particular the presence of very deep defect-induced occupied and unoccupied
energy levels, yielding optical transitions between defect levels with an
energy that can span over the visible and UV domains. Recently, the
demonstration that such centers could be responsible for room-temperature
single-photon emission sparkled further the interest of the community to
better identify the origin of such bright emission lines [9, 10, 12, 14, 15,
16].

As a paradigmatic defect-associated emission line, the so-called zero-phonon
4.1 eV bright photoemission band in the UV range is very widely observed in
h-BN systems [1, 2, 3, 4, 5, 6, 10, 13]. This sharp emission line appears at
much lower energy than the intense emission at 5.77 eV (215 nm) found in high-
purity h-BN samples [40, 39]. Spatially-resolved cathodoluminescence allowed
to characterize the localized character of the defect-associated excitons
[10], in contrast with the “free” high-purity h-BN excitons that yield a
homogeneous emission all over the sample [41]. The evolution of the
ultraviolet photoluminescent (PL) spectrum as a function of the h-BN synthesis
conditions and reactants led to early identifications of this feature with
either carbon-related defects [1, 2, 5] or nitrogen vacancies [6]. Very
recently, a periodic-boundary constrained density functional theory (DFT)
calculation suggested the carbon dimer defect as the likely source of
emission, with a calculated zero-phonon emission line value of 4.3 eV [32].

At the theoretical level, the electronic properties of defects in periodic
semiconductors or insulators are commonly explored using periodic boundary
conditions (PBC) with large supercells designed to minimize defect
interactions between neighboring cells [42, 43]. While standard density
functional theory (DFT) calculations can be efficiently used for neutral
defects, charged centers require to carefully cancel the long range Coulomb
interaction between cells [42, 43, 44, 45, 46, 47]. However, even in the case
of neutral defects, the study of their electronic properties using accurate
many-body approaches, such as the $GW$ [48] and Bethe-Salpeter [49, 50]
formalisms, run into similar difficulties. In particular, the $GW$ formalism
properly accounts for the effects associated with the addition or removal of a
charge in the system during a photo-emission experiment used to measure
electronic energy levels. As such, Coulomb interactions between localized
“charged” excitations repeated periodically must also be canceled carefully.
This was clearly illustrated in a $GW$ study of h-BN, showing that the
photoemission gap converges very slowly to the monolayer limit as a function
of interlayer spacing [51]. Due to computational cost, the study of defects in
extended systems is still in its infancy at the $GW$ and Bethe-Salpeter levels
[52, 18, 47, 21], and examples of the monitoring of the evolution of a
localized defect electronic properties as a function of unit cell size, or
using a proper scheme to truncate long-range Coulomb effects, remain very
scarce [21].

In the present work, we explore the electronic and optical properties
associated with the neutral carbon-dimer defect in h-BN using many-body
Green’s function perturbation theory (MBPT) in a finite size cluster approach.
Namely, we use the $GW$ [48, 53, 54, 55, 56, 57, 58, 59] and Bethe-Salpeter
equation (BSE) [49, 50, 60, 61, 62, 63, 57, 58, 64] formalisms to calculate
the in-gap quasiparticle energies and lowest optical excitations associated
with the carbon-dimer defect in h-BN flakes [1, 4, 5, 6, 11, 32]. We find that
the vertical BSE excitation energies, together with the excited-state
reorganization energy and zero-point-vibrational energy (ZPVE), converge
relatively fast with system size, leading to a zero-phonon line emission
energy $(E_{ZPL})$ of 4.36 eV in monolayer h-BN. A very similar value is
obtained using TD-DFT with an adequate hybrid exchange-correlations
functional. This is in good agreement with the 4.31 eV value stemming from a
recent theoretical study using an elegant periodic-boundary-conditions
constrained-DFT [32], but not including ZPVE correction calculated here to be
of the order of 0.13 eV. We further demonstrate that multilayer packing may
lower the defect absorption and emission energies by a few tenths of an eV,
bringing theoretical calculations in closer agreement with the 4.1 eV
experimental value [1, 4, 5, 6, 11]. Finally, extrapolating quasiparticle $GW$
data to the infinite monolayer size limit, we predict the HOMO-LUMO gap
associated with the defect to be of the order of 7.6 eV, yielding by
comparison with the optical gap a large excitonic binding energy of about 3 eV
in relation with the very localized nature of the defect-associated
excitation.

##  II  Formalisms and Technical details

The DFT calculations used to generate the needed Kohn-Sham (KS) orbitals that
enter the construction of the $GW$ and Bethe-Salpeter Hamiltonians are
performed with the ORCA code [65] with which we also perform our TD-DFT
calculations. Concerning the many-body Green’s functions perturbation
theories, the $GW$ calculations are performed with the BeDeft (Beyond-DFT)
code [66, 67] that is an ongoing rewriting and extension of the Fiesta code
[68, 69, 70, 71] that we use for Bethe-Salpeter (BSE) calculations. In a
nutshell, the $GW$ self-energy operator, that accounts for exchange and
correlation effects, is first built from input Kohn-Sham
$\\{\varepsilon_{n}^{\mathrm{KS}},\phi_{n}^{\mathrm{KS}}\\}$ one-body
eigenstates, with :

$\displaystyle G({\bf r},{\bf r}^{\prime};\omega)$
$\displaystyle=\sum_{n}\frac{\phi_{n}^{\mathrm{KS}}({\bf
r})\phi_{n}^{\mathrm{KS}}({\bf
r}^{\prime})}{\omega-\varepsilon_{n}^{\mathrm{KS}}+i\eta\times\text{sgn}(\varepsilon_{n}^{\mathrm{KS}}-\mu)}$
$\displaystyle\chi_{0}({\bf r},{\bf r}^{\prime};i\omega)$
$\displaystyle=2\sum_{ja}\frac{{\left(\phi_{j}^{\mathrm{KS}}({\bf
r})\right)}^{*}\phi_{a}^{\mathrm{KS}}({\bf
r}){\left(\phi_{a}^{\mathrm{KS}}({\bf
r}^{\prime})\right)}^{*}\phi_{j}^{\mathrm{KS}}({\bf
r}^{\prime})}{i\omega-(\varepsilon_{a}^{\mathrm{KS}}-\varepsilon_{i}^{\mathrm{KS}})}+c.c.$

where $G$ is the time-ordered Green’s function and $\chi_{0}$ the independent-
electron susceptibility used to build, within the random phase approximation,
the screened Coulomb potential $W=v+v\chi_{0}W$ with $v$ the bare Coulomb
potential. In the preceding equations, $\eta$ is a positive infinitesimal and
$\mu$ the chemical potential. Following quantum chemistry notations, (j/a)
index occupied/empty energy levels. The $GW$ exchange correlation self-energy
operator reads:

$\displaystyle\Sigma^{GW}({\bf r},{\bf r}^{\prime};E)=\frac{i}{2\pi}\int
d\omega\;e^{i\eta\omega}G({\bf r},{\bf r}^{\prime};E+\omega)W({\bf r},{\bf
r}^{\prime};\omega)$

where $W$ is obtained along the real axis by an analytic continuation (AC)
from the $W(i\omega)$ calculated at $n_{\omega}$=14 optimized frequencies
along the imaginary axis following a recently improved AC scheme [72, 66, 67].
The knowledge of $\Sigma^{GW}$ allows to obtained $GW$ quasiparticle energies:

$\displaystyle\varepsilon_{n}^{GW}=\varepsilon_{n}^{\mathrm{KS}}+\left\langle\phi_{n}^{\mathrm{KS}}\left|\Sigma^{GW}\left(\varepsilon_{n}^{GW}\right)-V_{\mathrm{XC}}^{\mathrm{DFT}}\right|\phi_{n}^{\mathrm{KS}}\right\rangle$

We explicitly correct 50 occupied and 50 unoccupied states around the gap,
states at lower/higher energy being corrected following the correction of the
lowest/highest explicitly corrected level.

Following extensive benchmarks on finite size systems [68, 73, 70], we adopt
for accuracy an eigenvalue self-consistent scheme (labeled ev$GW$) where the
corrected quasiparticle energies are reinjected self-consistently in place of
the input KS energy levels to calculate $G$ and $\chi_{0}$. Besides accuracy,
such a self-consistent scheme allows to significantly reduce the impact of the
arbitrariness associated with the choice of the DFT
$V_{\mathrm{XC}}^{\mathrm{DFT}}$ functional used to generate the input
$\\{\varepsilon_{n}^{\mathrm{KS}},\phi_{n}^{\mathrm{KS}}\\}$ eigenstates.
Following again recent benchmarks [68, 73, 70], we select the global hybrid
PBE0 functional for this initial step [74]. Finally, Coulomb integrals are
expressed using Coulomb-fitting resolution-of-the-identity techniques [75, 76,
77] with standard auxiliary basis sets [78] designed for the well-established
cc-pVTZ [79] and augmented _aug_ -cc-pVTZ basis sets used below, or the
auxiliary universal Weigend Coulomb fitting basis set [80] in conjunction with
the smaller 6-311G(d) triple-zeta plus polarization basis set [81] adopted for
exploring large size systems. In the calculation of the susceptibility and
self-energy, all unoccupied levels are considered.

Once the ev$GW$ quasiparticle energies and screened Coulomb potential are
obtained, singlet optical excitation energies $\\{\Omega_{\lambda}\\}$ can be
calculated within the Bethe-Salpeter equation (BSE) formalism that expresses
the two-body electron-hole excitation eigenstates on the basis of the
occupied-to-virtual level transitions:

$\displaystyle\psi^{\mathrm{BSE}}_{\lambda}({\bf r}_{\mathrm{e}},{\bf
r}_{\mathrm{h}})=\sum_{ia}\left[X_{ia}^{\lambda}\phi_{i}({\bf
r}_{\mathrm{h}})\phi_{a}({\bf r}_{\mathrm{e}})+Y_{ia}^{\lambda}\phi_{a}({\bf
r}_{\mathrm{h}})\phi_{i}({\bf r}_{\mathrm{e}})\right]$

where the $Y_{ia}^{\lambda}$ components indicate that we go beyond the Tamm-
Dancoff approximation. These eigenstates solve a Casida-like [82] equation :

$\begin{pmatrix}R&C\\\ -C^{*}&-R^{*}\end{pmatrix}\begin{pmatrix}X^{\lambda}\\\
Y^{\lambda}\end{pmatrix}=\Omega_{\lambda}\begin{pmatrix}X^{\lambda}\\\
Y^{\lambda}\end{pmatrix},$ (1)

where the leading resonant block reads:

$R_{ia,jb}=\left(\varepsilon_{a}^{GW}-\varepsilon_{i}^{GW}\right)\delta_{ij}\delta_{ab}+\kappa
V_{ia,jb}^{\mathrm{x}}-W_{ij,ab}^{\mathrm{d}}$ (2)

with $\kappa$=2 or 0 for singlet or triplet excitations, and :

$\displaystyle V_{ia,jb}^{\mathrm{x}}=(ia|jb)=\int\mathrm{d}{\bf
r}\mathrm{d}{\bf r}^{\prime}\;\phi_{i}({\bf r})\phi_{a}({\bf r})v({\bf r},{\bf
r}^{\prime})\phi_{i}({\bf r}^{\prime})\phi_{b}({\bf r}^{\prime})$
$\displaystyle W_{ij,ab}^{\mathrm{d}}=\int\mathrm{d}{\bf r}\mathrm{d}{\bf
r}^{\prime}\;\phi_{i}({\bf r})\phi_{j}({\bf r})W({\bf r},{\bf
r}^{\prime};\omega=0)\phi_{a}({\bf r}^{\prime})\phi_{b}({\bf r}^{\prime})$

Good convergence can be obtained by including all occupied-to-virtual
transitions such that
$\left(\varepsilon_{a}^{GW}-\varepsilon_{i}^{GW}\right)\leq E_{\mathrm{cut}}$
with $E_{\mathrm{cut}}$ = 30 eV (see SM [83] for convergence study). Otherwise
stated, this will be the cutoff adopted. Further, the convergence with basis
sets size is detailed in the SM [83] evidencing the slower convergence of the
BSE data as compared to TD-DFT.

For computing the zero-point vibrational energy (ZPVE) associated with the
ground and excited states, we take advantage of the analytic DFT and TD-DFT
Hessian available in Gaussian 16.A.03 [84]. For these calculations, performed
at the PBE0 [74] 6-311G(d) level in gas phase, we optimized the geometries
with _tight_ convergence criteria for both the SCF and residual forces, and
select the so-called _ultra-fine_ DFT integration grid. The second-order
coupled-cluster, CC2 [85], and algebraic diagrammatic construction, ADC(2)
[86], calculations used for accuracy validation, are performed with the
Turbomole 7.3 code [87], using various Dunning’s basis sets and applying the
RI approach with the corresponding auxiliary basis. Default Turbomole
parameters and algorithms are applied. The calculations of the vibronic
couplings and Huang-Rhys factors was achieved using the FCClasses-3 package
[88, 89, 90, 91], on the basis of the vibrational data associated with the
ground and excited states as obtained with the Gaussian code. These vibronic
calculations, performed in a time-independent (TI) approach, use the Franck-
Condon (FC) approximation and account for the so-called Duschinsky effects as
well as temperature (298 K), and are based on the most refined adiabatic
Hessian (AH) model to determine the vibronic couplings [90]. We used 108
integrals in the TI scheme, allowing a good convergence (FC recovery $>$ 0.9).

##  III Systems studied

Figure 1: Isocontour representation for the 6-311G(d) PBE0 (a) HOMO and (b)
LUMO levels in the C2@BN138 intermediate flake. In (c,d) the hole-averaged
electron (orange) and electron-averaged hole (purple) densities, respectively,
are given for the lowest singlet BSE excitation.

The finite size h-BN flakes with a central C-C dimer in substitution are
represented in the Supplemental Material (SM) [83]. One of them, of
intermediate size, is shown in Fig. 1. The systems studied comprise from 42 to
240 B,C,N atoms, or 58 to 278 atoms including the H atoms used to avoid
dangling bonds at the borders. These systems are labeled C2@BNX where X counts
the total number of atoms in the systems. The smallest flake is thus named
C2@BN58 and the largest C2@BN278. In such systems, the minimal number of bonds
between a central C atom and an edge B/N atom ranges from 4 to 11.

The geometries of all systems were optimized at the PBE0/6-311G(d) level [74].
The corresponding C-C bond length is found to range between 1.368 Å and 1.369
Å in the smallest and largest flakes, respectively, in close agreement with
the 1.361 Å periodic boundary condition PBEh(0.4)+D3 value of Ref. 32, where
PBEh(0.4) indicates a global hybrid in the PBE0 family with 40 % of exact
exchange (instead of the default 25 %), while the D3 points to the Grimme
empirical correction scheme for dispersion energies. For sake of illustration,
we plot in Fig. 1(a,b) the PBE0 HOMO and LUMO levels associated with the
intermediate size C2@BN138 system. As expected, independently of the system
size, the PBE0 KS HOMO and LUMO levels are localized on the C-C carbon dimer.
The bonding/anti-bonding character of the HOMO/LUMO along the C-C bond is
clearly seen in Fig. 1(a,b).

Anticipating on the analysis of the Bethe-Salpeter optical excitation
energies, we further plot in Fig. 1(c,d) the hole-averaged
$\rho_{\mathrm{e}}({\bf r})$ electron distribution and electron-averaged
$\rho_{\mathrm{h}}({\bf r})$ hole distribution, with e.g.,

$\rho_{\mathrm{e}}({\bf r}_{\mathrm{e}})=\int\mathrm{d}{\bf
r}_{\mathrm{h}}\;\left|\psi^{\mathrm{BSE}}_{\lambda=1}({\bf
r}_{\mathrm{e}},{\bf r}_{\mathrm{h}})\right|^{2}\,.$

The electron and hole localization allows to recognize that the lowest singlet
BSE excitation is mainly of HOMO-LUMO character. The weight $X_{ia}^{\lambda}$
squared of the BSE electron-hole wavefunction on the HOMO-LUMO transition is
found to be consistently of the order of 0.94, rather independently of the
system size beyond the C2@BN86 system, in close agreement with the 0.96 weight
at the TD-PBE0 level. The remaining contributions, even though small, explain
the evolution of the BSE excitation energy with respect to the
$E_{\mathrm{cut}}$ parameter (see SM [83]). While the two C atoms are non-
equivalent, being bonded either to B or N atoms with different
electronegativities, there is no sign of significant charge-transfer character
for the HOMO-LUMO transition, as confirmed below by the large value of the
oscillator strength.

##  IV Benchmark of methods

We first benchmark the accuracy of the present BSE/ev$GW$ scheme, together
with that of standard TD-DFT calculations using various exchange-correlation
functionals, for the present carbon substituted h-BN flakes. The study of
finite size systems allows in particular to use accurate quantum chemistry
wavefunction based methods, such as coupled-cluster (CC) or algebraic
diagrammatic techniques (ADC). A similar comparison between mean-field DFT/TD-
DFT approaches and wavefunction based methods on finite size systems was
already proposed in the case of a few defects in a cluster geometry [26, 31,
35]. Significantly more expensive than BSE and TD-DFT calculations,
calculations with the second-order wavefunction approaches, CC2 and ADC(2),
with a large converged basis set (_aug_ -cc-pVTZ) are only possible for the
smallest C2@BN58 system. Even though limited in size, this small system
presents already well defined defect states, with a lowest singlet excitation
energy within no more than 0.10–0.15 eV of its limit value obtained in larger
flakes (see below).

While very accurate “theoretical best estimates" (TBE), such as those obtained
from third order CC3 calculations [92, 93], cannot be obtained for the present
systems, CC2 and ADC(2) calculations have been shown, thanks to very large
benchmarks on molecular families, to present mean absolute errors of the order
of 0.10 eV to 0.15 eV as compared to TBE or experimental data for various type
of excitations, from localized (Frenkel) to Rydberg or charge-transfer type
transitions [94, 73, 95]. In particular, even though an approximation to CCSD,
CC2 provides a similar level of accuracy for localized (Frenkel) excitations
with a reduced cost scaling [92, 93].

Our data are presented in Table 1. In the converged _aug_ -cc-pVTZ limit, we
observe that the Bethe-Salpeter vertical excitation energy relying on self-
consistent ev$GW$@PBE0 input quasiparticle energies and screened Coulomb
potential, is in close agreement with the CC2 value, even though 0.1 eV larger
than its ADC(2) analogue. In other words, BSE falls within the error bar of
second-order wavefunction schemes. The present results confirm the accuracy of
the BSE/ev$GW$ approach, already demonstrated in large benchmarks on molecular
systems [68, 70] .

Table 1: Lowest singlet ($S_{1}$) excitation energy for the C2@BN058 system. We compare quantum chemistry wavefunction methods [CC2 and ADC(2)], the many-body BSE/ev$GW$@PBE0 method, and TDDFT with the PBE, PBE0 and PBEh(0.4) functionals at the _aug_ -cc-pVTZ level. Energies are in eV and oscillator strength in parenthesis. CC2 | ADC(2) | BSE/evGW@PBE0
---|---|---
4.761 (0.494) | 4.660 (0.468) | 4.765 (0.407)
PBE0 | PBE | PBEh(0.4)
4.716 (0.406) | 4.024 (0.119) | 4.966 (0.443)

Among the TD-DFT calculations, we further observe that the standard TD-PBE0,
with its 25 % ratio of exact exchange, offers an excellent accuracy, standing
in between ADC(2) and CC2 data. TD-PBE0 and BSE/ev$GW$ excitation energies
fall within 50 meV of each other, while their oscillator strengths are also in
excellent agreement. Such values are also in good agreement with the 4.78 eV
CAM-B3LYP energy of Ref. [31] obtained with a minimally augmented cc-pVDZ
basis for a cluster equivalent in size to our C2@BN138 system. In contrast,
TD-PBE calculations yield an excitation energy that is significantly too
small, together with an oscillator strength that is a factor 4 smaller than
that of the CC2 or ADC(2) references. On the contrary, adopting the PBEh(0.4)
functional with 40 % of exact exchange yields an excitation energy that is
somehow too large. The balanced performance of the PBE0 functional is
consistent again with large TD-DFT benchmark on molecular systems, proving
that PBE0 is one of the most accurate functional for localized excited states
[96]. We adopt thus the BSE/ev$GW$ and TD-PBE0 schemes for exploring now the
convergence with system size. Further, since analytic forces in the excited
states have not been developed at the BSE level, reorganization energy and
vibrational-spectrum calculations in the excited state are conducted below at
the TD-PBE0 level.

##  V Vertical excitation energies

We now explore the evolution of the lowest singlet ($S_{1}$) vertical
excitation energy with respect to system size. While TD-DFT and BSE data for
the smallest C2@BN58 system have been explicitly calculated at the _aug_ -cc-
pVTZ level, calculations on the largest flakes have been performed with the
6-311G(d) Pople atomic basis set level, stopping at the C2@BN202 system for
BSE calculations. As shown in the SM [83], the evolution with system size of
the excitation energy is weakly dependent on the basis size, allowing to
safely extrapolate the 6-311G(d) data to their _aug_ -cc-pVTZ value. Further,
the small evolution of these energies beyond the C2@BN86 system insures that
errors in the extrapolation are very limited.

Figure 2: Evolution of the lowest vertical singlet ($S_{1}$) excitation energy
(in eV) with respect to the number of non-H atoms ($N_{\mathrm{at}}$). TD-PBE0
(orange) and BSE (blue) data are presented. Plotted values have been
extrapolated to the _aug_ -cc-pVTZ level (see text and the SM [83]). The
largest flakes data points (out of the shaded area) are fitted by an
asymptotic $\;f(x)=a\cdot x^{3}+c\;$ functional with
$x=1/\sqrt{N_{\mathrm{at}}}$ where $N_{\mathrm{at}}\sim{R^{2}}$ with $R$ the
flake average radius. The infinite radius extrapolated energy value is
indicated in the Inset.

As shown in Fig. 2, both the TD-PBE0 and BSE vertical excitation energies
decrease with system size, converging relatively fast, with the C2@BN138
system being within $\sim$20 meV of the extrapolated infinite size limit. The
decay rate within the TD-PBE0 and BSE formalisms are similar. Several effects
can be invoked to understand the reduction of the excitation energy with
system size. A first possible explanation could hinge on the confinement of
the defect states, resulting in an increase of the HOMO-LUMO gap in the small
size limit. The analysis of the KS HOMO and LUMO levels indicates a 75 meV
reduction of the gap between the C2@BN58 and C2@BN278 systems. Even though not
directly related to the optical gap, this reduction of the Kohn-Sham gap, that
fully accounts for confinement effects, can contribute to the $\sim$0.11 eV
decrease of the optical gap with system size. Confinement is expected to yield
an evolution of the band gap that scales as $1/R^{2}$ where $R$ is the flake
radius.

Figure 3: Size evolution with respect to the smallest C2@BN58 system of the
ev$GW$ HOMO-LUMO gap (blue), the lowest vertical singlet ($S_{1}$) excitation
energy (orange), the
$\langle\psi^{\mathrm{BSE}}_{1}|2V^{\mathrm{x}}|\psi^{\mathrm{BSE}}_{1}\rangle$
(green) and
$\langle\psi^{\mathrm{BSE}}_{1}|W^{\mathrm{d}}|\psi^{\mathrm{BSE}}_{1}\rangle$
(red) matrix elements. Calculations performed with the 6-311G(d) atomic basis
set. Energies are in eV.

Another source of reduction of the optical gap with system size is the
polarization energy, or the increase of environmental screening upon increase
of the system size. It is interesting to analyze in this respect the evolution
with size of the various contributions to the BSE absorption energy, including
in particular the ev$GW$ HOMO-LUMO gap together with the expectation values
$\langle\psi^{\mathrm{BSE}}_{1}|2V^{\mathrm{x}}|\psi^{\mathrm{BSE}}_{1}\rangle$
and
$\langle\psi^{\mathrm{BSE}}_{1}|W^{\mathrm{d}}|\psi^{\mathrm{BSE}}_{1}\rangle$
of the bare and screened Coulomb potentials over the $S_{1}$ BSE eigenstate
(see Eqn. 2). As shown in Fig. 3, the drop of the ev$GW$ gap is much larger
than the 75 meV calculated at the KS level. The $GW$ formalism properly mimics
a photoemission experiment that measures charging energies. Upon charging the
defect levels, the generated Coulomb field strongly polarizes the surrounding
medium. The environment charge rearrangement generates a reaction field that
stabilizes the added charge. This effect, that closes the gap, is a strong
long range polarization effect that slowly converges with system size.

The much faster convergence of the BSE excitation energy can be explained by
the evolution of the screened electron-hole
$\langle\psi^{\mathrm{BSE}}_{1}|W^{\mathrm{d}}|\psi^{\mathrm{BSE}}_{1}\rangle$
term that decays with system size similarly as the ev$GW$ energy levels,
compensating largely the drop in the gap. In other words, the evolution of the
screened Coulomb potential with respect to the size of the polarizable
environment affects in a similar fashion both the electron-electron (exchange-
correlation) interaction that enters the ev$GW$ gap, and the electron-hole
interaction that is subtracted in the BSE Hamiltonian. The unscreened bare
electron
$\langle\psi^{\mathrm{BSE}}_{1}|2V^{\mathrm{x}}|\psi^{\mathrm{BSE}}_{1}\rangle$
contribution varies much less dramatically with system size, confirming that
it is the change in long-range screening that determines the evolution of the
ev$GW$ gap and the counterbalancing electron-hole interaction. Overall, the
optical excitation is a neutral excitation that polarizes the environment much
more weakly as compared to charged excitations.

To analyze more precisely the scaling properties of the optical excitation
energy, we briefly recall the so-called linear-response (LR) and state-
specific (SS) contributions to the energy shift induced by a polarizable
environment on a localized excitation [97, 71, 98]. The localized excitation
is here the dimer-defect excitation and the polarizable environment is the
h-BN host. Assuming for simplicity a pure HOMO-LUMO transition, the LR
contribution to the shift of the absorption energy induced by the environment
reads $\langle\rho_{\mathrm{HL}}|v^{\mathrm{reac}}|\rho_{\mathrm{HL}}\rangle$,
with
$\left|\rho_{\mathrm{HL}}\right\rangle=\left|\phi_{\mathrm{H}}\phi_{\mathrm{L}}\right\rangle$
the product of the HOMO and LUMO one-body orbitals. The reaction field
$v^{\mathrm{reac}}({\bf r},{\bf r}^{\prime})$ describes the potential
generated in ${\bf r}^{\prime}$ by the variation of charge induced in the
environment by a point charge added in ${\bf r}$. In the LR contribution, the
environment reacts to the $\rho_{HL}$ transition codensity, and the LR effect
is closely related to the magnitude of the oscillator strength. In the SS
contribution, the energy shift reads
$\langle\Delta\rho_{\mathrm{ES}}|v^{\mathrm{reac}}|\Delta\rho_{\mathrm{ES}}\rangle/2$,
with $\Delta\rho_{\mathrm{ES}}$ the variation of the total charge density
between the ground-state and the excited state, i.e., the SS contribution is
large when the density reorganization is large, such as in CT excitations.
Both contributions scale as the reaction field that we now show to vary as
$1/R^{3}$ with $R$ the radius of the flake as expected in 2D systems.

Both the $\rho_{\mathrm{HL}}$ and $\Delta\rho_{\mathrm{ES}}$ density
variations have no monopolar component. However, the rather large oscillator
strength associated with the transition indicates a large $\rho_{\mathrm{HL}}$
dipolar component. The $1/r^{3}$ electric field associated with these dipolar
components induces dipoles proportional to $\alpha_{\mathrm{at}}/r^{3}$ on B/N
atoms located at a distance $r$ of the defect, where $\alpha_{\mathrm{at}}$ is
some effective average atomic polarizability. Each induced dipole generates a
reaction potential that stabilizes the defect excited-state dipole.
Integrating over the reaction field generated by all B/N atoms in the flake
leads indeed to an overall reaction potential that scales as $1/R^{3}$, with
$R\simeq\sqrt{N_{\mathrm{at}}}$, taking $N_{\mathrm{at}}$ to be the number of
B/C/N atoms. Fitting with such a functional form the TD-PBE0 data for the 3
largest flakes (Fig. 2, dashed blue line) leads to a fit that smoothly goes
close to the 2 smallest systems for which such an asymptotic behavior may be
less accurate, accounting in particular for the $1/R^{2}$ behavior associated
with confinement. Fitting the 2 largest systems for the BSE data leads to a
fit somehow less satisfactory in the small size limit (Fig. 2, dashed orange
line). We must note that such an asymptotic behavior relies on the assumption
that the defect excitation $\rho_{\mathrm{HL}}$ and $\Delta\rho_{\mathrm{ES}}$
charge variations do not change with system size, a condition that may not be
fully validated in the small size limit.

Overall, the large system size extrapolated TD-PBE0 vertical $S_{1}$
excitation energy is found to be 4.61 eV, within 5 meV of the value obtained
for the largest C2@BN278 system. Similarly, the infinite-system size
extrapolated BSE value of 4.64 eV lies within 10 mev of the value obtained for
the C2@BN202 system. We adopt in the following these extrapolated values. The
present results confirm the conclusions of Ref. 35 that optical excitations
converge rather fast with cluster size.

##  VI Reorganization energies and ZPVE corrections

As a next step, we explore the convergence with system size of the
reorganization energies in the ground and excited states, namely the
differences
$\;\lambda^{GS}=E_{\mathrm{GS}}[R_{\mathrm{ES}}]-E_{\mathrm{GS}}[R_{\mathrm{GS}}]\;$
and
$\;\lambda^{ES}=E_{\mathrm{ES}}[R_{\mathrm{GS}}]-E_{\mathrm{ES}}[R_{\mathrm{ES}}]\;$
where $R_{\mathrm{GS}}$ and $R_{\mathrm{ES}}$ indicate the $S_{0}$ and $S_{1}$
geometries, respectively. Analytic forces in the ground and excited states are
not available within the BSE formalism so that we use TD-PBE0 for this
exploration. As shown in Fig. 4, the reorganization energies slightly decrease
with system size with variations that are of the same magnitude as those
evidenced above for the vertical excitation energy. While an elastic theory
analysis of the deformation energy induced by a point defect in 2D finite-size
systems is beyond the present study, we find that the ground-state
reorganization energy can be accurately fitted by a $1/N_{\mathrm{at}}^{2}$ or
$1/R$ law, where $N_{\mathrm{at}}$ in the number of atoms and $R$ the
effective radius of the flake. We extrapolate the ground-state reorganization
energy to 0.164 eV. For the ES, the reorganization energy is found to be
somehow smaller with a decay law that seems to rely on a larger exponent.
Whatever the reliability of the fit, the extrapolated value and the one
explicitly calculated for the largest C2@BN202 system are within a very few
meV from each other, indicating again that any extrapolation inaccuracy can
hardly affect the conclusions. We adopt the values of 0.16 eV and 0.15 eV for
the reorganization energy in the ground and excited state, respectively. The
smaller value of the reorganization energy in the excited state can be
associated with smaller vibrational energies as compared to the ground state.
This is what we now explore by considering the zero-point vibrational energy
(ZPVE).

Figure 4: Evolution of the ground and excited state reorganization energy (in
eV) at the TD-PBE0/6-311G(d) level. The solid lines are a
1/$N_{\mathrm{at}}^{2}$ fit to the data points. The dashed lines are a "free"
exponent 1/$N_{\mathrm{at}}^{\alpha}$ with $\alpha\simeq 2.08$ for the ground-
state and $\alpha\simeq 3.71$ for the excited-state. Free fit infinite size
extrapolated values amount to 0.164 eV for the ground-state and 0.150 eV for
the excited-state, respectively.

To further account for the ZPVE correction, we perform calculations of the
vibrational spectrum in the ground-state at the ground-state geometry, and in
the $S_{1}$ excited state at the excited-state geometry. The difference in
ZPVE between the excited and ground states amounts to 0.147 eV for the
smallest C2@BN58 system and 0.131 eV in the larger C2@BN86 system. As expected
[99, 100], the vibrational modes in the excited states are “softer”,
consistent with the antibonding character of the LUMO level (see Fig. 1b).
Such vibrational energy differences are rather large yet not exceptional [99,
100]. This is likely related to a significant tensile stress for the C-C bond
embedded in the h-BN lattice with larger bond length, a constraint that is
essentially released upon populating an antibonding state. In fact, the C-C
bond length changes from 1.37 Å to 1.45 Å at the $S_{1}$ excited state
geometry $R_{\mathrm{ES}}$, a bond length that can be compared to 1.44 Å, the
PBE0/6-311G(d) B-N bond length at the center of an undoped pristine h-BN
flake. In the excited state, the C-C bond length is thus much closer to that
of the host network, confirming the release of the tensile stress. Assuming
that the ZPVE correction will further decrease slightly in the limit of larger
systems, we adopt in the following the value of 0.13 eV.

Finally, taking the 4.64 eV calculated BSE vertical absorption energy,
subtracting the 0.15 eV reorganization energy in the $S_{1}$ excited state and
the ZPVE difference contribution of 0.13 eV, we obtain a zero-phonon-line
(ZPL) energy of about 4.36 eV (see summary in Fig. 5). This is close to the
4.31 eV value found in Ref. 32 using constrained DFT calculations in a PBC
geometry, even though this previous value does not contain the ZPVE
contribution.

##  VII Huang-Rhys and Debye-Waller factors

Figure 5:  Symbolic representation of the lowest singlet vertical and zero-
phonon transitions, with the reorganization energy in the ground
($\lambda^{\mathrm{GS}}$) and excited ($\lambda^{\mathrm{ES}}$) states, ZPVE
contribution and Huang-Rhys factor $S$.

We now discuss the Huang-Rhys and Debye-Waller factors as an additional mean
to connect to experimental data [13]. Closely related to the reorganization
energy from the ground-state to the excited-state geometry, the Huang-Rhys
(HR) factor is a measure of the strength of coupling between the excitation
and the vibrational modes. In our case, we computed the HR factors for all
vibrational modes using the adiabatic Hessian vibronic model relying on the
(TD-)PBE0/6-311G(d) vibrational data of the two states. The HR factor reported
below, $S$, is obtained as the sum of the HR terms determined for all
vibrational modes.

The $S$ factor associated with the smallest C2@BN58 system is found to be
$S=1.37$, decreasing to $S=1.23$ in the case of the C2@BN86 system. This 10 %
decrease with system size is consistent with the 10 % decrease of the
reorganization energy between these two systems in the excited state (see Fig.
4). The very limited decrease of the reorganization energy at larger sizes
indicates that the value obtained for the C2@BN86 system is a very good
estimate of the converged $S$ in the large size limit. Our calculated value
falls within the 1–2 range of the experimental estimate [13]. The related
$W_{ZPL}\approx e^{-S}$ Debye-Waller coefficient is found to be 0.29, in good
agreement with the estimated 0.26 experimental value [11].

In the C2@BN58 system, the largest single mode contribution to the total $S$
factor is a clear C-C stretching mode located at 1347 cm-1 (168 meV) in the
excited state, that according to the Duschinsky analysis, couples to several
ground state mode, the strongest being the alike C-C stretching mode at 1644
cm-1 (204 meV) that was shown in Ref. [31] to stand higher in energy as
compared to the optical modes of pristine h-BN. The second and third
contributions originate in couplings between low-frequency breathing modes:
the excited state modes appearing at 338 cm-1 (42 meV) and 236 cm-1 (29 meV)
being respectively coupled to the ground state modes at 254 cm-1 (31 meV) and
242 cm-1 (39 meV). In the C2@BN86 flake, the largest single mode contributions
to $S$ are coming from, on the one hand, a breathing mode (272 cm-1 or 33 meV
in the excited state), and on the other hand, the excited-state C-C stretching
at 1260 cm-1 (156 meV). These modes are represented in the SM [83]. In the
computed vibrationally-resolved emission spectrum (see the SM [83]), the
strongest “stick” contributions also involve these vibrations. We note that
the energy of the stretching and breathing modes can be associated with the
174 meV and 60 meV energy spacing between phonon replica of the 4.08 eV zero-
phonon line as observed in the experimental PL spectrum of Ref. 13.

##  VIII Multi-layer polarization effects

The theoretical $\sim$4.36 eV zero-phonon emission energy calculated for the
defected monolayer system remains larger than the 4.08 eV experimental value.
However, experimental data were not acquired for defected single-layer h-BN
sheet but multilayer systems. As such, environmental screening phenomena are
expected to affect the optical transition energies, an effect that may be
labeled a solvatochromic shift by analogy with the change in absorption or
emission energies for solvated species as compared to the gas phase. Such
environmental screening phenomena were already shown in Refs. [28, 33] to
stabilize the defect states, namely reducing the HOMO-LUMO photoemission gap.
To study this effect in the present case of photoluminescence, we explore
again the optical properties of our smallest C2@BN58 system in a bilayer and
trilayer geometry. For these model systems, the layers are individually
relaxed and then stacked without further relaxation. We adopt the standard AA’
geometry (B atoms on top of N atoms) with a 3.33 Å interlayer spacing [101,
102]. In this geometry, the C-C dimer is located on top of a B-N dimer of
neighboring layers. As shown in Fig 6, the defect-associated HOMO and LUMO
levels remain strongly localized on the defected layer.

Figure 6: Isocontour representation of the defect-associated (a) HOMO and (b)
LUMO levels in a C2@BN58 trilayer system with the defect localized in the
inner layer. The stacking is of AA’ type and the interlayer spacing has been
set to 3.33 Å. In c) and d), the hole-averaged electron density and electron-
averaged hole density for the lowest defect-associated BSE singlet excitation
are represented.

The most salient feature is that the defect-associated BSE absorption energy
is red-shifted by 0.17 eV and 0.27 eV from the monolayer to the bi- and
trilayer systems, respectively, taking the central layer in its ground-state
geometry. The very same 0.27 eV redshift is found for the $S_{1}$ energy of a
trilayer system taking the central layer in its $S_{1}$ excited state
geometry. As such, the effect on the excitation energy of the neighboring
layers seems rather independent of the central layer geometry, suggesting that
the reorganization energy and ZPVE contributions calculated for the defected
monolayer remain good estimates in the multilayer case.

To better analyze this evolution, we selected the central monolayer in its
ground-state geometry and found that between the monolayer and the trilayer
systems, the KS gap closes by 0.17 eV while the $GW$ gap closes by 0.78 eV. As
such, direct hybridization between layers is not negligible, but is marginal
as compared to polarization effects, namely long-range screening by adjacent
layers. This confirms the stabilization of the defect states, namely the
reduction of the ionization potential and increase of the electronic affinity,
induced by environmental screening. The same effect was studied in the case of
h-BN layers on a substrate [33]. Again, interlayer screening reduces as well
the electron-hole interaction, leading to a reduction of the optical
excitation that is significantly smaller than that of the (photoemission) gap.
This remains however a non-negligible effect when attempting to corroborate
theoretical with experimental data, a message that applies to the present
carbon-dimer defect but certainly to any other defects.

We can speculate that most of the screening effects, inducing a redshift of
the optical lines as compared to the defected monolayer system, originates
from the screening by a very few neighbouring layers, tentatively explaining
why experiments performed on bulk or few-layers h-BN provide a very stable ZPL
emission line at about 4.1 eV. Studying the evolution of the emission line
with respect to the number of layers would be prohibitively expensive and lies
beyond the scope of the present study. The main message however is that
calculations of absorption or luminescence energies on monolayer systems
should overestimate the experimental values by a few tenths of an eV. We adopt
the plausible value of 0.3 eV, namely slightly larger than the 0.27 eV value
explicitly calculated for the trilayer system.

Reducing the calculated 4.36 eV ZPL energy for the monolayer system by a
redshift of the order of 0.3 eV brings that calculated value to close
agreement with the 4.08 eV experimental line. Such a close agreement is
certainly fortuitous given the $\simeq$0.1 eV difference on the vertical
absorption energies between high-level methods [BSE, CC2, ADC(2)], the rough
estimate of the neighboring layers screening effect (rounded to 0.3 eV) and
the extrapolations needed to reach the infinite system size and complete basis
set limits. However, the present estimate certainly comforts the carbon-dimer
defect as the origin of this emission line as proposed by the authors of Ref.
32.

##  IX Photoemission gap and excitonic binding energy

We close this study by a brief discussion of the dopant-associated HOMO-LUMO
gap as obtained within the ev$GW$ formalism and the resulting excitonic
binding energy obtained from subtracting the optical gap. Since no
experimental data are available for these quantities, the resulting values
will stand as predictions. More details are available in the SM [83].

Starting with the smallest C2@BN58 system, our most converged _aug_ -cc-pVTZ
ev$GW$ value yields 8.28 eV which is very close to the extrapolated 8.3 eV to
8.4 eV values in the complete basis set (CBS) limit. Adapting the discussion
provided above for the scaling laws associated with long-range polarization
effects, but treating here the response to charged (photoemission)
excitations, the evolution with size of the gap scales as 1/$R^{2}$ (see e.g.
Ref. 103 in the case of molecular 2D systems). Fitting our $GW$ data points
with such a functional form, the extrapolated gap in the infinite system size
limit is found to close by 0.77 eV as compared to the smallest C2@BN58 system,
leading to $\simeq$7.6 eV for the carbon dimer defect HOMO-LUMO gap in an
infinite h-BN monolayer. Subtracting the 4.64 eV vertical optical absorption
energy, one obtains a very large $\simeq$3 eV electron-hole binding energy.
This is much larger than the estimated 0.7-0.8 eV [104, 105, 106] and 1.9 eV -
2.1 eV [51, 107, 108, 104] excitonic binding energy for layered and single-
sheet undefected h-BN systems. The clear localization of the Frenkel electron-
hole pair on the defect, as shown in Fig. 1(c,d), can certainly explain such a
large increase of the excitonic binding energy.

Finally, the knowledge of the defect $GW$ ionization potential and electronic
affinity allows estimating the defect charge transition levels (CTLs) as an
indication of the chemical potential energy range where the neutral defect is
more stable than the charged ones (see e.g. Refs. [23, 34]). Exploiting the
fast convergence with system size of structural reorganization energies, we
estimate from the large C2@BN138 system the reorganization energy from the
neutral to the cationic state geometry to be 0.35 eV [PBE0 6-311G(d) value].
Further, in the infinite size limit (see SM [83], Section VIII), we find the
position of the occupied defect level to be $\simeq$0.64 eV above the
monolayer h-BN occupied band edge. These data bring the (+1/0) CTL at about
1.0 eV above the h-BN top of the valence band, in close agreement with the DFT
calculations of Ref. [32]. Further, with a reorganization energy to the
anionic geometry of 0.25 eV and a defect HOMO-LUMO gap of 7.6 eV, we find the
(0/-1) CTL to lie at about 8.0 eV above the valence band edge. These
calculations confirm the outcomes of Ref. [32] that for the isolated defected
monolayer, the neutral state should be the most stable for a chemical
potential located throughout most of the gap. We will not explore here the
effect of multilayer screening on the CTL positions [34].

##  X Conclusion

We have studied the optoelectronic properties of the recently proposed carbon-
dimer defect in hexagonal boron-nitride using the $GW$ and Bethe-Salpeter
equation many-body Green’s function formalisms in a finite size cluster
(flake) approach. More precisely, we calculate vertical excitation energies
within the BSE formalism, while relaxation in the excited state and zero-point
vibrational energies (ZPVE) are calculated at the TD-PBE0 level. Comparison
with wavefunction CC2 and ADC(2) calculations on a small C2@BN58 atoms system
validated the accuracy of both the BSE formalism and of the TD-PBE0 approach
in the family of TD-DFT techniques.

The study of all relevant quantities with respect to system size allows to
find that the vertical optical absorption energy, the reorganization energy,
and the ZPVE contributions, converge rather quickly with system size, being
close to their extrapolated value for systems containing between 100 and 200
B/C/N atoms. As a result, we find a BSE ZPL emission energy of 4.36 eV for the
defected monolayer, accounting for a 0.15 eV reorganization energy in the
excited state and a significant 0.13 eV ZPVE reduction from the ground to the
excited state. We show further that additional screening associated with
neighboring layers reduce the defect emission energy, with a 0.27 eV redshift
in the simple case of a model trilayer system, bringing the calculated ZPL
energy in close agreement with the 4.08 eV experimental value. Further, the
S$\simeq$1.37 Huang-Rhys and $W_{ZPL}\simeq$0.25 Debye-Waller coefficients are
found to be in excellent agreement with experimental estimates. The present
results thus confirm the identification proposed in Ref. 32 that the carbon-
dimer defect is at the origin of the 4.08 eV ZPL luminescence line. On a more
general ground, the present study suggests that many-body Bethe-Salpeter
calculations, performed on finite size systems containing a unique defect,
stand as a valuable alternative to periodic boundary approaches for the study
of the optical properties of defects in the dilute limit.

In contrast to the neutral optical excitations, the $GW$ quasiparticle
energies (HOMO, LUMO, etc.), namely charged excitations, converge very slowly
with system size. Extrapolation of quasiparticle energies to the infinite size
limit can however be safely obtained with a 1/$R^{2}$ law, with $R$ the radius
of the 2D flake containing the defect. This stems from a simple polarization
scaling law in insulating 2D systems. As a result, the gap is found to close
by $\simeq$0.77 eV from the smallest C2@BN58 system to the infinite size
limit, leading to an ev$GW$ HOMO-LUMO gap of 7.6 eV in the monolayer complete-
basis set and infinite system size limits. Subtracting the vertical absorption
energy, one obtains a very large excitonic binding energy of 3 eV. The very
localized Frenkel character of the excitation can justify that this electron-
hole interaction is much larger than what is found in pristine h-BN systems
where the excitons are much more delocalized.

###### Acknowledgements.

X.B., I.D., and D.J. thank the ANR for financial support in the framework of
the BSE-forces grant (ANR-20-CE29-0005). Calculations have been performed on
the French TGCC Joliot-Curie supercomputing center thanks to a GENCI
allocation under contract noA0090910016. M.H.E.B. and D.J. are indebted to the
CCIPL computational center installed in Nantes for (the always very) generous
allocation of computational time. M.W. thanks the RRZE for allocation of
computational time and support, the ERASMUS+ traineeship program for their
financial support and P. Hansmann as co-supervisor of his master’s period.

## References

  * Katzir  _et al._ [1975] A. Katzir, J. T. Suss, A. Zunger, and A. Halperin, Point defects in hexagonal boron nitride. I. EPR, thermoluminescence, and thermally-stimulated-current measurements, Phys. Rev. B 11, 2370 (1975).
  * Era  _et al._ [1981] K. Era, F. Minami, and T. Kuzuba, Fast luminescence from carbon-related defects of hexagonal boron nitride, J. Lumin. 24-25, 71 (1981).
  * Wu  _et al._ [2004] J. Wu, W.-Q. Han, W. Walukiewicz, J. W. Ager, W. Shan, E. E. Haller, and A. Zettl, Raman spectroscopy and time-resolved photoluminescence of BN and $B_{x}C_{y}N_{z}$ nanotubes, Nano Lett. 4, 647 (2004), https://doi.org/10.1021/nl049862e .
  * Silly  _et al._ [2007] M. G. Silly, P. Jaffrennou, J. Barjon, J.-S. Lauret, F. Ducastelle, A. Loiseau, E. Obraztsova, B. Attal-Tretout, and E. Rosencher, Luminescence properties of hexagonal boron nitride: Cathodoluminescence and photoluminescence spectroscopy measurements, Phys. Rev. B 75, 085205 (2007).
  * Museur  _et al._ [2008] L. Museur, E. Feldbach, and A. Kanaev, Defect-related photoluminescence of hexagonal boron nitride, Phys. Rev. B 78, 155204 (2008).
  * Du  _et al._ [2015] X. Z. Du, J. Li, J. Y. Lin, and H. X. Jiang, The origin of deep-level impurity transitions in hexagonal boron nitride, Appl. Phys. Lett. 106, 021110 (2015), https://doi.org/10.1063/1.4905908 .
  * Meuret  _et al._ [2015] S. Meuret, L. H. G. Tizei, T. Cazimajou, R. Bourrellier, H. C. Chang, F. Treussart, and M. Kociak, Photon bunching in cathodoluminescence, Phys. Rev. Lett. 114, 197401 (2015).
  * Wong  _et al._ [2015] D. Wong, J. Velasco, L. Ju, J. Lee, S. Kahn, H.-Z. Tsai, C. Germany, T. Taniguchi, K. Watanabe, A. Zettl, F. Wang, and M. F. Crommie, Characterization and manipulation of individual defects in insulating hexagonal boron nitride using scanning tunnelling microscopy, Nature Nanotech. 10 (2015).
  * Tran  _et al._ [2016] T. Tran, K. Bray, M. Ford, M. Toth, and I. Aharonovich, Quantum emission from hexagonal boron nitride monolayers, Nature Nanotech 11, 37–41 (2016).
  * Bourrellier  _et al._ [2016] R. Bourrellier, S. Meuret, A. Tararan, O. Stéphan, M. Kociak, L. H. G. Tizei, and A. Zobelli, Bright UV single photon emission at point defects in h-BN, Nano Lett. 16, 4317 (2016), pMID: 27299915, https://doi.org/10.1021/acs.nanolett.6b01368 .
  * Vuong  _et al._ [2016] T. Q. P. Vuong, G. Cassabois, P. Valvin, A. Ouerghi, Y. Chassagneux, C. Voisin, and B. Gil, Phonon-photon mapping in a color center in hexagonal boron nitride, Phys. Rev. Lett. 117, 097402 (2016).
  * Jungwirth  _et al._ [2016] N. R. Jungwirth, B. Calderon, Y. Ji, M. G. Spencer, M. E. Flatté, and G. D. Fuchs, Temperature dependence of wavelength selectable zero-phonon emission from single defects in hexagonal boron nitride, Nano Lett. 16, 6052 (2016), pMID: 27580074, https://doi.org/10.1021/acs.nanolett.6b01987 .
  * Vokhmintsev  _et al._ [2019] A. Vokhmintsev, I. Weinstein, and D. Zamyatin, Electron-phonon interactions in subband excited photoluminescence of hexagonal boron nitride, J. Lumin. 208, 363 (2019).
  * Gottscholl  _et al._ [2020] A. Gottscholl, M. Kianinia, V. Soltamov, S. Orlinskii, G. Mamin, C. Bradac, C. Kasper, K. Krambrock, A. Sperlich, M. Toth, I. Aharonovich, and V. Dyakonov, Initialization and read-out of intrinsic spin defects in a van der waals crystal at room temperature, Nat. Mater. 19, 540 (2020).
  * Hayee  _et al._ [2020] F. Hayee, L. Yu, J. L. Zhang, C. J. Ciccarino, M. Nguyen, A. F. Marshall, I. Aharonovich, J. Vučkovič, P. Narang, T. F. Heinz, and J. A. Dionne, Revealing multiple classes of stable quantum emitters in hexagonal boron nitride with correlated optical and electron microscopy, Nat. Mater. 19, 534 (2020).
  * Mendelson  _et al._ [2021] N. Mendelson, D. Chugh, J. R. Reimers, T. S. Cheng, A. Gottscholl, H. Long, C. J. Mellor, A. Zettl, V. Dyakonov, P. H. Beton, S. V. Novikov, C. Jagadish, H. H. Tan, M. J. Ford, M. Toth, C. Bradac, and I. Aharonovich, Identifying carbon as the source of visible single-photon emission from hexagonal boron nitride, Nat. Matter. 20, 321 (2021).
  * Azevedo  _et al._ [2009] S. Azevedo, J. R. Kaschny, C. M. de Castilho, and F. de Brito Mota, Electronic structure of defects in a boron nitride monolayer, Eur. Phys. J. B 67, 507 (2009).
  * Attaccalite  _et al._ [2011] C. Attaccalite, M. Bockstedte, A. Marini, A. Rubio, and L. Wirtz, Coupling of excitons and defect states in boron-nitride nanostructures, Phys. Rev. B 83, 144115 (2011).
  * Huang and Lee [2012] B. Huang and H. Lee, Defect and impurity properties of hexagonal boron nitride: A first-principles calculation, Phys. Rev. B 86, 245406 (2012).
  * Berseneva  _et al._ [2013] N. Berseneva, A. Gulans, A. V. Krasheninnikov, and R. M. Nieminen, Electronic structure of boron nitride sheets doped with carbon from first-principles calculations, Phys. Rev. B 87, 035404 (2013).
  * Wu  _et al._ [2017a] F. Wu, A. Galatas, R. Sundararaman, D. Rocca, and Y. Ping, First-principles engineering of charged defects for two-dimensional quantum technologies, Phys. Rev. Materials 1, 071001 (2017a).
  * Tawfik  _et al._ [2017] S. A. Tawfik, S. Ali, M. Fronzi, M. Kianinia, T. T. Tran, C. Stampfl, I. Aharonovich, M. Toth, and M. J. Ford, First-principles investigation of quantum emission from hBN defects, Nanoscale 9, 13575 (2017).
  * Wu  _et al._ [2017b] F. Wu, A. Galatas, R. Sundararaman, D. Rocca, and Y. Ping, First-principles engineering of charged defects for two-dimensional quantum technologies, Phys. Rev. Materials 1, 071001 (2017b).
  * Cheng  _et al._ [2017] G. Cheng, Y. Zhang, L. Yan, H. Huang, Q. Huang, Y. Song, Y. Chen, and Z. Tang, A paramagnetic neutral $C_{B}V_{N}$ center in hexagonal boron nitride monolayer for spin qubit application, Comput. Mater. Sci. 129, 247 (2017).
  * Jungwirth and Fuchs [2017] N. R. Jungwirth and G. D. Fuchs, Optical absorption and emission mechanisms of single defects in hexagonal boron nitride, Phys. Rev. Lett. 119, 057401 (2017).
  * Reimers  _et al._ [2018] J. R. Reimers, A. Sajid, R. Kobayashi, and M. J. Ford, Understanding and calibrating density-functional-theory calculations describing the energy and spectroscopy of defect sites in hexagonal boron nitride, J. Chem. Theory Comput. 14, 1602 (2018), pMID: 29412670, https://doi.org/10.1021/acs.jctc.7b01072 .
  * Weston  _et al._ [2018] L. Weston, D. Wickramaratne, M. Mackoit, A. Alkauskas, and C. G. Van de Walle, Native point defects and impurities in hexagonal boron nitride, Phys. Rev. B 97, 214104 (2018).
  * Smart  _et al._ [2018] T. J. Smart, F. Wu, M. Govoni, and Y. Ping, Fundamental principles for calculating charged defect ionization energies in ultrathin two-dimensional materials, Phys. Rev. Materials 2, 124002 (2018).
  * Sajid  _et al._ [2018] A. Sajid, J. R. Reimers, and M. J. Ford, Defect states in hexagonal boron nitride: Assignments of observed properties and prediction of properties relevant to quantum computation, Phys. Rev. B 97, 064101 (2018).
  * Turiansky  _et al._ [2019] M. E. Turiansky, A. Alkauskas, L. C. Bassett, and C. G. Van de Walle, Dangling bonds in hexagonal boron nitride as single-photon emitters, Phys. Rev. Lett. 123, 127401 (2019).
  * Korona and Chojecki [2019] T. Korona and M. Chojecki, Exploring point defects in hexagonal boron-nitrogen monolayers, Int. J. Quant. Chem. 119, e25925 (2019), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.25925 .
  * Mackoit-Sinkevičiené  _et al._ [2019] M. Mackoit-Sinkevičiené, M. Maciaszek, C. G. Van de Walle, and A. Alkauskas, Carbon dimer defect as a source of the 4.1 eV luminescence in hexagonal boron nitride, Appl. Phys. Lett. 115, 212101 (2019), https://doi.org/10.1063/1.5124153 .
  * Wang and Sundararaman [2019] D. Wang and R. Sundararaman, Substrate effects on charged defects in two-dimensional materials, Phys. Rev. Materials 3, 083803 (2019).
  * Wang and Sundararaman [2020] D. Wang and R. Sundararaman, Layer dependence of defect charge transition levels in two-dimensional materials, Phys. Rev. B 101, 054103 (2020).
  * Reimers  _et al._ [2020] J. R. Reimers, J. Shen, M. Kianinia, C. Bradac, I. Aharonovich, M. J. Ford, and P. Piecuch, Photoluminescence, photophysics, and photochemistry of the $V_{B}^{-}$ defect in hexagonal boron nitride, Phys. Rev. B 102, 144105 (2020).
  * Sajid and Thygesen [2020] A. Sajid and K. S. Thygesen, $V_{N}C_{B}$ defect as source of single photon emission from hexagonal boron nitride, 2D Materials 7, 031007 (2020).
  * Jara  _et al._ [2021] C. Jara, T. Rauch, S. Botti, M. A. L. Marques, A. Norambuena, R. Coto, J. E. Castellanos-Águila, J. R. Maze, and F. Munoz, First-principles identification of single photon emitters based on carbon clusters in hexagonal boron nitride, J. Phys. Chem. A 125, 1325 (2021), pMID: 33554602, https://doi.org/10.1021/acs.jpca.0c07339 .
  * Attaccalite  _et al._ [2013] C. Attaccalite, L. Wirtz, A. Marini, and A. Rubio, Efficient gate-tunable light-emitting device made of defective boron nitride nanotubes: from ultraviolet to the visible, Sci. Rep. 3 (2013).
  * Cassabois  _et al._ [2016] G. Cassabois, P. Valvin, and B. Gil, Hexagonal boron nitride is an indirect bandgap semiconductor, Nature Photon. 10 (2016).
  * Watanabe  _et al._ [2004] K. Watanabe, T. Taniguchi, and H. Kanda, Direct-bandgap properties and evidence for ultraviolet lasing of hexagonal boron nitride single crystal, Nature Mater. 3 (2004).
  * Jaffrennou  _et al._ [2007] P. Jaffrennou, J. Barjon, J.-S. Lauret, B. Attal-Trétout, F. Ducastelle, and A. Loiseau, Origin of the excitonic recombinations in hexagonal boron nitride by spatially resolved cathodoluminescence spectroscopy, J. Appl. Phys. 102, 116102 (2007), https://doi.org/10.1063/1.2821413 .
  * Leslie and Gillan [1985] M. Leslie and N. J. Gillan, The energy and elastic dipole tensor of defects in ionic crystals calculated by the supercell method, J. Phys., C, Solid State Phys. 18, 973 (1985).
  * Makov and Payne [1995] G. Makov and M. C. Payne, Periodic boundary conditions in ab initio calculations, Phys. Rev. B 51, 4014 (1995).
  * Freysoldt  _et al._ [2009] C. Freysoldt, J. Neugebauer, and C. G. Van de Walle, Fully ab initio finite-size corrections for charged-defect supercell calculations, Phys. Rev. Lett. 102, 016402 (2009).
  * Taylor and Bruneval [2011] S. E. Taylor and F. Bruneval, Understanding and correcting the spurious interactions in charged supercells, Phys. Rev. B 84, 075155 (2011).
  * Freysoldt  _et al._ [2014] C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti, and C. G. Van de Walle, First-principles calculations for point defects in solids, Rev. Mod. Phys. 86, 253 (2014).
  * Chen and Pasquarello [2015] W. Chen and A. Pasquarello, First-principles determination of defect energy levels through hybrid density functionals and $GW$, J. Phys. : Condens. Matter 27, 133202 (2015).
  * Hedin [1965] L. Hedin, New method for calculating the one-particle green’s function with application to the electron-gas problem, Phys. Rev. 139, A796 (1965).
  * Salpeter and Bethe [1951] E. E. Salpeter and H. A. Bethe, A relativistic equation for bound-state problems, Phys. Rev. 84, 1232 (1951).
  * Strinati [1982] G. Strinati, Dynamical shift and broadening of core excitons in semiconductors, Phys. Rev. Lett. 49, 1519 (1982).
  * Wirtz  _et al._ [2006] L. Wirtz, A. Marini, and A. Rubio, Excitons in boron nitride nanotubes: Dimensionality effects, Phys. Rev. Lett. 96, 126104 (2006).
  * Rinke  _et al._ [2009] P. Rinke, A. Janotti, M. Scheffler, and C. G. Van de Walle, Defect formation energies without the band-gap problem: Combining density-functional theory and the $GW$ approach for the silicon self-interstitial, Phys. Rev. Lett. 102, 026402 (2009).
  * Strinati  _et al._ [1980] G. Strinati, H. J. Mattausch, and W. Hanke, Dynamical correlation effects on the quasiparticle bloch states of a covalent crystal, Phys. Rev. Lett. 45, 290 (1980).
  * Hybertsen and Louie [1986] M. S. Hybertsen and S. G. Louie, Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies, Phys. Rev. B 34, 5390 (1986).
  * Godby  _et al._ [1988] R. W. Godby, M. Schlüter, and L. J. Sham, Self-energy operators and exchange-correlation potentials in semiconductors, Phys. Rev. B 37, 10159 (1988).
  * von der Linden and Horsch [1988] W. von der Linden and P. Horsch, Precise quasiparticle energies and hartree-fock bands of semiconductors and insulators, Phys. Rev. B 37, 8351 (1988).
  * Onida  _et al._ [2002] G. Onida, L. Reining, and A. Rubio, Electronic excitations: density-functional versus many-body green’s-function approaches, Rev. Mod. Phys. 74, 601 (2002).
  * Martin  _et al._ [2016] R. Martin, L. Reining, and D. Ceperley,  _Interacting Electrons: Theory and Computational Approaches_ (Cambridge University Press, 2016).
  * Golze  _et al._ [2019] D. Golze, M. Dvorak, and P. Rinke, The $GW$ compendium: A practical guide to theoretical photoemission spectroscopy, Front. Chem. 7, 377 (2019).
  * Strinati [1984] G. Strinati, Effects of dynamical screening on resonances at inner-shell thresholds in semiconductors, Phys. Rev. B 29, 5718 (1984).
  * Benedict  _et al._ [1998] L. X. Benedict, E. L. Shirley, and R. B. Bohn, Optical absorption of insulators and the electron-hole interaction: An ab initio calculation, Phys. Rev. Lett. 80, 4514 (1998).
  * Rohlfing and Louie [1998] M. Rohlfing and S. G. Louie, Electron-hole excitations in semiconductors and insulators, Phys. Rev. Lett. 81, 2312 (1998).
  * Albrecht  _et al._ [1998] S. Albrecht, L. Reining, R. Del Sole, and G. Onida, Ab initio calculation of excitonic effects in the optical spectra of semiconductors, Phys. Rev. Lett. 80, 4510 (1998).
  * Ping  _et al._ [2013] Y. Ping, D. Rocca, and G. Galli, Electronic excitations in light absorbers for photoelectrochemical energy conversion: first principles calculations based on many body perturbation theory, Chem. Soc. Rev. 42, 2437 (2013).
  * Neese [2012] F. Neese, The orca program system, WIREs Comput. Mol. Sci. 2, 73 (2012).
  * Duchemin and Blase [2020] I. Duchemin and X. Blase, Robust analytic-continuation approach to many-body $GW$ calculations, J. Chem. Theory Comput. 16, 1742 (2020), pMID: 32023052, https://doi.org/10.1021/acs.jctc.9b01235 .
  * Duchemin and Blase [2021] I. Duchemin and X. Blase, Cubic-scaling all-electron $GW$ calculations with a separable density-fitting space–time approach, J. Chem. Theory Comput. 17, 2383 (2021), pMID: 33797245, https://doi.org/10.1021/acs.jctc.1c00101 .
  * Jacquemin  _et al._ [2015a] D. Jacquemin, I. Duchemin, and X. Blase, Benchmarking the bethe–salpeter formalism on a standard organic molecular set, J. Chem. Theory Comput. 11, 3290 (2015a), pMID: 26207104, https://doi.org/10.1021/acs.jctc.5b00304 .
  * Li  _et al._ [2016] J. Li, G. D’Avino, I. Duchemin, D. Beljonne, and X. Blase, Combining the many-body $GW$ formalism with classical polarizable models: Insights on the electronic structure of molecular solids, J. Phys. Chem. Lett. 7, 2814 (2016), pMID: 27388926, https://doi.org/10.1021/acs.jpclett.6b01302 .
  * Jacquemin  _et al._ [2017] D. Jacquemin, I. Duchemin, and X. Blase, Is the bethe–salpeter formalism accurate for excitation energies? comparisons with TD-DFT, CASPT2, and EOM-CCSD, J. Phys. Chem. Lett. 8, 1524 (2017), pMID: 28301726, https://doi.org/10.1021/acs.jpclett.7b00381 .
  * Duchemin  _et al._ [2018] I. Duchemin, C. A. Guido, D. Jacquemin, and X. Blase, The bethe–salpeter formalism with polarisable continuum embedding: Reconciling linear-response and state-specific features, Chem. Sci. 9, 4430 (2018).
  * Friedrich [2019] C. Friedrich, Tetrahedron integration method for strongly varying functions: Application to the $GT$ self-energy, Phys. Rev. B 100, 075142 (2019).
  * Jacquemin  _et al._ [2015b] D. Jacquemin, I. Duchemin, and X. Blase, 0–0 energies using hybrid schemes: Benchmarks of TD-DFT, CIS(D), ADC(2), CC2, and BSE/$GW$ formalisms for 80 real-life compounds, J. Chem. Theory Comput. 11, 5340 (2015b), pMID: 26574326, https://doi.org/10.1021/acs.jctc.5b00619 .
  * Adamo and Barone [1999] C. Adamo and V. Barone, Toward reliable density functional methods without adjustable parameters: the PBE0 model, J. Comp. Phys. 110, 6158 (1999).
  * Whitten [1973] J. L. Whitten, Coulombic potential energy integrals and approximations, J. Chem. Phys. 58, 4496 (1973), https://doi.org/10.1063/1.1679012 .
  * Ren  _et al._ [2012] X. Ren, P. Rinke, V. Blum, J. Wieferink, A. Tkatchenko, A. Sanfilippo, K. Reuter, and M. Scheffler, Resolution-of-identity approach to hartree–fock, hybrid density functionals, RPA, MP2 and $GW$ with numeric atom-centered orbital basis functions, New J. Phys. 14, 053020 (2012).
  * Duchemin  _et al._ [2017] I. Duchemin, J. Li, and X. Blase, Hybrid and constrained resolution-of-identity techniques for coulomb integrals, J. Chem. Theory Comput. 13, 1199 (2017), pMID: 28094983, https://doi.org/10.1021/acs.jctc.6b01215 .
  * Weigend  _et al._ [2002] F. Weigend, A. Köhn, and C. Hättig, Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations, J. Chem. Phys. 116, 3175 (2002), https://doi.org/10.1063/1.1445115 .
  * Dunning [1989] T. H. Dunning, Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen, J. Chem. Phys. 90, 1007 (1989), https://doi.org/10.1063/1.456153 .
  * Weigend [2008] F. Weigend, Hartree–fock exchange fitting basis sets for H to Rn, J. Comp. Chem. 29, 167 (2008), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.20702 .
  * Krishnan  _et al._ [1980] R. Krishnan, J. S. Binkley, R. Seeger, and J. A. Pople, Self-consistent molecular orbital methods. XX. A basis set for correlated wave functions, J. Chem. Phys. 72, 650 (1980), https://doi.org/10.1063/1.438955 .
  * Casida [1995] M. E. Casida, Time-dependent density functional response theory for molecules, in  _Recent Advances in Density Functional Methods_ (1995) pp. 155–192, https://www.worldscientific.com/doi/pdf/10.1142/9789812830586_0005 .
  * [83] See Supplemental Material at [Link to be added upon publication] for the representation of the C2@BNX structures, the convergence with basis set and atomic structure sizes of TD-DFT and BSE calculations, the complete basis set (CBS) analysis of the smallest flake $GW$ energy levels, the extrapolation to infinite system size of the $GW$ HOMO-LUMO gap, the illustration of the vibrational modes contributing to the Huang-Rhys factor, the theoretical emission spectrum at room temperature and finally a study of the defect levels position with respect to the h-BN monolayer valence band edge.
  * Frisch  _et al._ [2016] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, Gaussian 16 Revision A.03 (2016), gaussian Inc. Wallingford CT.
  * Koch  _et al._ [1995] H. Koch, O. Christiansen, P. Jorgensen, and J. Olsen, Excitation energies of BH, CH2 and Ne in full configuration interaction and the hierarchy CCS, CC2, CCSD and CC3 of coupled cluster models, CPL 244, 75 (1995).
  * Dreuw and Wormit [2015] A. Dreuw and M. Wormit, The algebraic diagrammatic construction scheme for the polarization propagator for the calculation of excited states, WIREs Comput. Mol. Sci. 5, 82 (2015).
  * Tur [2018] Turbomole v7.3 2018, a development of university of karlsruhe and forschungszentrum karlsruhe gmbh, 1989-2007, turbomole gmbh, since 2007; available from http://www.turbomole.com (accessed 13 june 2016). (2018).
  * Santoro  _et al._ [2007a] F. Santoro, R. Improta, A. Lami, J. Bloino, and V. Barone, Effective method to compute franck-condon integrals for optical spectra of large molecules in solution, J. Comp. Phys. 126, 084509 (2007a).
  * Santoro  _et al._ [2007b] F. Santoro, R. Improta, A. Lami, J. Bloino, and V. Barone, Effective method to compute vibrationally resolved optical spectra of large molecules at finite temperature in the gas phase and in solution, J. Comp. Phys. 126, 184102 (2007b).
  * Santoro and Jacquemin [2016] F. Santoro and D. Jacquemin, Going beyond the vertical approximation with time-dependent density functional theory, WIREs Comput. Mol. Sci. 6, 460 (2016).
  * Cerezo and Santoro [2020] J. Cerezo and F. Santoro, Fcclasses 3.0, http://www.pi.iccom.cnr.it/fcclasses (2020).
  * Koch  _et al._ [1997] H. Koch, O. Christiansen, P. Jorgensen, A. M. Sanchez de Merás, and T. Helgaker, The cc3 model: An iterative coupled cluster approach including connected triples, J. Chem. Phys. 106, 1808 (1997).
  * Véril  _et al._ [2021] M. Véril, A. Scemama, M. Caffarel, F. Lipparini, M. Boggio-Pasqua, D. Jacquemin, and P.-F. Loos, Questdb: a database of highly-accurate excitation energies for the electronic structure community, WIREs Comput. Mol. Sci. 11, e1517 (2021).
  * Winter  _et al._ [2013] N. O. C. Winter, N. K. Graf, S. Leutwyler, and C. Hättig, Benchmarks for 0–0 transitions of aromatic organic molecules: DFT/B3LYP, ADC(2), CC2, SOS-CC2 and SCS-CC2 compared to high-resolution gas-phase data, Phys. Chem. Chem. Phys. 15, 6623 (2013).
  * Oruganti  _et al._ [2016] B. Oruganti, C. Fang, and B. Durbeej, Assessment of a composite CC2/DFT procedure for calculating 0–0 excitation energies of organic molecules, Mol. Phys. 114, 3448 (2016), http://dx.doi.org/10.1080/00268976.2016.1235736 .
  * Laurent and Jacquemin [2013] A. D. Laurent and D. Jacquemin, TD-DFT benchmarks: A review, Int. J. Quantum Chem. 113, 2019 (2013).
  * Cammi  _et al._ [2005] R. Cammi, S. Corni, B. Mennucci, and J. Tomasi, Electronic excitation energies of molecules in solution: State specific and linear response methods for nonequilibrium continuum solvation models, J. Chem. Phys. 122, 104513 (2005), https://doi.org/10.1063/1.1867373 .
  * Guido and Caprasecca [2019] C. A. Guido and S. Caprasecca, On the description of the environment polarization response to electronic transitions, Int. J. Quantum Chem. 119, e25711 (2019), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.25711 .
  * Goerigk  _et al._ [2009] L. Goerigk, J. Moellmann, and S. Grimme, Computation of accurate excitation energies for large organic molecules with double-hybrid density functionals, Phys. Chem. Chem. Phys. 11, 4611 (2009).
  * Jacquemin  _et al._ [2012] D. Jacquemin, A. Planchat, C. Adamo, and B. Mennucci, A td-dft assessment of functionals for optical 0-0 transitions in solvated dyes, J. Chem. Theory Comput. 8, 2359 (2012).
  * Lin  _et al._ [2010] Y. Lin, T. V. Williams, and J. W. Connell, Soluble, exfoliated hexagonal boron nitride nanosheets, J. Phys. Chem. Lett. 1, 277 (2010), https://doi.org/10.1021/jz9002108 .
  * Shi  _et al._ [2010] Y. Shi, C. Hamsen, X. Jia, K. K. Kim, A. Reina, M. Hofmann, A. L. Hsu, K. Zhang, H. Li, Z.-Y. Juang, M. S. Dresselhaus, L.-J. Li, and J. Kong, Synthesis of few-layer hexagonal boron nitride thin film by chemical vapor deposition, Nano Lett. 10, 4134 (2010), pMID: 20812716, https://doi.org/10.1021/nl1023707 .
  * D’Avino  _et al._ [2016] G. D’Avino, L. Muccioli, F. Castet, C. Poelking, D. Andrienko, Z. G. Soos, J. Cornil, and D. Beljonne, Electrostatic phenomena in organic semiconductors: fundamentals and implications for photovoltaics, J. Phys.: Cond. Mat. 28, 433002 (2016).
  * Paleari  _et al._ [2018] 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 Materials 5, 045017 (2018).
  * Aggoune  _et al._ [2018] W. Aggoune, C. Cocchi, D. Nabok, K. Rezouali, M. A. Belkhir, and C. Draxl, Dimensionality of excitons in stacked van der waals materials: The example of hexagonal boron nitride, Phys. Rev. B 97, 241114 (2018).
  * Kolos and Karlický [2019] M. Kolos and F. Karlický, Accurate many-body calculation of electronic and optical band gap of bulk hexagonal boron nitride, Phys. Chem. Chem. Phys. 21, 3999 (2019).
  * Galvani  _et al._ [2016] T. Galvani, F. Paleari, H. P. C. Miranda, A. Molina-Sánchez, L. Wirtz, S. Latil, H. Amara, and F. m. c. Ducastelle, Excitons in boron nitride single layer, Phys. Rev. B 94, 125303 (2016).
  * Cudazzo  _et al._ [2016] 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).