††thanks: These two authors contributed equally.††thanks: These two authors
contributed equally.

# Hybridized defects in solid-state materials as artificial molecules

Derek S. Wang  Harvard John A. Paulson School of Engineering and Applied
Sciences, Harvard University, Cambridge, MA 02138, USA  Christopher J.
Ciccarino  Harvard John A. Paulson School of Engineering and Applied Sciences,
Harvard University, Cambridge, MA 02138, USA  Department of Chemistry and
Chemical Biology, Harvard University, Cambridge, MA 02138, USA  Johannes Flick
Center for Computational Quantum Physics, Flatiron Institute, New York, NY
10010, USA  Prineha Narang<EMAIL_ADDRESS>Harvard John A. Paulson
School of Engineering and Applied Sciences, Harvard University, Cambridge, MA
02138, USA

###### Abstract

Two-dimensional materials can be crafted with structural precision approaching
the atomic scale, enabling quantum defects-by-design. These defects are
frequently described as “artificial atoms” and are emerging optically-
addressable spin qubits. However, interactions and coupling of such artificial
atoms with each other, in the presence of the lattice, is remarkably
underexplored. Here we present the formation of “artificial molecules” in
solids, introducing a new degree of freedom in control of quantum
optoelectronic materials. Specifically, in monolayer hexagonal boron nitride
as our model system, we observe configuration- and distance-dependent
dissociation curves and hybridization of defect orbitals within the bandgap
into bonding and antibonding orbitals, with splitting energies ranging from
$\sim$ 10 meV to nearly 1 eV. We calculate the energetics of cis and trans
out-of-plane defect pairs CH${}_{\textrm{B}}$-CH${}_{\textrm{B}}$ against an
in-plane defect pair C${}_{\textrm{B}}$-C${}_{\textrm{B}}$ and find that in-
plane defect pair interacts more strongly than out-of-plane pairs. We
demonstrate an application of this chemical degree of freedom by varying the
distance between C${}_{\textrm{B}}$ and V${}_{\textrm{N}}$ of
C${}_{\textrm{B}}$V${}_{\textrm{N}}$ and observe changes in the predicted peak
absorption wavelength from the visible to the near-infrared spectral band. We
envision leveraging this chemical degree of freedom of defect complexes to
precisely control and tune defect properties towards engineering robust
quantum memories and quantum emitters for quantum information science.

††preprint: APS/123-QED

##  I Introduction

Point defects in both 2D and 3D solids have seen wide applicability in quantum
information science [1, 2, 3, 4, 5, 6], especially as quantum memories,
because they have the potential to combine favorable coherence and non-
classical emission properties of isolated atoms [7, 8] with the scalability
and stability of solid-state technologies [5, 4, 9, 10]. These defects often
introduce localized electronic states that can be effectively decoupled from
the host lattice, thereby creating an “artificial atom”. Such artificial atom
systems can include simple substitutional or vacancy defects, as well as more
complicated defect complexes [11, 12, 13] that are comprised of multiple
imperfections but effectively act as a single defect. To harness the
electronic and optical properties of these defects, defects can be coupled to
external fields, including electric, magnetic, and also strain, as well as to
waveguides and cavity environments [14, 15, 16, 17, 18, 19, 20, 21]. Despite
these advances, demands to the properties of defect systems are ever-
increasing in complexity, such as specific level structures for emission of
entangled photonic states [6, 22] or implementation of multi-qubit photonic
gates [23]. Therefore, the need for defect centers with engineered electronic,
optical, and spin properties has motivated the search for novel defect centers
[24, 25, 9, 26]. Theory and computationally-led searches have proved
particularly valuable here, resulting, for instance, in the discovery of new
group IV and prediction of group III emitters in diamond [27, 28, 29, 30, 31].

Discoveries in synthesis and patterning of atomically-architectured materials
have so far yielded a diverse portfolio of photonic and optoelectronic quantum
materials, including a variety of two-dimensional layered architectures that
can be crafted with structural precision approaching the atomic scale,
enabling quantum defects-by-design. Inspired by advances in experimental
techniques for atom-by-atom manipulation, creation, and observation of
individual defects [32, 33, 34, 35, 36, 37, 18, 24, 38, 39, 40, 41, 42], we
add the crucial chemical degree of freedom to the state-of-the-art in quantum
defects by demonstrating how these “artificial atoms” become “artificial
molecules” when placed within a few angstroms of each other and investigate
the properties of such coupled defects. We study pairs of substitutional
defects within a monolayer of hexagonal boron nitride (hBN), which has a wide
bandgap well-suited for hosting well-defined defect orbitals.

In our calculations, we observe features analogous to molecules, such as the
formation of delocalized, bonding and anti-bonding “molecular” orbitals
arising from the hybridization of the localized, “atomic” orbitals of the
single defects with splitting energies ranging from $\sim 10$ meV to nearly 1
eV. Notably, we determine the energetics of paired in-plane defects and cis
and trans paired out-of-plane defects, showing that the studied paired in-
plane defects interact more strongly than out-of-plane ones. With these
energies, we generate dissociation curves, and find that they are not merely a
function of the inter-defect distance, as is the case in chemical bonding
between centrosymmetric atoms, but also the configuration of the defects
within the host lattice. Importantly, we find that at these length scales, the
host lattice is not only a passive dielectric that screens interactions
between defects in complexes, but actively participates in structure
relaxation and determines configuration-dependent interactions.

Armed with this physical understanding, we demonstrate how this chemical
degree of freedom can be used to tune quantum emitters. Specifically, we vary
the defect-defect distance in C${}_{\textrm{B}}$V${}_{\textrm{N}}$ comprised
of C${}_{\textrm{B}}$ and V${}_{\textrm{N}}$, a previously-proposed candidate
for emission of visible light from hBN, and observe changes in the predicted
peak absorption wavelength from the visible to the near-infrared spectral
range. We anticipate that further exploring this chemical degree of freedom of
defect complexes with higher order predictive methods will enable new
approaches to precisely engineer defect properties as quantum memories and
quantum emitters for quantum information science.

Figure 1: Two defects $A_{\alpha}$-$B_{\beta}@X$ are placed in (a) a large
nanoflake of hBN, where boron is gray and nitrogen is blue. (b) The first
defect substituent $A$ replaces the $\alpha$ boron in position 0, and the
second defect substituent $B$ replaces the $\beta$ boron atom at position $X$
from 1-4. Out-of-plane and in-plane views of three model double defects (c)
cis CHB-CHB@1, (d) trans CHB-CHB@1, and (e) CB-CB@1, where hydrogen is bright
blue and carbon is yellow.

##  II Results and Discussion

We model the interactions between two defects spaced up to four unit cells
apart in a nanoflake of monolayer hBN large enough to simulate bulk behavior
and terminated with hydrogen atoms. For delineating the double defect
structure, we adopt the notation to be generally of the form
$A_{\alpha}$-$B_{\beta}@X$. Here, substituents $A$ and $B$ correspond to the
defect substituent or vacancy being introduced. The indices $\alpha$ and
$\beta$ represent the atomic sites on which the substituents are placed, which
in our case is always boron B for $\alpha$ and either B or nitrogen N for
$\beta$. Finally, $X$ corresponds to the location of the second defect marked
in the schematic in Fig. 1(b). To theoretically model these systems, we use
the pseudopotential, real-space density functional theory method Octopus [43,
44, 45] (see Methods section for computational details). We choose a real-
space as opposed to a periodic density functional approach to avoid artificial
interactions across neighboring supercells which could obscure the defect-
defect interactions of interest. Furthermore, it has been shown that real-
space electronic structure calculations on hexagonal boron nitride (hBN)
nanoflakes can be extrapolated onto periodic calculations [46, 47]. The size
of the nanoflake, around 220 atoms in total, is chosen such that there are at
least 4 unit cells of hBN beyond the boundaries of the localized defect
orbitals before termination with passivating hydrogen atoms. The calculated
energy difference between the valence and conduction bands is found to be
$E_{\rm bg}=4.28$ eV $\pm 1\%$ for all calculations presented here.

To develop a general understanding of the artificial molecule picture, we
first study substituents where $A$ and $B$ are either a carbon atom C or a
carbon bonded to a hydrogen CH. Both present a single defect orbital within
the band gap, but upon relaxation the $\rm CH_{B}$ defect adopts an out-of-
plane trigonal pyramidal geometry, while the $\rm C_{B}$ defect adopts an in-
plane trigonal planar geometry. From these two single defects $\rm CH_{B}$ and
$\rm C_{B}$, we compare three distinct pairs of defects. Using the
nomenclature defined above, these are cis CHB-CH${}_{\rm B}@X$, trans CHB-
CH${}_{\rm B}@X$, and CB-C${}_{\rm B}@X$. Here, $X=1$ to 4 corresponds to
lattice sites; see Fig. 1(a). There is always one defect located at $X=0$. The
corresponding in-plane distances for the defect pairs for $X=1$ to 4 are 2.62,
4.54, 5.25 and 6.93 Å, respectively. The out-of-plane and in-plane geometries
of these three distinct pairs of defects are shown in Fig. 1(c)-(e). These
defect pairs allow us to generate qualitative understandings of wavefunction
overlap in both in- and out-of-plane molecular configurations.

Figure 2: Orbital of single defects within the band gap for (a) and (d) CHB
and (g) CB. Bonding (middle) and anti-bonding (right) defect orbitals,
corresponding to lower and higher energy, within the bandgap for the nearest
defect-defect positions calculated for (b)-(c) cis CHB-CH${}_{\rm B}@1$,
(e)-(f) trans CHB-CH${}_{\rm B}@1$, and (h)-(i) CB-C${}_{\rm B}@1$. Red and
purple colors of the wavefunction isosurfaces correspond to positive and
negative values, respectively.

The single defects CHB and CB discussed above each contribute one orbital in
the band gap between the highest occupied (“valence band”) and lowest
unoccupied delocalized Kohn-Sham states (“conduction band”), plotted in Fig.
2(a)/(d) and Fig. 2(g), respectively. When a second defect is added, we find
two orbitals collectively that are separated in energy from each other. In the
spirit of molecular orbital theory, we denote the lower energy one the
“bonding orbital” and the higher one the “antibonding orbital.” Figure 2
compares the orbital character of these defect orbitals of the most strongly
interacting configuration of each defect pair, where $X=1$. For all three
defect pairs, both the bonding and antibonding orbitals are delocalized over
both defect sites, supporting the “artificial molecules” picture where the
defects have hybridized.

We further understand the implications of these hybridized defect orbitals on
the electronic structure by analyzing the orbital splitting energies and
occupations of the single defects vs. the paired defects in molecular orbital-
like diagrams at $X=1$ in Fig. 3. The energies and occupations of the orbitals
of the single defects CHB and CB within the bandgap are plotted on the left
and right in Fig. 3(a)/(b) and Fig. 3(c), respectively. Their energetic
positions are in qualitative agreement with results in Ref. 11. Upon
hybridization of two single CHB or CB defect orbitals in the double defects,
at $X=1$, we see that the double defect orbitals are symmetrically split by
$\Delta=70$, 290, and 860 meV around the energy of the single defect orbital
for cis $\rm CH_{B}$-CH${\rm{}_{B}}@1$, trans $\rm CH_{B}$-CH${\rm{}_{B}}@1$,
and $\rm C_{B}$-C${\rm{}_{B}}@1$ in the middle of Fig. 3(a), (b), and (c),
respectively. The splitting energy $\Delta$ of the in-plane paired defects
$\rm C_{B}$-C${\rm{}_{B}}@1$ is several times higher than the out-of-plane
defect pairs, cis and trans $\rm CH_{B}$-C${\rm{}_{B}}@1$. The occupation of
these artificial molecular orbitals is based on the total number of electrons
contributed by the individual defects, where the single defect orbital of $\rm
CH_{B}$ or $\rm C_{B}$ is doubly or singly occupied for the charge neutral
configuration in the absence of spin-polarization. The paired defects then
result in full occupation of the artificial molecular orbitals of cis and
trans CHB-CH${\rm{}_{B}}@1$ and occupation of only the bonding orbital of
CB-C${\rm{}_{B}}@1$.

Figure 3: Energy, hybridization, and occupation of single defect orbitals and
hybridized double defect orbitals for (a) cis CHB-CH${}_{\rm B}@1$, (b) trans
CHB-CH${}_{\rm B}@1$, (c) CB-C${}_{\rm B}@1$, and (d) unoptimized CB-C${}_{\rm
B}@1$. The splitting energy is $\Delta$. The bottom and top gray bars
correspond to the valence and conduction bands, respectively. All values shown
are in eV.

The hybridization and occupation of the defect orbitals can be analogized to
the canonical examples of bonding between two hydrogen or helium atoms to form
H2 or He2, where H2 forms a covalent bond, while He2 interacts via weaker
dispersion interactions. In these cases, orbital hybridization similarly
results in a lower energy bonding orbital and a higher energy antibonding
orbital. In the gas phase, only H2 is stable, however, due to the full
occupation of only the lower energy bonding orbital. We can therefore
anticipate the $\rm C_{B}$-C${\rm{}_{B}}@X$ defect pair to become more stable
upon defect orbital hybridization. Furthermore, as the atom-atom distance
increases beyond the equilibrium bond position, we can also expect the
potential energy of the the $\rm C_{B}$-C${\rm{}_{B}}@X$ defect pair to
increase as the orbital overlap decreases.

Figure 4: (a) Wavefunction overlaps, (b) splittings $\Delta$ of bonding-
antibonding orbitals, and (c) relative energies of formation
$E^{r}_{A_{\alpha}\textnormal{-}B_{\beta}@X}$ for cis CHB-CH${}_{\rm B}@X$,
trans CHB-CH${}_{\rm B}@X$, CB-C${}_{\rm B}@X$, and unoptimized CB-C${}_{\rm
B}@X$.

We test this analogy by calculating important signatures of orbital
hybridization for increasing defect-defect distances in Fig. 4. In Fig. 4(a),
we plot the wavefunction overlaps, calculated as the modulus squared of the
inner product of single defect orbitals placed at $X=0$ and $X=1$ to 4. Based
on the molecular bonding analogy, since we observe larger $\Delta$ for the in-
plane defect pair CB-C${\rm{}_{B}}@1$ than for the out-of-plane defect pairs
cis and trans CHB-CH${\rm{}_{B}}@X$, we should observe larger overlaps for the
in-plane paired defects than the out-of-plane ones, and indeed we do for all
$X$. Based further on the molecular analogy of bonding centrosymmetric atoms,
we should also observe lower overlap with increasing defect-defect distance.
From $X=1$ to 3 where the defect-defect distances increases from 2.62 to 5.25
Å, as expected, the overlap decreases. However, at $X=4$, we see for all
defect pairs that the overlap does not continue to decrease, but rather
increases. This trend reversal can be attributed to the non-centrosymmetric
nature of the defect orbitals. Evidently the $X=4$ position is more conducive
than $X=3$ for overlap between the lobes of the three-fold symmetric defect
orbitals, despite the larger defect-defect distance. The trends in splitting
energies $\Delta$ indicating the degree or orbital hybridization generally
mirror those of the overlaps: the in-plane defect pairs have larger overlaps
and thus larger $\Delta$ than the out-of-plane ones, and $\Delta$ decreases
with increasing defect-defect distances until $X=3$ to $X=4$, where both the
overlap and splitting are higher for the latter. The only notable exception is
cis CHB-CH${\rm{}_{B}}@X$, where the overlap increases but splitting decreases
from from $X=1$ to $X=2$. Recalling that overlap is calculated for single
defects at positions 0 and $X=1$ to 4, steric effects may prevent effective
overlap in this paired defect, resulting in lower $\Delta$ at $X=1$.

In Fig. 4(c), we plot dissociation curves for the paired defects, computing
the relative energy of formation $E_{A_{\alpha}\textnormal{-}B_{\beta}@X}^{\rm
r}$ relative to the energy $E_{A_{\alpha}\textnormal{-}B_{\beta}@X}$ of the
geometry corresponding to the farthest defect-defect distance at $X=4$. We
also take into account systematic energy differences for single defects placed
at positions $X=1$ to 4. The overall expression for
$E_{A_{\alpha}\textnormal{-}B_{\beta}@X}^{\rm r}$ is:

$\displaystyle E_{A_{\alpha}\textnormal{-}B_{\beta}@X}^{\rm r}$
$\displaystyle=(E_{A_{\alpha}\textnormal{-}B_{\beta}@X}-E_{A_{\alpha}\textnormal{-}B_{\beta}@4})$
$\displaystyle+(E_{B_{\beta}@X}-E_{B_{\beta}@4}).$ (1)

Based on the molecular analogy, for the $\rm CH_{B}$-CH${\rm{}_{B}}@X$ defects
with fully occupied bonding and antibonding orbitals symmetrically split in
energy around the energy of the single defect orbitals, barring any other
effects, we would expect low-magnitude energies of formation
$E^{r}_{A_{\alpha}\textnormal{-}B_{\beta}@X}$. Indeed, the energy of formation
of cis $\rm CH_{B}$-CH${\rm{}_{B}}@X$ is barely impacted by the defect-defect
interactions. However, its trans counterpart in Fig. 4(b) experiences a
substantial energy barrier of $\sim$ 140 meV for the nearest possible
interaction at $X=1$. Beyond $X=1$, the defect-defect interaction quickly
drops off. The large formation energy of trans $\rm CH_{B}$-$\rm H_{B}@1$ may
therefore be due to a locally strained environment, demonstrating the
importance of the host lattice in determining energetics. These trends for the
out-of-plane defects are in sharp contrast to that of $\rm
C_{B}$-C${\rm{}_{B}}@X$ shown in Fig. 4(c), where the energy of formation is
opposite in sign and three times larger in magnitude at $X=1$ compared to the
energy of formation of trans $\rm CH_{B}$-CH${\rm{}_{B}}@X$, suggesting that
$\rm C_{B}$ defects are more stable when paired closely together because the
lower energy bonding orbital is occupied while the higher energy antibonding
orbital is left unoccupied. We also note that sign of
$E^{r}_{A_{\alpha}\textnormal{-}B_{\beta}@X}$ flips from $X=2$ to $X=3$ before
returning to 0 at $X=4$ by definition, whereas in a traditional molecular
dissociation curve, the potential simply approaches 0 with no sign flip. We
can attribute this feature to the same trend reversal in overlap and $\Delta$
discussed previously, overall illustrating that overlap, $\Delta$, and
$E^{r}_{A_{\alpha}\textnormal{-}B_{\beta}@X}$ depend not only on the defect-
defect distance but also the geometry. We note that we also investigate the
impact of spin-polarization on the open-shell CB single defect and determine
that the SCF energy is uniformly lowered by $\sim 270$ meV for CB@X where
$X=1$ to 4, resulting in no change to the relative formation energy of the
double defects $E_{A_{\alpha}-B_{\beta}@X}^{r}$.

While the analogy of defects hybridizing as atoms do to form molecules can be
convenient, the host lattice can non-trivially impact the electronic structure
of the defects. To determine the quantitative effect of host lattice
relaxation, we show the orbital energies, occupations, and relative formation
energies for the paired defect CB-CB@$X$ and the single defect $\rm C_{B}$ in
Fig. 3(d) and Fig. 4(d) where the CB defects are inserted into a relaxed,
pristine nanoflake of hBN and not geometry optimized again afterward. Similar
to its optimized counterpart, the single defect orbital is singly occupied,
but its energy relative to the valence band is shifted downward by 180 meV.
The orbitals of the double defect are also split generally evenly around this
energy, as in the three, optimized model double defects. Notably, however, the
splitting is in general smaller than in the relaxed CB-CB@$X$ double defects,
suggesting that host lattice relaxation enables more efficient overlap of the
single defect orbitals. The relative decrease in splitting then causes the
magnitude of the relative formation energy to be smaller for close defect-
defect distances at $X=1$ and $X=2$. These results highlight the importance of
including both the defect complex and host lattice in computational studies
and not simply computing the hybridization of the orbitals of single defects
to predict properties of defect complexes.

Figure 5: Double defect and emitter C${}_{\textrm{B}}$-V${}_{\mathrm{N}}@Y$.
(a) Geometry of the double defect C${}_{\textrm{B}}$-V${}_{\mathrm{N}}@5$, or
C${}_{\textrm{B}}$VN, in the same sized nanoflake of hBN as shown in Fig. 1.
For all double defects C${}_{\textrm{B}}$-V${}_{\mathrm{N}}@Y$,
C${}_{\textrm{B}}$ is at position 0 and V${}_{\textrm{N}}$ is at position $Y$
from 5 to 8. Boron is gray, nitrogen is blue, and carbon is yellow. (b)
Molecular orbital diagram of C${}_{\textrm{B}}$-V${}_{\mathrm{N}}@5$ and its
single defect constituents C${}_{\textrm{B}}$ and VN, where all localized
Kohn-Sham orbitals are labelled, as well as the lowest-lying level in the
“conduction” band $\phi_{\mathrm{CB}}$. (c) Splittings $\Delta_{12}$ and
$\Delta_{23}$ between the defect orbitals of
C${}_{\textrm{B}}$-V${}_{\mathrm{N}}@Y$. The lower, dashed and upper, solid
black lines correspond to energy differences between
$\psi^{\mathrm{V_{N}}}_{1}\leftrightarrow\psi^{\mathrm{C_{B}}}$ and
$\psi^{\mathrm{C_{B}}}\leftrightarrow\psi^{\mathrm{V_{N}}}_{2}$, respectively,
in the limit of zero hybridization. (d) Excitation energy of the lowest-lying
excitation in C${}_{\textrm{B}}$-V${}_{\mathrm{N}}@Y$ for $Y=5$ to 8.

We leverage our general understanding of the picture of bonding artificial
atoms into artificial molecules to design new quantum emitters with the
chemical degree of freedom. As an example, we study derivatives of the neutral
C${}_{\textrm{B}}$VN defect. This well-studied defect in hBN is hypothesized
to be a bright emitter of visible light [48, 49] with potential applications
as a single photon emitter in technologies. C${}_{\textrm{B}}$VN is
conventionally treated as a single defect, but here, we consider it a defect
complex C${}_{\textrm{B}}$-V${}_{\mathrm{N}}@Y$ of single defects
C${}_{\textrm{B}}$ and VN, where C${}_{\textrm{B}}$ is at $X=0$ and VN is
placed at $Y=5$ to 8, as labelled in Fig. 5(a). These geometries correspond to
in-plane distances ranging from 1.46 Å to 5.46 Å. In our notation, therefore,
C${}_{\textrm{B}}$VN is equivalent to C${}_{\textrm{B}}$-V${}_{\mathrm{N}}@5$.
As shown in Fig. 5(b), C${}_{\textrm{B}}$-V${}_{\mathrm{N}}@Y$ has three
defect orbitals within the band gap, labelled $\phi_{1}$, $\phi_{2}$, and
$\phi_{3}$ in order of increasing energy. One is from C${}_{\textrm{B}}$
($\psi^{\mathrm{C_{B}}}$) and two are from V${}_{\mathrm{N}}@Y$
($\psi^{\mathrm{V_{N}}}_{1}$ and $\psi^{\mathrm{V_{N}}}_{2}$). We also label
the lowest-lying conduction orbital $\phi_{\mathrm{CB}}$. Because
$\phi_{\mathrm{CB}}$ and $\phi_{3}$ lie so closely in energy (a few meV
apart), numerical error results in some mixing in these two orbitals between
the localized character of a defect orbital and the delocalized character of a
bulk-like band state. The single defect constituents C${}_{\textrm{B}}$ and VN
each contribute one electron to the double defect
C${}_{\textrm{B}}$-V${}_{\mathrm{N}}@Y$, resulting in full occupation of the
lowest energy defect orbital. As with the previously discussed in- and out-of-
plane defect pairs, we find that adjusting the relative geometry from $Y=5$ to
8 of the single defects CB and VN results in energy splitting $\Delta_{12}$
between the localized defect orbitals of CB-V${}_{\mathrm{N}}@Y$, up to 1.48
eV for $Y=5$ as plotted in Fig. 5(c) and closely matching the results of Ref.
48. Notably, however, $\Delta_{23}$ is largely unchanged by varying $Y$, as
$\psi^{\mathrm{C_{B}}}$ and $\psi_{2}^{\mathrm{V_{N}}}$ are well-separated in
energy, resulting in less splitting for a given perturbation. Finally, as
expected, as the defects move farther apart and hybridize less, $\Delta_{12}$
and $\Delta_{23}$ approach the energy differences between
$\psi^{\mathrm{V_{N}}}_{1}\leftrightarrow\psi^{\mathrm{C_{B}}}$ and
$\psi^{\mathrm{C_{B}}}\leftrightarrow\psi^{\mathrm{V_{N}}}_{2}$, respectively,
in the limit of zero hybridization as the lower, dashed and upper, solid black
lines in Fig. 5(c), which correspond to the energy differences between the
single defect orbitals in Fig. 5(b).

To demonstrate the impact of bonding on the optical activity of defects, in
Fig. 5(d), we plot the excitation energies for the lowest-lying excitation of
C${}_{\textrm{B}}$-V${}_{\textrm{N}}@Y$ for $Y=5$ to 8 calculated within the
linear-response framework (Casida method) of the time-dependent density
functional theory [50, 51, 52, 53] and discussed further in Appendix A. From
$Y=5$ to $Y=8$, the energy of this transition follows the same trend in
splitting $\Delta_{12}$, decreasing from 1.89 eV to 1.19 eV, or from red of
the visible spectrum to the near-infrared. These results indicate that
angstrom-level control of defect geometry may be a powerful method for tuning
the excitation energies of defect emitters.

##  III Conclusions and Outlook

The presented results support the picture of artificial molecules formed by
bonding artificial atoms comprised of defects with localized orbitals within
the band gap of a low-dimensional semiconductor or insulator. We theoretically
describe specific defect pairs with varying defect-defect distances in
monolayer flakes of hexagonal boron nitride that remain either in-plane or
move out-of-plane. We show that orbital overlap-induced orbital splitting
ranges from the order of tens of millielectronvolts to nearly 1 eV, realized
via a pair of in-plane defects, $\rm C_{B}$-$\rm C_{B}$, at the closest
defect-defect distance $X=1$ corresponding to $2.62$ Å. In this particular
case, due to full occupation of the lower energy bonding orbital and empty
occupation of the higher energy antibonding orbital, the bonded defect pair is
more energetically favorable than when the defects are separated farther apart
and interact weakly by several hundreds of millielectronvolts. In essence, we
observe an effective chemical bond between defect orbitals. These concepts can
be used to design novel quantum emitters with custom optical properties. As a
example, we show that the absorption wavelength of
C${}_{\textrm{B}}$-V${}_{\mathrm{N}}@Y$ can be tuned from the visible to the
near-infrared range at different relative defect-defect geometries $Y$.

The results of this study drive two major conclusions. The first is that this
study further illustrates how the formation of artificial molecules, or multi-
defect complexes, could be responsible for unassigned zero phonon lines in the
emission of defect centers in hBN [54], and the observation that dramatic
shifts in the emission line of a given defect center cannot be explained with
local strain alone [55]. The second conclusion is that this chemical degree of
freedom, coupled with recent experimental innovations that enable
deterministic placement of defect centers, offers a novel approach for
precisely engineering the optical properties of solid-state quantum emitters,
properties that cannot be readily accessed via coupling to cavities for
example.

We anticipate that our predictions will motivate experimental studies towards
systematic synthesis and characterization of artificial molecules in solid-
state materials via techniques with atomic precision. Promising directions for
further theoretical studies include inclusion of electron-phonon coupling [56,
57, 58], charge state dynamics and charge transition levels [59, 60, 28],
Jahn-Teller effects [31, 27, 61, 62] and application of group theoretical
methods [49], especially to describe the properties of optically excited
states. We expect these results to spur further fundamental and technological
advances in the use of quantum defects-by-design in solid-state materials for
applications in quantum technologies.

##  Appendix A Computational methods

To theoretically model these systems, we use a pseudopotential, real-space
density functional theory (DFT) code Octopus [43, 44, 45]. Optimized molecular
geometries and ground state electronic densities are calculated with SG15
optimized norm-conserving Vanderbilt pseudopotentials [63, 64] and the Perdew,
Burke, and Ernzerhof (PBE) generalized gradient approximation exchange-
correlation functional [65], which has been used in previous studies on the
electronic and optical properties of defects in hBN [11, 48, 66, 67, 68]. Note
that while the PBE functional is known to underestimate the bandgap and
systematically mischaracterize certain optical properties [49], results from
PBE often qualitatively match those from the more accurate and computationally
expensive Heyd–Scuseria–Ernzerhof functional (HSE) functional [69]. The real-
space simulation box is parameterized with a mesh spacing of 0.20 Å, as in
Ref. 68, and consists of spheres with 4 Å radius around each atom. The ground
state self-consistent field (SCF) energies are converged below 1 meV/atom, and
the energy difference $E_{\mathrm{bg}}$ between the highest occupied and
lowest unoccupied bulk-like Kohn-Sham orbitals is converged within 1% of its
mean value of 4.28 eV for nanoflakes with either a single or a pair of
defects. $E_{\mathrm{bg}}$ agrees with the range of bandgaps calculated for
nanoflakes of hBN with approximately 100 boron and nitrogen atoms with either
nitrogen or boron edges [70], in contrast to a calculated bandgap of 4.50 eV
with periodic DFT codes and the PBE functional [68].

Figure 6: (a) Absorption spectra and (b) ionization energies of the double
defect C${}_{\textrm{B}}$-V${}_{\mathrm{N}}@Y$ from $Y=5$ to 8.

The energies of the excited states of interest converged with 100 extra states
in the Casida method, a standard linear-response time-dependent density-
functional theory (TDDFT) [50]. The resulting absorption spectrum of the
lowest-lying excitation is plotted in Fig. 6(a), where the energies correspond
to those plotted in Fig. 5(d) and the absorption is proportional to the
squared transition dipole moment times the transition frequency [51, 52, 53].
We confirm that the adiabatic LDA functional is a valid approximation within
the Casida formalism by checking that the ionization energies plotted in Fig.
6(b) and determined with a $\Delta$SCF calculation are several eV higher than
the excitation energies in the absorption spectra. Notably, we find that the
relative absorption intensities for varying $Y$ depend on the level of linear-
response time-dependent density functional theory applied. Specifically, the
Petersilka [71] and Tamm-Dancoff [72] approximations to the full Casida method
result in different absorption orderings, a well-known issue [73]. Resolving
these discrepancies motivates the development of a group theoretical approach
to predict the optical activity of artificial molecules whose symmetry
strongly depends on the host lattice, defect-defect interactions, and in-
plane/out-of-plane movements. Despite these differences in absorption
intensities, all three methods predict similar transition energy orderings.
Structures and Kohn-Sham orbitals are visualized with XCrySDen [74].

## Acknowledgements

D.S.W. and C.J.C contributed equally to this work. We acknowledge fruitful
discussions with Tomáš Neuman, John Philbin, and Chitraleema Chakraborty. This
work was supported by the Department of Energy ‘Photonics at Thermodynamic
Limits’ Energy Frontier Research Center under grant DE-SC0019140 (theoretical
and computational approaches in light-matter interactions), by the U.S.
Department of Energy, Office of Science, Basic Energy Sciences (BES),
Materials Sciences and Engineering Division under FWP ERKCK47 ‘Understanding
and Controlling Entangled and Correlated Quantum States in Confined Solid-
state Systems Created via Atomic Scale Manipulation’ (spin qubits), and by the
the Army Research Office MURI (Ab Initio Solid-State Quantum Materials) grant
number W911NF-18-1-0431 (computational methods for low dimensional materials).
D.S.W. is an NSF Graduate Research Fellow. P.N. is a Moore Inventor Fellow
through Grant GBMF8048 from the Gordon and Betty Moore Foundation. The
Flatiron Institute is a division of the Simons Foundation. This research used
resources of the National Energy Research Scientific Computing Center, a DOE
Office of Science User Facility supported by the Office of Science of the U.S.
Department of Energy under Contract No.DE-AC02-05CH11231. Additional
calculations were performed using resources from the Department of Defense
High Performance Computing Modernization program as well as resources at the
Research Computing Group at Harvard University.

## References

  * Wrachtrup  _et al._ [2001] J. Wrachtrup, S. Ya Kilin, and A. P. Nizovtsev, “Quantum computation using the 13C nuclear spins nearby the single NV defect center in diamond,” Opt. Spectrosc. 91, 429–437 (2001).
  * Weber  _et al._ [2010] J R Weber, W F Koehl, J B Varley, A Janotti, B B Buckley, C G Van de Walle, and D D Awschalom, “Quantum computing with defects.” Proc. Natl. Acad. Sci. U.S.A 107, 8513–8 (2010).
  * Wrachtrup [2010] Joerg Wrachtrup, “Defect center room-temperature quantum processors,” Proc. Natl. Acad. Sci. U.S.A. 107, 9479–9480 (2010).
  * Aharonovich  _et al._ [2016] Igor Aharonovich, Dirk Englund, and Milos Toth, “Solid-state single-photon emitters,” Nat. Photonics 10, 631–641 (2016).
  * Degen  _et al._ [2017] C. L. Degen, F. Reinhard, and P. Cappellaro, “Quantum sensing,” Rev. Mod. Phys. 89, 035002 (2017).
  * Wang  _et al._ [2020a] Derek S. Wang, Tomá š Neuman, and Prineha Narang, “Dipole-coupled emitters as deterministic entangled photon-pair sources,” Phys. Rev. Research 2, 043328 (2020a).
  * Kurtsiefer  _et al._ [2000] Christian Kurtsiefer, Sonja Mayer, Patrick Zarda, and Harald Weinfurter, “Stable solid-state source of single photons,” Phys. Rev. Lett. 85, 290–293 (2000).
  * Ye  _et al._ [2019] Meng Ye, Hosung Seo, and Giulia Galli, “Spin coherence in two-dimensional materials,” npj Comput. Mater. 5, 1–6 (2019).
  * Atatüre  _et al._ [2018] Mete Atatüre, Dirk Englund, Nick Vamivakas, Sang Yun Lee, and Joerg Wrachtrup, “Material platforms for spin-based photonic quantum technologies,” Nat. Rev. Mater. 3, 38–51 (2018).
  * Childress  _et al._ [2014] Lilian Childress, Ronald Walsworth, and Mikhail Lukin, “Atom-like crystal defects: From quantum computers to biological sensors,” Phys. Today 67, 38–43 (2014).
  * McDougall  _et al._ [2017] Nicholas L. McDougall, Jim G. Partridge, Rebecca J. Nicholls, Salvy P. Russo, and Dougal G. McCulloch, “Influence of point defects on the near edge structure of hexagonal boron nitride,” Phys. Rev. B 96, 144106 (2017).
  * MacKoit-Sinkevičiene  _et al._ [2019] M. MacKoit-Sinkevičiene, 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).
  * Czelej and Majewski [2020] Kamil Czelej and Jacek A Majewski, “Electronic structure and magneto-optical properties of silicon-nitrogen-vacancy complexes in diamond,” Phys. Rev. B 102, 115102 (2020).
  * Rogers  _et al._ [2008] L. J. Rogers, S. Armstrong, M. J. Sellars, and N. B. Manson, “Infrared emission of the NV centre in diamond: Zeeman and uniaxial stress studies,” New. J. Phys. 10, 103024 (2008).
  * Momenzadeh  _et al._ [2015] S. Ali Momenzadeh, Rainer J. Stöhr, Felipe Favaro De Oliveira, Andreas Brunner, Andrej Denisenko, Sen Yang, Friedemann Reinhard, and Jörg Wrachtrup, “Nanoengineered diamond waveguide as a robust bright platform for nanomagnetometry using shallow nitrogen vacancy centers,” Nano Lett. 15, 165–169 (2015).
  * Faraon  _et al._ [2012] Andrei Faraon, Charles Santori, Zhihong Huang, Victor M. Acosta, and Raymond G. Beausoleil, “Coupling of nitrogen-vacancy centers to photonic crystal cavities in monocrystalline diamond,” Phys. Rev. Lett. 109, 033604 (2012).
  * Chakraborty  _et al._ [2019] Chitraleema Chakraborty, Nicholas R. Jungwirth, Gregory D. Fuchs, and A. Nick Vamivakas, “Electrical manipulation of the fine-structure splitting of WSe2 quantum emitters,” Phys. Rev. B 99, 045308 (2019).
  * Zhang  _et al._ [2018] Jingyuan Linda Zhang, Shuo Sun, Michael J. Burek, Constantin Dory, Yan Kai Tzeng, Kevin A. Fischer, Yousif Kelaita, Konstantinos G. Lagoudakis, Marina Radulaski, Zhi Xun Shen, Nicholas A. Melosh, Steven Chu, Marko Lončar, and Jelena Vučković, “Strongly Cavity-Enhanced Spontaneous Emission from Silicon-Vacancy Centers in Diamond,” Nano Lett. 18, 1360–1365 (2018).
  * Machielse  _et al._ [2019] B. Machielse, S. Bogdanovic, S. Meesala, S. Gauthier, M. J. Burek, G. Joe, M. Chalupnik, Y. I. Sohn, J. Holzgrafe, R. E. Evans, C. Chia, H. Atikian, M. K. Bhaskar, D. D. Sukachev, L. Shao, S. Maity, M. D. Lukin, and M. Lončar, “Quantum Interference of Electromechanically Stabilized Emitters in Nanophotonic Devices,” Phys. Rev. X 9, 031022 (2019).
  * Neuman  _et al._ [2020] Tomáš Neuman, Derek S. Wang, and Prineha Narang, “Nanomagnonic cavities for strong spin-magnon coupling,” Phys. Rev. Lett. 125, 247702 (2020).
  * Wang  _et al._ [2020b] Derek S. Wang, Tomáš Neuman, and Prineha Narang, “Spin emitters beyond the point dipole approximation in nanomagnonic cavities,” arXiv:2012.04662 (2020b).
  * Trivedi  _et al._ [2020] Rahul Trivedi, Kevin A. Fischer, Jelena Vučković, and Kai Müller, “Generation of Non‐Classical Light Using Semiconductor Quantum Dots,” Adv. Quantum Technol. 3, 1900007 (2020).
  * Dai  _et al._ [2020] David D. Dai, Derek S. Wang, and Prineha Narang, “Passive controlled-variable phase gate on photonic qubits via cascade emitter,” arXiv:2011.09302 (2020).
  * Narang  _et al._ [2019] Prineha Narang, Christopher J. Ciccarino, Johannes Flick, and Dirk Englund, “Quantum Materials with Atomic Precision: Artificial Atoms in Solids: Ab Initio Design, Control, and Integration of Single Photon Emitters in Artificial Quantum Materials,” Adv. Funct. Mater. 29, 1–7 (2019).
  * Awschalom  _et al._ [2018] David D. Awschalom, Ronald Hanson, Jörg Wrachtrup, and Brian B. Zhou, “Quantum technologies with optically interfaced solid-state spins,” Nat. Photonics 12, 516–527 (2018).
  * Alkauskas  _et al._ [2019] Audrius Alkauskas, Lee C. Bassett, Annemarie L. Exarhos, and Kai-Mei C. Fu, “Defects by design: Quantum nanophotonics in emerging materials,” Nanophotonics 8, 1863–1865 (2019).
  * Thiering and Gali [2018] Gergő Thiering and Adam Gali, “Ab Initio Magneto-Optical Spectrum of Group-IV Vacancy Color Centers in Diamond,” Phys. Rev. X 8, 021063 (2018).
  * Harris  _et al._ [2020] Isaac Harris, Christopher J. Ciccarino, Johannes Flick, Dirk R. Englund, and Prineha Narang, “Group-iii quantum defects in diamond are stable spin-1 color centers,” Phys. Rev. B 102, 195206 (2020).
  * Trusheim  _et al._ [2019] Matthew E. Trusheim, Noel H. Wan, Kevin C. Chen, Christopher J. Ciccarino, Johannes Flick, Ravishankar Sundararaman, Girish Malladi, Eric Bersin, Michael Walsh, Benjamin Lienhard, Hassaram Bakhru, Prineha Narang, and Dirk Englund, “Lead-related quantum emitters in diamond,” Phys. Rev. B 99, 075430 (2019).
  * Ferrenti  _et al._ [2020] Austin M. Ferrenti, Nathalie P. de Leon, Jeff D. Thompson, and Robert J. Cava, “Identifying candidate hosts for quantum defects via data mining,” npj Comput. Mater. 6, 2–7 (2020).
  * Ciccarino  _et al._ [2020] Christopher J. Ciccarino, Johannes Flick, Isaac B. Harris, Matthew E. Trusheim, Dirk R. Englund, and Prineha Narang, “Strong spin–orbit quenching via the product Jahn–Teller effect in neutral group IV qubits in diamond,” npj Quantum Mater. 5, 1–6 (2020).
  * Dyck  _et al._ [2018a] Ondrej Dyck, Songkil Kim, Sergei V. Kalinin, and Stephen Jesse, “E-beam manipulation of Si atoms on graphene edges with an aberration-corrected scanning transmission electron microscope,” Nano Res. 11, 6217–6226 (2018a).
  * Dyck  _et al._ [2018b] Ondrej Dyck, Songkil Kim, Elisa Jimenez-Izal, Anastassia N. Alexandrova, Sergei V. Kalinin, and Stephen Jesse, “Building Structures Atom by Atom via Electron Beam Manipulation,” Small 14, 1–9 (2018b).
  * Fuechsle  _et al._ [2012] Martin Fuechsle, Jill A. Miwa, Suddhasatta Mahapatra, Hoon Ryu, Sunhee Lee, Oliver Warschkow, Lloyd C.L. Hollenberg, Gerhard Klimeck, and Michelle Y. Simmons, “A single-atom transistor,” Nat. Nanotechnol. 7, 242–246 (2012).
  * Kalinin and Dyck [2019] Sergei V Kalinin and Ondrej Dyck, “Materials and Devices with Probes and Beams: Down to the Atomic Level and Back Up,” Adv. Funct. Mater. 29, 33–35 (2019).
  * Ziatdinov  _et al._ [2019] Maxim Ziatdinov, Ondrej Dyck, Xin Li, Bobby G Sumpter, Stephen Jesse, Rama K Vasudevan, and Sergei V Kalinin, “Building and exploring libraries of atomic defects in graphene: Scanning transmission electron and scanning tunneling microscopy study,” Sci. Adv. 5, 1–10 (2019).
  * Dyck  _et al._ [2018c] Ondrej Dyck, Songkil Kim, Sergei V Kalinin, and Stephen Jesse, “Placing single atoms in graphene with a scanning transmission electron microscope microscope,” Appl. Phys. Lett. 111, 113104 (2018c).
  * Susi  _et al._ [2019] Toma Susi, Jannik C Meyer, and Jani Kotakoski, “Quantifying transmission electron microscopy irradiation effects using two-dimensional materials,” Nat. Rev. Phys. 1, 397–405 (2019).
  * Schuler  _et al._ [2019a] Bruno Schuler, Diana Y. Qiu, Sivan Refaely-Abramson, Christoph Kastl, Christopher T. Chen, Sara Barja, Roland J. Koch, D. Frank Ogletree, Shaul Aloni, Adam M. Schwartzberg, Jeffrey B. Neaton, Steven G. Louie, and Alexander Weber-Bargioni, “Large spin-orbit splitting of deep in-gap defect states of engineered sulfur vacancies in monolayer $\text{WS}_{2}$,” Phys. Rev. Lett. 123, 076801 (2019a).
  * Schuler  _et al._ [2019b] Bruno Schuler, Jun-Ho Lee, Christoph Kastl, Katherine A. Cochrane, Christopher T. Chen, Sivan Refaely-Abramson, Shengjun Yuan, Edo van Veen, Rafael Roldán, Nicholas J. Borys, Roland J. Koch, Shaul Aloni, Adam M. Schwartzberg, D. Frank Ogletree, Jeffrey B. Neaton, and Alexander Weber-Bargioni, “How Substitutional Point Defects in Two-Dimensional WS2 Induce Charge Localization, Spin–Orbit Splitting, and Strain,” ACS Nano 13, 10520–10534 (2019b).
  * Tian  _et al._ [2020] Xuezeng Tian, Dennis S. Kim, Shize Yang, Christopher J. Ciccarino, Yongji Gong, Yongsoo Yang, Yao Yang, Blake Duschatko, Yakun Yuan, Pulickel M. Ajayan, Juan Carlos Idrobo, Prineha Narang, and Jianwei Miao, “Correlating the three-dimensional atomic defects and electronic properties of two-dimensional transition metal dichalcogenides,” Nat. Mater. 19, 867–873 (2020).
  * Scherpelz and Galli [2017] Peter Scherpelz and Giulia Galli, “Optimizing surface defects for atomic-scale electronics: Si dangling bonds,” Phys. Rev. Mater. 1, 021602 (2017).
  * Marques  _et al._ [2003] M. A. L. Marques, A. Castro, G. F. Bertsch, and A. Rubio, “octopus: a first-principles tool for excited electron-ion dynamics,” Comput. Phys. Commun. 151, 60–78 (2003).
  * Andrade  _et al._ [2015] X. Andrade, Strubbe. D. A., U. De Giovannini, A. H. Larsen, M. J. T. Oliveira, J. Alberdi-Rodriguez, A. Varas, I. Theophilou, N. Helbig, M. J. Verstraete, L. Stella, F. Nogueira, A. Aspuru-Guzik, A. Castro, M. A. L. Marques, and A. Rubio, “Real-space grids and the Octopus code as tools for the development of new simulation approaches for electronic systems,” Phys. Chem. Chem. Phys. 17, 31371 (2015).
  * Tancogne-Dejean  _et al._ [2020] Nicolas Tancogne-Dejean, Micael J. T. Oliveira, Xavier Andrade, Heiko Appel, Carlos H. Borca, Guillaume Le Breton, Florian Buchholz, Alberto Castro, Stefano Corni, Alfredo A. Correa, Umberto De Giovannini, Alain Delgado, Florian G. Eich, Johannes Flick, Gabriel Gil, Adrián Gomez, Nicole Helbig, Hannes Hübener, René Jestädt, Joaquim Jornet-Somoza, Ask H. Larsen, Irina V. Lebedeva, Martin Lüders, Miguel A. L. Marques, Sebastian T. Ohlmann, Silvio Pipolo, Markus Rampp, Carlo A. Rozzi, David A. Strubbe, Shunsuke A. Sato, Christian Schäfer, Iris Theophilou, Alicia Welden, and Angel Rubio, “Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems,” J. Chem. Phys. 152, 124119 (2020).
  * Barcza  _et al._ [2020] Gergely Barcza, Viktor Ivády, Tibor Szilvási, Márton Vörös, Libor Veis, Adám Gali, and Ors Legeza, “DMRG on top of plane-wave Kohn-Sham orbitals: case study of defected boron nitride,” arXiv:2006.04557 (2020).
  * Reimers  _et al._ [2020] Jeffrey R. Reimers, Jun Shen, Mehran Kianinia, Carlo Bradac, Igor Aharonovich, Michael J. Ford, and Piotr Piecuch, “Photoluminescence, photophysics, and photochemistry of the ${{\mathrm{V}}_{\mathrm{B}}}^{-}$ defect in hexagonal boron nitride,” Phys. Rev. B 102, 144105 (2020).
  * Tawfik  _et al._ [2017] Sherif Abdulkader Tawfik, Sajid Ali, Marco Fronzi, Mehran Kianinia, Toan Trong Tran, Catherine Stampfl, Igor Aharonovich, Milos Toth, and Michael J. Ford, “First-principles investigation of quantum emission from hBN defects,” Nanoscale 9, 13575–13582 (2017).
  * Abdi  _et al._ [2018] Mehdi Abdi, Jyh Pin Chou, Adam Gali, and Martin B. Plenio, “Color Centers in Hexagonal Boron Nitride Monolayers: A Group Theory and Ab Initio Analysis,” ACS Photonics 5, 1967–1976 (2018).
  * Casida [1995] Mark E Casida, “Time-dependent density functional response theory for molecules,” in  _Recent Advances In Density Functional Methods: (Part I)_ (World Scientific, 1995) pp. 155–192.
  * Flick  _et al._ [2019] Johannes Flick, Davis M. Welakuh, Michael Ruggenthaler, Heiko Appel, and Angel Rubio, “Light–matter response in nonrelativistic quantum electrodynamics,” ACS Photonics 6, 2757–2778 (2019).
  * Flick and Narang [2020] Johannes Flick and Prineha Narang, “Ab initio polaritonic potential-energy surfaces for excited-state nanophotonics and polaritonic chemistry,” J. Chem. Phys. 153, 094116 (2020).
  * Wang  _et al._ [2020c] Derek S. Wang, Tomáš Neuman, Johannes Flick, and Prineha Narang, “Weak-to-Strong Light-Matter Coupling and Dissipative Dynamics from First Principles,” arXiv:2002.10461 (2020c).
  * Bommer and Becher [2019] Alexander Bommer and Christoph Becher, “New insights into nonclassical light emission from defects in multi-layer hexagonal boron nitride,” Nanophotonics 8, 2041 (2019).
  * Hayee  _et al._ [2020] Fariah Hayee, Leo Yu, Jingyuan Linda Zhang, Christopher J. Ciccarino, Minh Nguyen, Ann F. Marshall, Igor Aharonovich, Jelena Vučković, Prineha Narang, Tony F. Heinz, and Jennifer A. Dionne, “Revealing multiple classes of stable quantum emitters in hexagonal boron nitride with correlated optical and electron microscopy,” Nat. Mater. 19, 534 (2020).
  * Norambuena  _et al._ [2016] Ariel Norambuena, Sebastián A. Reyes, José Mejía-Lopéz, Adam Gali, and Jerónimo R. Maze, “Microscopic modeling of the effect of phonons on the optical properties of solid-state emitters,” Phys. Rev. B 94, 134305 (2016).
  * Alkauskas  _et al._ [2014] Audrius Alkauskas, Bob B. Buckley, David D. Awschalom, and Chris G. Van de Walle, “First-principles theory of the luminescence lineshape for the triplet transition in diamond NV centres,” New J. Phys. 16, 073026 (2014).
  * Jara  _et al._ [2020] Cesar Jara, T. Rauch, Silvana Botti, Miguel A. L. Marques, A. Norambuena, R. Coto, J. R. Maze, and F. Munoz, “First-principles identification of single photon emitters based on carbon clusters in hexagonal boron nitride,” arXiv:2007.15990 (2020).
  * Freysoldt  _et al._ [2014] Christoph Freysoldt, Blazej Grabowski, Tilmann Hickel, Jörg Neugebauer, Georg Kresse, Anderson Janotti, and Chris G. Van de Walle, “First-principles calculations for point defects in solids,” Rev. Mod. Phys. 86, 253–305 (2014).
  * Dreyer  _et al._ [2018] Cyrus E. Dreyer, Audrius Alkauskas, John L. Lyons, Anderson Janotti, and Chris G. Van de Walle, “First-Principles Calculations of Point Defects for Quantum Technologies,” Annu. Rev. Mater. Res. 48, 1–26 (2018).
  * Thiering and Gali [2019] Gergő Thiering and Adam Gali, “The $(e_{g}\otimes e_{u})\otimes E_{g}$ product Jahn–Teller effect in the neutral group-IV vacancy quantum bits in diamond,” npj Comput. Mater. 5, 1–6 (2019).
  * Ivády  _et al._ [2020] Viktor Ivády, Gergely Barcza, Gergő Thiering, Song Li, Hanen Hamdi, Jyh-Pin Chou, Örs Legeza, and Adam Gali, “Ab initio theory of the negatively charged boron vacancy qubit in hexagonal boron nitride,” npj Comput. Mater. 6, 1–6 (2020).
  * Schlipf and Gygi [2015] Martin Schlipf and François Gygi, “Optimization algorithm for the generation of ONCV pseudopotentials,” Comput. Phys. Commun. 196, 36–44 (2015).
  * Hamann [2013] D. R. Hamann, “Optimized norm-conserving Vanderbilt pseudopotentials,” Phys. Rev. B 88, 085117 (2013).
  * Perdew  _et al._ [1996] John P. Perdew, Kieron Burke, and Matthias Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett. 77, 3865–3868 (1996).
  * Tran  _et al._ [2016] Toan Trong Tran, Kerem Bray, Michael J. Ford, Milos Toth, and Igor Aharonovich, “Quantum emission from hexagonal boron nitride monolayers,” Nat. Nanotechnol. 11, 37–41 (2016).
  * Wu  _et al._ [2017] Feng Wu, Andrew Galatas, Ravishankar Sundararaman, Dario Rocca, and Yuan Ping, “First-principles engineering of charged defects for two-dimensional quantum technologies,” Phys. Rev. Mater. 1, 071001 (2017).
  * Tancogne-Dejean and Rubio [2018] Nicolas Tancogne-Dejean and Angel Rubio, “Atomic-like high-harmonic generation from two-dimensional materials,” Sci. Adv. 4, 1–6 (2018).
  * Ernzerhof [2003] Matthias Ernzerhof, “Hybrid functionals based on a screened coulomb potential,” J. Chem. Phys. 118, 8207 (2003).
  * Maruyama and Okada [2018] Mina Maruyama and Susumu Okada, “Energetics and Electronic Structure of Triangular Hexagonal Boron Nitride Nanoflakes,” Sci. Rep. 8, 1–10 (2018).
  * Petersilka  _et al._ [1996] M. Petersilka, U. J. Gossmann, and E. K.U. Gross, “Excitation energies from time-dependent density-functional theory,” Phys. Rev. Lett. 76, 1212–1215 (1996).
  * Hirata and Head-Gordon [1999] So Hirata and Martin Head-Gordon, “Time-dependent density functional theory within the Tamm-Dancoff approximation,” Chem. Phys. Lett. 314, 291–299 (1999).
  * Hu  _et al._ [2014] Chunping Hu, Osamu Sugino, and Kazuyuki Watanabe, “Performance of Tamm-Dancoff approximation on nonadiabatic couplings by time-dependent density functional theory,” J. Chem. Phys. 140, 054106 (2014).
  * Kokalj [1999] Anton Kokalj, “XCrySDen-a new program for displaying crystalline structures and electron densities,” J. Mol. Graph. Model. 17, 176–179 (1999).


# Calibrating Bayesian UNet++ for Sub-seasonal Forecasting

Büşra Asan Abdullah Akgül Department of Mathematics and Computer Science,
University of Southern Denmark  Alper Ünal Eurasia Institute of Earth
Sciences, Istanbul Technical University  Melih Kandemir Department of
Mathematics and Computer Science, University of Southern Denmark  Gözde Ünal
Department of Computer Engineering, Istanbul Technical University

###### Abstract

Seasonal forecasting is a crucial task when it comes to detecting the extreme
heat and colds that occur due to climate change. Confidence in the predictions
should be reliable since a small increase in the temperatures in a year has a
big impact on the world. Calibration of the neural networks provides a way to
ensure our confidence in the predictions. However, calibrating regression
models is an under-researched topic, especially in forecasters. We calibrate a
UNet++ based architecture, which was shown to outperform physics-based models
in temperature anomalies. We show that with a slight trade-off between
prediction error and calibration error, it is possible to get more reliable
and sharper forecasts. We believe that calibration should be an important part
of safety-critical machine learning applications such as weather forecasters.

##  1 Introduction

Seasonal forecasting is a crucial task when it comes to foreseeing the effects
of climate change, especially in making predictions and decisions based on
these effects. Generating accurate seasonal and sub-seasonal forecasts demands
substantial resources, such as the curation of Coupled Model Intercomparison
Projects (CMIP) datasets. These datasets combine outputs from over a hundred
climate models worldwide, facilitating top-tier climate simulations.
Leveraging the vast data reservoirs from CMIP6 (Eyring et al., 2016), the
latest phase of CMIP, there are ongoing efforts to harness deep learning
methodologies for enhanced climate forecasting. For instance, Luo et al.
(2022) use Bayesian Neural Networks (BNN) with CMIP6 for climate prediction in
the North Atlantic, and Anochi et al. (2021) use CMIP6 to assess
precipitation. Andersson et al. (2021) forecast the change in Arctic sea ice
area with the same dataset. In this work, we also utilize the CMIP6 dataset to
produce well-calibrated and sharp forecasts which are crucial for climate
sciences (Gneiting et al., 2007).

We expand the capabilities of the forecast model introduced in Unal et al.
(2023) using the calibration approach proposed by Kuleshov et al. (2018). This
model is shown to achieve better performance than physics-based methods on
sub-seasonal forecasting, especially at predicting temperature anomalies that
indicate extreme hot and cold temperatures that are crucial for climate
change.

In this work, we calibrate the Bayesian version of the forecaster in a
regression setting. We show that the BNNs produce the most calibrated and
sharp forecasts. We compare the performance of BNNs, Monte Carlo Dropout (MC-
Dropout) Gal and Ghahramani (2016), and Deep Ensemble Lakshminarayanan et al.
(2017) methods for assessing climate forecast uncertainty and sharpness,
exploring the potential for improved reliability through calibration. Our
contributions can be listed as follows:

  1. 1.

We apply calibration to a sub-seasonal forecaster that is able to predict
extreme events better than simulations. We show that calibrating deep learning
models should be a crucial step while applying deep learning to climate
sciences.

  2. 2.

We show that well-calibrated forecasters not only produce better confidence
intervals but may also improve the sharpness of the forecasts.

  3. 3.

This method may be generalized to any other application in climate sciences
that gives critical importance to the reliability of the results such as
extreme events, precipitation, and natural disasters such as earthquakes,
floods, and drought.

##  2 Methodology

We formulate the problem as predicting the monthly average air temperature at
2 meters above the earth surface ($2m$) for each coordinate in a 2D
temperature grid which we will name as the temperature map. Our aim is to
construct a reliable confidence interval for each coordinate since Bayesian
methods often produce uncalibrated results Kuleshov et al. (2018).

We first train the model on CMIP6 climate simulations, then fine-tune it with
ERA5 reanalysis data based on real climate measurements. We denote the 2D
temperature map at time $t$ as $x_{t}$. Train set
$D=\\{X_{t},Y_{t}\\}^{T}_{t=1}$ consists of stacked monthly time-series
temperature maps as the input. The input $X_{t}$ refers to $x_{t-1:t-k-m}$
which denotes the range of the stacked months and $Y_{t}$ corresponds to
$x_{t}$. The periodical month selection process from the given range is
described in Unal et al. (2023) and the same setting is used to make a fair
comparison.

###  2.1 Bayesian UNet++ for Temperature Prediction

For temperature forecasting, we convert UNet++ into a BNN (Goan and Fookes,
2020). BNNs are highly regarded for their capability in quantifying
uncertainty, offering robust insights into predictive models Kristiadi et al.
(2020). Thus, we converted the final three layers of the UNet++ architecture
into Bayesian convolutional layers, where we model the weights $\theta$ of the
neural network as a Gaussian distribution. Letter, we maximize evidence lower
bound Kingma and Welling (2022).

###  2.2 Uncertainty Analysis

Confidence Intervals are used for measuring the uncertainty. Quantiles are
calculated from the predictions and checked whether the correct portion of the
predictions actually conforms to those intervals.

Calibration in neural networks refers that if the confidence interval is
chosen as $95\%$, then the intervals should capture around $95\%$ of the
observed outcomes $Y_{t}$. To measure calibration, we count observations that
stay below the predicted upper bound for the quantile $p$ of the sample $t$,
then normalize with the size of the dataset. A neural network is said to be
calibrated if it satisfies the following

$\frac{1}{T}\sum_{t=1}^{T}(\mathbbm{1}\\{Y_{t}\leq
F_{t}^{-1}(p)\\})\longrightarrow p,$ (1)

as $T\rightarrow\infty$ (Gneiting et al., 2007), where $F_{t}$ refers to the
Cumulative Distribution Function (CDF) of the output of the neural network $H$
for the input $X_{t}$.

Calibrated Regression. We need to match empirical and predicted CDFs according
to Equation 1. Therefore, the training partition of the ERA5 dataset is used
to construct a calibration dataset to map the predicted CDF to the empirical
CDF. We train a new regressor $R:[0,1]\rightarrow[0,1]$ on the calibration
dataset. Thus, we expect $R\circ F_{k}$ to be calibrated. CDFs are
monotonically increasing functions, hence the choice of $R$ is an Isotonic
Regressor (Niculescu-Mizil and Caruana, 2005).

From the training partition $S=\\{X_{t-1:t-k-m},Y_{t}\\}^{T^{\prime}}_{t=1}$
of ERA5, calibration dataset $C=\\{c_{t},y_{t}\\}^{T^{\prime}}_{t=1}$ is
constructed where $c_{t}$ refers to $F_{t}(Y_{t})$ and $y_{t}$ refers to
$\hat{P}(F_{t}(Y_{t}))$ using the dataset generation method in Kuleshov et al.
(2018). $\hat{P}$ is formulated as

$\hat{P}=\frac{1}{T}|\left\\{Y_{t}\left|F_{t}\left(Y_{t}\right)<p,t=1,...,T\right\\}\right|$
(2)

where $|A|$ refers to the cardinality of the set A. It calculates empirical
CDF from the predicted CDF by normalizing the count of output $Y_{t}$ staying
below $p^{th}$ quantile of $F_{t}$.

Table 1: Metrics (lower the better) for calibrated and uncalibrated versions of the models. | Uncalibrated | Calibrated
---|---|---
Metrics | CE | MAE | Sharpness | CE | MAE | Sharpness
UNet++ | N/A | 0.975 | N/A | N/A | N/A | N/A
Bayesian UNet++ | 0.023 | 2.237 | 0.291 |  0.015 $(\downarrow 34.8\%)$ | 2.298 $(\uparrow 2.7\%)$ |  0.274 $(\downarrow 6.9\%)$
Dropout (40%) | 0.131 | 0.993 | 0.853 | 0.035 $(\downarrow 73.2\%)$ |  0.990 $(\downarrow 0.3\%)$ | 0.847 $(\downarrow 0.7\%)$
Deep Ensemble | 0.086 | 1.548 | 0.789 | 0.024 $(\downarrow 70.0\%)$ | 1.366 $(\downarrow 11.8\%)$ | 0.799 $(\uparrow 1.3\%)$

As a result, the estimated $\mathbb{P}(Y\leq F_{X}^{-1}(p))$ by the regressor
$R$ provides the calibrated probability that a random $Y$ falls into the
credible interval so that we can adjust the predicted probability to the
empirical probability.

###  2.3 Training & Evaluation

Figure 1: $50\%$ confidence interval (Top) and $90\%$ confidence interval
(Bottom) of the Bayesian UNet++ for a sample in the North West Coast of
America are given. The mean coverage percentages for confidence intervals are
$63\%$ and $91\%$ for the calibrated, and $66\%$ and $98\%$ for the
uncalibrated models. Thus, we choose a representative sample. Uncalibrated
confidence interval plots are shown on the left (green), and calibrated plots
are on the right (blue). Grey dots refer to the average temperature values for
each month in the given time period (2016-2021). The percentage of the values
falling within the intervals aligns more closely with the expected confidence
levels, both at $50\%$ and $90\%$ in the calibrated model’s plot.

We train our model using the $2m$ temperature variable from 9 ensembles of the
CMIP6 dataset. 1700 samples are separated for training and 100 for validation.
400 samples from the ERA5 dataset are used for fine-tuning and the
construction of the calibration dataset. 116 samples from ERA5 are used in the
evaluation of all methods as in Unal et al. (2023).

Metrics. To measure accuracy, Mean Absolute Error (MAE) is used. MAE for
calibrated models in Table 1 is recalculated using the mid-quantile values
from the calibrated forecaster.

Sharpness (Gneiting et al., 2007) is a metric which is widely used in climate
forecasting. It measures the concentration of the forecasts as

$sharpness(F_{1},F_{2},...,F_{T})=\frac{1}{T}\sum_{t=1}^{T}var(F_{t}).$ (3)

Calibration error (CE) proposed by Kuleshov et al. (2018) is used for
assessing the quality of the calibration of the forecasts as

$\text{CE}(F_{1},Y_{1},...,F_{T},Y_{T})=\frac{1}{m}\sum_{i=1}^{m}w_{j}(p_{j}-\hat{p_{j}})$
(4)

where $m$ refers to the number of confidence levels $0\leq p_{1}\leq...\leq
p_{j}\leq...\leq p_{m}\leq 1$ and $\hat{p_{j}}$ is empirical frequency. In
this setting, $w_{j}$ is chosen as 1.

##  3 Results

We use the experimental settings of Unal et al. (2023). Table 1 illustrates
the impact of calibration on uncertainty quantification methods. The Bayesian
model demonstrates the highest sharpness and calibration as we expected.
However, there exists a trade-off between MAE and CE, with Bayesian
demonstrating the lowest CE, followed by Ensemble and Dropout. Apart from the
reduction in CE, we observe a decrease in MAE for Dropout and Deep Ensemble
which suggests that calibration not only improves the accuracy of the network
performance but also enhances the capture probability percentages of
confidence intervals around point estimates. MAE is calculated using the
actual $50\%$ quantile values predicted by the Isotonic Regressor.

Figure 1 demonstrates that for the calibrated case, roughly $90\%$ of the
$90\%$ confidence intervals capture the true temperature values in the test
dataset. We also observe the same result for $50\%$ interval. Thus, the
proposed calibration produced results in line with the expected proportion of
confidence intervals capturing the true outcome at the given confidence level,
suggesting that the model is well-calibrated.

Figure 2: Calibration plot suggested by Kuleshov et al. (2018) given for a
sample in the grid in Figure 1 to evaluate the calibration of the forecasts.
Each predicted confidence level is plotted against its corresponding expected
confidence level. Predictions illustrate the frequency of observing an outcome
$Y_{t}$ at each level. We expect calibrated models to be closer to $y=x$.

CE is visualized in Figure 2. Equation 4 is applied to the values calculated
for the calibration plot for each quantile, and the mean is used as the CE of
that sample. After the calibration, our model converges to the $y=x$ line
which indicates that the predicted confidences for the samples are closer to
expected confidences, especially for quantiles larger than $50\%$.

##  4 Conclusion

We proposed a method to enhance the sharpness and reliability of weather
forecasts by calibrating them using a CDF-based calibration approach. This
involved transforming the final layers of UNet++ to Bayesian. Periodically
stacked multi-dimensional time-series data used as input. As we designed the
output of the network to produce a CDF, we trained an isotonic regressor to
calibrate the confidence intervals. We benchmarked the calibrated and
uncalibrated results of three uncertainty quantification methods. Furthermore,
we show that calibrating Dropout and Deep Ensemble might increase the accuracy
of the network along with improving the uncertainty quantification.

This work emphasizes the significance of calibrating neural networks while
suggesting potential improvements for forecast reliability. Various fields in
climate sciences can benefit from calibration since uncertainties arise from
incomplete modeling of the earth and the inherent complexity of climate
systems. While our focus was on temperature forecasting, this approach can be
extended to predicting other essential climate variables such as
precipitation, pressure, and wind components.

## References

  * Andersson et al. (2021) Tom R Andersson, J Scott Hosking, María Pérez-Ortiz, Brooks Paige, Andrew Elliott, Chris Russell, Stephen Law, Daniel C Jones, Jeremy Wilkinson, Tony Phillips, et al.  Seasonal arctic sea ice forecasting with probabilistic deep learning.  _Nature Communications_ , 12(1):5124, 2021.
  * Anochi et al. (2021) Juliana Aparecida Anochi, Vinícius Albuquerque de Almeida, and Haroldo Fraga de Campos Velho.  Machine learning for climate precipitation prediction modeling over south america.  _Remote Sensing_ , 13(13):2468, 2021.
  * Eyring et al. (2016) V. Eyring, S. Bony, G. A. Meehl, C. A. Senior, B. Stevens, R. J. Stouffer, and K. E. Taylor.  Overview of the coupled model intercomparison project phase 6 (cmip6) experimental design and organization.  _Geoscientific Model Development_ , 9:1937–1958, 2016.
  * Gal and Ghahramani (2016) Yarin Gal and Zoubin Ghahramani.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In _International Conference on Machine Learning_ , pages 1050–1059. PMLR, 2016.
  * Gneiting et al. (2007) Tilmann Gneiting, Fadoua Balabdaoui, and Adrian E Raftery.  Probabilistic forecasts, calibration and sharpness.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 69(2):243–268, 2007.
  * Goan and Fookes (2020) Ethan Goan and Clinton Fookes.  Bayesian neural networks: An introduction and survey.  In _Case Studies in Applied Bayesian Data Science_ , pages 45–87. Springer International Publishing, 2020.  doi: 10.1007/978-3-030-42553-1˙3.
  * Kingma and Welling (2022) Diederik P Kingma and Max Welling.  Auto-encoding variational bayes, 2022.
  * Kristiadi et al. (2020) Agustinus Kristiadi, Matthias Hein, and Philipp Hennig.  Being bayesian, even just a bit, fixes overconfidence in relu networks, 2020.
  * Kuleshov et al. (2018) Volodymyr Kuleshov, Nathan Fenner, and Stefano Ermon.  Accurate uncertainties for deep learning using calibrated regression.  In _International Conference on Machine Learning_ , pages 2796–2804. PMLR, 2018.
  * Lakshminarayanan et al. (2017) Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell.  Simple and scalable predictive uncertainty estimation using deep ensembles.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Luo et al. (2022) Xihaier Luo, Balasubramanya T Nadiga, Ji Hwan Park, Yihui Ren, Wei Xu, and Shinjae Yoo.  A bayesian deep learning approach to near-term climate prediction.  _Journal of Advances in Modeling Earth Systems_ , 14(10):e2022MS003058, 2022.
  * Niculescu-Mizil and Caruana (2005) Alexandru Niculescu-Mizil and Rich Caruana.  Predicting good probabilities with supervised learning.  In _International Conference on Machine Learning_ , pages 625–632, 2005.
  * Unal et al. (2023) Alper Unal, Busra Asan, Ismail Sezen, Bugra Yesilkaynak, Yusuf Aydin, Mehmet Ilicak, and Gozde Unal.  Climate model-driven seasonal forecasting approach with deep learning.  _Environmental Data Science_ , 2:e29, 2023.  doi: 10.1017/eds.2023.24.