# Towards ab initio identification of paramagnetic substitutional carbon
defects in hexagonal boron nitride acting as quantum bits

Philipp Auburger  Wigner Research Centre for Physics, Hungarian Academy of
Sciences, P.O. Box 49, 1525 Budapest, Hungary  Adam Gali<EMAIL_ADDRESS>Wigner Research Centre for Physics, Hungarian Academy of Sciences, P.O. Box
49, 1525 Budapest, Hungary  Department of Atomic Physics, Budapest University
of Technology and Economics, Budafoki út 8, 1111 Budapest, Hungary

###### Abstract

Paramagnetic substitutional carbon (C${}_{\text{B}}$, C${}_{\text{N}}$)
defects in hexagonal boron nitride (hBN) are discussed as candidates for
quantum bits. Their identification and suitability are approached by means of
photoluminescence (PL), charge transitions, electron paramagnetic resonance,
and optically detected magnetic resonance (ODMR) spectra. Several clear trends
in these are revealed by means of an efficient plane wave periodic supercell
ab initio density functional theory approach. In particular, this yields
insight into the role of the separation between C${}_{\text{B}}$ and
C${}_{\text{N}}$. In most of the cases the charge transition between the
neutral and a singly charged ground state of a defect is predicted to be
experimentally accessible, since the charge transition level (CTL) position
lies within the band gap. A posteriori charge corrections are also discussed.
A near-identification of an experimentally isolated single spin center as the
neutral C${}_{\text{B}}$ point defect was found via comparison of results to
recently observed PL and ODMR spectra.

##  I  Introduction

Single-photon emitters (SPEs) in semiconductors are the most attractive
fundamental building blocks for novel quantum technologies such as quantum
computing [1, 2, 3, 4, 5, 6], quantum information distribution [7, 1, 5],
quantum sensing [1, 5], quantum photonics [7, 8, 9], quantum cryptography [2,
3, 4], hybrid spin-photon interfaces [8] and spin-mechanics interfaces [8].

To create more versatile structures an embedding material with a large band
gap, such as diamond [10], is desirable, because atomic like defect states and
levels isolated from the environment can be controlled. In particular, this
control over the electronic configuration also concerns charge and spin.
Another issue is to engineer the suitable atomic like defects for scalable
quantum architecture where twodimensional (2D) materials offer a natural
platform to this end.

Hence hexagonal boron nitride (hBN), an sp2-bonded strongly covalent [8]
layered van-der-Waals solid with indirect band gap at 5.955 eV [11] and
excellent chemical and thermal stability [2, 3] retaining the wide band gap in
its exfoliated form [8], has recently emerged as a promising alternative [12,
7, 1, 8, 9, 2, 3, 5, 6]. It hosts a plethora of stable optically active
defects with a broad emission range (1.2-5.3 eV [8]) and very high brightness,
also at room temperature [5]. Moreover they exhibit additional favorable
quantum-optical properties including narrow linewidth, high emission into the
zero-phonon line (ZPL) and addressability via spin-selective optical
transitions [7].

Recently, a breakthrough has been achieved by observing single spin resonance
in hBN as spin-dependent fluorescence intensity, i.e., continuous-wave (cw)
optically detected magnetic resonance (ODMR) measurements, where the electron
spin resonance (ESR) spectrum was broadened by the interaction with the
electron spin and the proximate nuclear spins of the host [12]. The origin of
the single spin centers is unknown, in stark contrast to the other ODMR center
in hBN [13] which is the negatively charged boron-vacancy defect [14, 15]. The
negatively charged boron-vacancy was proposed earlier as a quantum bit
candidate by theory [16] which implies a strong predictive power of first
principles calculations. With using the same theoretical approach, this study
now focuses on the recently observed single spin centers with using the
advantage that the observed ESR signals provide direct information about the
spin density distribution of the defect in the paramagnetic ground state [12].

Generally the atomic structure of isolated optical emitters in hBN cannot be
determined directly in experiments. This is due to the vast number of
candidates, e.g. different complexes of substitutional carbon or oxygen atoms
possibly with adjacent vacancies or antisites. However, simple models can be
developed to capture the fundamental motifs before the identification of
specific defects. Substitutional carbon defects are popular candidates, since
they exhibit low formation energies [2] as well as high migration barriers [8]
and carbon interstitials are highly mobile [1]. Indeed Mendelson et al. [7]
recently found a direct correlation between the photoluminescence of visible
SPEs ($E_{\text{ZPL}}$=2.08-2.16 eV upon 532-nm illumination) and implanted
carbon surveying various samples grown via metalorganic vapor-phase epitaxy
(MOVPE), molecular beam epitaxy (MBE) and highly oriented pyrolytic graphite
(HOPG) conversion. Additionally they also associated the spin readout of room-
temperature ODMR performed on MOVPE samples with the density of carbon-related
quantum emitters which was varied from singles to ensembles. Their X-ray
photospectroscopy (XPS) experiments yielded even more evidence for the
presence of carbon point defects and complexes thereof. The neutral
C${}_{\text{B}}$C${}_{\text{N}}$ dimer has been proposed to explain the
ubiquitous 4.1-eV photoluminescence [9]. Since a stepwise change of ZPL
energies has been reported for the single spin centers [12], not only isolated
or adjacent point defects, but also structures with increased separation
between C${}_{\text{B}}$ and C${}_{\text{N}}$ are analyzed. At a later stage
the results might also be helpful for target-oriented engineering of charge
transition level (CTL) positions in quantum technologies. We note that,
parallel to our work, the optical signatures of neutral C${}_{\text{B}}$ and
C${}_{\text{N}}$ pair and trimer configurations have been considered by first
principles studies [1, 8]. In our study, beside the optical signals, we
systematically investigate CTLs associated with the photostability upon
illumination and the ESR signals associated with ODMR spectra too, including
charged C${}_{\text{B}}$ and C${}_{\text{N}}$ pair configurations that are
paramagnetic in their ground state.

This paper is organized as follows. Section II elucidates the theoretical
model. Section III discusses the results for substitutional carbon
(C${}_{\text{B}}$, C${}_{\text{N}}$) defects focusing on paramagnetic
structures. It successively elaborates on photoluminescence (Sec. III.1),
charge transition levels (Sec. III.2) and cw EPR/ODMR spectra (Sec. III.3).
The results are summarized in Sec. III.4. Finally, Sect. IV consolidates the
implications for defect identification and quantum bit engineering, and it
provides an outlook. The Supplemental Material [17] is a collection of data
regarding hyperfine coupling and electronic electric field gradients (EFGs) at
the atomic nuclei.

##  II Methodology

The investigation is conducted within the framework of plane wave periodic
supercell ab initio density functional theory (DFT) using the VASP package
https://www.vasp.at/ [18, 19]. The systems are intralayer defects and
respective defect complexes consisting of one to three substitutional carbon
atoms embedded in bulk hexagonal boron nitride (hBN). The bulk cells are
hexagonal ($8\times 8\times 1$) cells. All calculations are performed for the
$\Gamma$ point.

Screened hybrid DFT, the modified Heyd-Scuzeria-Ernzerhof (HSE) [20, 21] with
dispersion corrections, 32% exact exchange (instead of the original HSE06 at
25%), screening parameter $\omega=0.2/\mathring{\text{A}}$) is used to
describe electronic structures and to calculate total energies. For the given
supercell size it yields an indirect (direct) band gap of 5.950 eV (6.439 eV),
which is almost identical to the experimental value of 5.955 eV [11] (6.42 eV
[22]). Ions are relaxed only in the defective molecular layer as other
deviations from the ideal bulk crystal geometry turned out to be negligible.
It was found that it is more efficient to optimize the geometry of the systems
with semilocal generalized gradient approximation DFT (PBE) [23, 24] including
dispersion corrections. Only the hybrid DFT lattice constants $a=2.487$ Å and
$c=6.459$ Å have to be adjusted afterwards. They are $\leq 2\%$ smaller than
the experimental values [25]. Dispersion corrections (Grimme-D3 [26, 27] with
Becke-Johnson (BJ) damping [28]) are indispensable for an adequate treatment
of the interlayer interaction and hence the determination of the distance
($c/2$) between adjacent layers.

The photoluminescence (PL) spectra are determined after the additionally
required data have been obtained from two separate calculations [29]. The
total energy and geometry of the excited state are obtained from a $\Delta$SCF
calculation [30]. We note that the $\Delta$SCF method inherently contains the
electron-hole interaction which is strong in hBN [31]. For singlets the ZPL
energy is corrected as given in Ref. 9. The ground state phonon modes are
calculated on the aforementioned semilocal gradient DFT level, which is
considered a reliable approximation drastically reducing the computational
effort [14, 1, 8]. Again only the ions in the defective molecular layer are
allowed to move. Furthermore the cw ODMR defect line broadening can also be
accessed, since the hyperfine coupling for the most abundant nuclear spin
active isotopes (11B, 14N, 13C) is analyzed with the aforementioned screened
hybrid DFT [32, 33]. The cw ODMR spectra are calculated based on this data
employing the MATLAB toolbox EasySpin [34, 35, 36]. EasySpin can also take the
nuclear Zeeman and quadrupole interaction into account. The electronic
electric field gradients (EFGs) at the positions of the (quadrupolar) atomic
nuclei ($I\geq 1$) required for the latter are calculated with VASP [37]. The
EFG matrix $V$ is converted into the quadrupole matrix $Q$ according to

$Q=\frac{eq_{n}}{2I(2I-1)h}V\text{,}$ (1)

where $q_{n}$ and $I$ are the nuclear quadrupole moment and nuclear spin
quantum number, respectively [38], and $h$ is the Planck-constant. This leads
to $Q$(MHz)=0.01636$V$(V/Å2) for 11B and $Q$(MHz)=0.02471$V$(V/Å2) for 14N.
Carbon atom has only non-quadrupolar isotopes.

The plane wave energy cutoff is 500 eV throughout. The energies are considered
converged for differences of less than 0.01 meV. In the semilocal generalized
gradient approximation DFT calculations the Hellmann-Feynman force components
within the defective molecular layer are reduced to less than 5
meV/$\mathring{\text{A}}$. The PL spectra were calculated with a Gaussian
smearing of 5 meV.

The total energy of the charged systems, and hence the position of the charge
transition levels (CTLs), is affected by a spurious contribution due to the
interaction between the charge and its periodic images. Hence a posteriori
charge corrections according to Freysoldt et al. [39, 40] are applied for the
charged bulk defects.

##  III  Results

Chejanovsky et al. [12] recorded the PL spectrum of an isolated optical
emitter (cf. Fig. 10) with $E_{\text{ZPL}}$=1.71 eV (cf. Tab. 2) upon
illumination by 633-nm (1.96-eV) laser. Its ground state has to be
paramagnetic with non-integer spin as concluded based on ODMR (c.f. Fig. 12)
and photokinetics. Furthermore, the electronic spin density is essentially
out-of-plane and $\pi$-like as indicated by angular-dependence and magnitude
of the hyperfine coupling. We concluded that paramagnetic substitutional
carbon defects with $S=1/2$ spin states can satisfy these criteria (see Fig. 5
for electronic structure). Since the abundance of $I=1/2$ 13C is low no
hyperfine splitting due to the carbon impurity is expected for the vast
majority of the spin centers.

With up to three adjacent carbon atoms there are six basic ways to construct a
paramagnetic defect (complex) exhibiting paramagnetic $S=1/2$ states: two
neutral isolated point defects (C${}_{\text{B}}$, C${}_{\text{N}}$), two
singly charged dimers (C${}_{\text{B}}$C${}_{\text{N}}$(+),
C${}_{\text{B}}$C${}_{\text{N}}$(-)) and two neutral trimers
(C2C${}_{\text{B}}$, C2C${}_{\text{N}}$) (cf. Figs. 1, 2, and 4,
respectively). Starting from the C${}_{\text{B}}$C${}_{\text{N}}$ dimer
structures with larger separation between the point defects are also
investigated. These are called “donor-acceptor pairs” (DAPs), because for the
neutral case an electron transfer resulting in a non-paramagnetic structure
(C${}_{\text{B}}^{+}$ and C${}_{\text{N}}^{-}$) is energetically favorable
(see the text below as well as Secs. III.1 and III.2). The DAPs are labeled
according to their separation in the ideal lattice with “1” being the dimer.
The five smallest separations are 1, 2, $\sqrt{7}$, $\sqrt{13}$, 4 (cf. Figs.
2, 3).

Figure 1:  Substitutional carbon defects in hBN (part 1, neutral paramagnetic,
from left to right): C${}_{\text{B}}$, C${}_{\text{N}}$; (B in ochre, N in
cyan, C in purple) Figure 2:  Substitutional carbon defects in hBN (part 2,
singly charged paramagnetic, from left to right):
C${}_{\text{B}}$C${}_{\text{N}}$-DAP-1 (dimer),
C${}_{\text{B}}$C${}_{\text{N}}$-DAP-2,
C${}_{\text{B}}$C${}_{\text{N}}$-DAP-$\sqrt{7}$; (B in ochre, N in cyan, C in
purple, cf. text for DAPs) Figure 3:  Substitutional carbon defects in hBN
(part 3, singly charged paramagnetic, from left to right):
C${}_{\text{B}}$C${}_{\text{N}}$-DAP-$\sqrt{13}$,
C${}_{\text{B}}$C${}_{\text{N}}$-DAP-4; (B in ochre, N in cyan, C in purple,
cf. text for DAPs) Figure 4:  Substitutional carbon defects in hBN (part 4,
neutral paramagnetic, from left to right):
C${}_{\text{2}}$C${}_{\text{B}}$-trimer,
C${}_{\text{2}}$C${}_{\text{N}}$-trimer; (B in ochre, N in cyan, C in purple)

We plot the in-gap Kohn-Sham levels of the neutral defects in Fig. 5. The
isolated neutral C${}_{\text{B}}$ and C${}_{\text{N}}$ defects introduce a
single defect level into the fundamental band gap occupied by a single
electron. The occupied level of C${}_{\text{B}}$ lies at higher energy than
the empty level of C${}_{\text{N}}$, thus the electron transfer occurs from
C${}_{\text{B}}$ towards C${}_{\text{N}}$ when the defects are relatively
close to each forming DAP complexes as noted above. This reflects in the
calculated defect levels of DAP defects which show a completely filled low-
energy level (localized around C${}_{\text{N}}$) and an empty high-energy
level (localized around C${}_{\text{B}}$). By combining three carbon defects,
i.e., the considered trimer complexes in our study, the electronic structure
shows the characteristic of the constituent defects. In
C${}_{\text{2}}$C${}_{\text{N}}$-trimer, two low-energy defect levels appear
localized around C${}_{\text{N}}$ atoms and a high-energy level localized
around C${}_{\text{B}}$ atom where the two low-energy defect levels are
occupied by three electrons. In C${}_{\text{2}}$C${}_{\text{B}}$-trimer,
rather one low-energy defect level and two high-energy defect levels appear in
the gap occupied by three electrons.

Figure 5:  Calculated Kohn-Sham defect levels for neutral carbon defects in
the fundamental band gap of hBN. DAP refers to the
C${}_{\text{B}}$C${}_{\text{N}}$-DAP donor acceptor pair defects. The
spinpolarized calculation introduces a gap between the occupied (filled
triangle) and unoccupied (empty triangle) defect levels.

The calculated binding energies are listed in Table 1. The donor-acceptor
pairs are stabilized through charge transfer from C${}_{\text{B}}$ donor to
C${}_{\text{N}}$ acceptor which involves Coulomb-interaction. As the
attractive Coulomb-interaction increases with shorter distances the binding
between C${}_{\text{B}}$ and C${}_{\text{N}}$ is stronger for pairs with
shorter distances. The calculated binding energies clearly show this trend. We
note that this looks to hold for the immediate C-C dimer defect but a C-C bond
is created in that defect, and it cannot be described by the simple donor-
acceptor pair model. Since the isolated carbon impurities introduce defect
wavefunctions localized inside the host sheet and interlayer diffusion of
carriers is hindered in hBN, complexes of carbon substitutional defects within
a single sheet of hBN layer are considered in which we expect direct electron
transfer between the constituting carbon defects. Certainly, different
configurations of DAP and trimer complexes may occur in each hBN layer in
multilayered hBN but we do not expect direct interaction between carbon
defects residing in different layers of hBN.

Table 1: The calculated binding energy of neutral carbon complexes within a layer of hBN. The negative binding energy refers to the favor of complex formation. Defect | Binding energy (eV)
---|---
C${}_{\text{B}}$C${}_{\text{N}}$-DAPs |
1 | -3.93
2 | -2.04
$\sqrt{7}$ | -1.90
$\sqrt{13}$ | -1.52
4 | -1.49
trimers |
C2C${}_{\text{B}}$ | -5.31
C2C${}_{\text{N}}$ | -5.28

###  III.1  Photoluminescence

The ZPL energy $E_{\text{ZPL}}$ is the first criterion for potentially
identifying an experimentally isolated emitter. It has been calculated for all
neutral systems and additionally also for the singly charged DAPs (cf. Tab.
2). The theoretical values range from 0.44 eV to 4.13 eV. $E_{\text{ZPL}}$
depends on the electronic structure of the defect (complex), in particular its
states inside the band gap, and its ionic relaxation triggered via
deexcitation resp. emission.

The first optical excitation of a neutral point defect either transfers the
unpaired electron from the band gap to the conduction band minimum
(C${}_{\text{B}}$, $E_{\text{ZPL}}$=1.695 eV, cf. Fig. 2) or instead transfers
another electron from the valence band maximum into its orbital
(C${}_{\text{N}}$, $E_{\text{ZPL}}$=2.468 eV, cf. Fig. 2). The former electron
transfer is within the majority spin channel, while it is the minority spin
channel for the latter. This behavior was obtained earlier by Jara et al. [1],
but they gave no values for $E_{\text{ZPL}}$.

Although the neutral point defects (C${}_{\text{B}}$, C${}_{\text{N}}$) are
paramagnetic, their interaction leads to non-paramagnetic structures
(C${}_{\text{B}}$C${}_{\text{N}}$-DAPs with “charge structure”
C${}_{\text{B}}^{+}$-C${}_{\text{N}}^{-}$), because transferring the unpaired
electron from C${}_{\text{B}}$ to C${}_{\text{N}}$ is energetically favorable.
The CTLs C${}_{\text{B}}$(+1$\mid$0) and C${}_{\text{N}}$(0$\mid$-1) (cf. Tab.
3) explain this for the limiting case of large separations. All numerical
calculations at finite separations yielded the same behavior. Each
C${}_{\text{B}}$C${}_{\text{N}}$-DAP exhibits two orbitals inside the band
gap. The lower one is fully occupied and localized at C${}_{\text{N}}$, while
the higher one is empty and localized at C${}_{\text{B}}$. The energy
difference between them diminishes with increasing separation corroborating
the results of Linderälv et al. [8].

$E_{\text{ZPL}}$ of the neutral C${}_{\text{B}}$C${}_{\text{N}}$-DAPs
diminishes with separation (cf. Tab. 2) as expected in view of the previous
findings. The neutral C${}_{\text{B}}$C${}_{\text{N}}$ dimer has been proposed
to explain the ubiquitous 4.1-eV photoluminescence [9]. The calculated
$E_{\text{ZPL}}$=4.13 eV is consistent with that and agrees well with the
hybrid DFT (HSE) results of Mackoit-Sinkevičienė et al. (4.31 eV in Ref. 9),
who also employed a bulk model, and Jara et al. (4.12 eV in Ref. 1), who
employed a single layer model instead. On the contrary the underlying gradient
DFT (PBE) is not sufficient to reproduce this result as found by Linderälv et
al. (3.34 eV in Ref. 8).

The “charge structure” of the neutral C${}_{\text{B}}$C${}_{\text{N}}$-DAPs is
C${}_{\text{B}}^{+}$-C${}_{\text{N}}^{-}$. Hence a singly charged DAP can be
created either by removing the excess electron from C${}_{\text{N}}^{-}$
(positive DAP, C${}_{\text{B}}^{+}$-C${}_{\text{N}}^{0}$) or by adding an
electron to C${}_{\text{B}}^{+}$ (negative DAP,
C${}_{\text{B}}^{0}$-C${}_{\text{N}}^{-}$). In both cases the unpaired
electron is localized at one point defect and the excess charge at the other
one.

$E_{\text{ZPL}}$ of both singly charged C${}_{\text{B}}$C${}_{\text{N}}$-DAPs
increases with separation (cf. Tab. 2). It seems to approach a value equal, or
at least similar, to that of the respective isolated neutral (paramagnetic)
point defect. This convergence behavior had to be expected, because it turned
out that the lowest optical excitation is linked to the same electronic
transition as for the isolated neutral point defect. Hence the influence of
the charged (non-paramagnetic) point defect vanishes with separation. The
aforementioned Jara et al. [1] also found that $E_{\text{ZPL}}$ of the singly
charged dimers is too low for the visible range.

Placing a point defect (C${}_{\text{B}}$ or C${}_{\text{N}}$) adjacent to the
dimer (C${}_{\text{B}}$C${}_{\text{N}}$, abbreviated as C2) creates a basic
trimer (C2C${}_{\text{B}}$ or C2C${}_{\text{N}}$ with the desired $S=1/2$
paramagnetic ground state, cf. Fig. 4). $E_{\text{ZPL}}$ of the neutral basic
trimers (1.36 eV for C2C${}_{\text{B}}$, 1.62 eV for C2C${}_{\text{N}}$) is
lower than that of their neutral “building blocks” point defect and dimer (cf.
Tab. 2). This is essentially due to the (overlap of) additional states inside
the band. Jara et al. [1] obtained a similar value for C2C${}_{\text{N}}$
(1.62 eV), but a clearly higher one for C2C${}_{\text{B}}$ (1.65 eV). The
origin of the discrepancy in the latter is unclear.

Beside the characteristic ZPL emission line (cf. Figs. 6, 7, 8, and 9), the
most intense feature in the phonon sideband (PSB) always exhibited an energy
detuning of 0.15-0.17 eV with respect to ZPL peak, except for the neutral
dimer (0.20 eV). The intensity of its replica decreases with the number of
quanta. Linderälv et al. [8] obtained similar results including a mostly
somewhat larger PSB energy detuning of 0.18-0.20 eV for the neutral points
defects and DAPs. According to them, these phonons are slightly distorted
high-frequency in-plane modes of the ideal lattice, which would be at 0.15 eV,
with the only exception (local mode) being the dimer (C-C bond). The present
results are corroborated by Jara et al. [1], who report 0.20 eV resp. 0.16 eV
for neutral C${}_{\text{B}}$C${}_{\text{N}}$ resp. C2C${}_{\text{N}}$ and
ascribe both to a local vibration mode. They also mention values in the same
range for other substitutional carbon defects. Experimental PL spectra with
similar features have been recorded by various groups [12, 7, 41, 6]. E.g.
Mendelson et al. [7] recently observed several emitters with
$E_{\text{ZPL}}$=2.08-2.16 eV and 0.16-0.20 eV PSB energy detuning upon
illumination by 532-nm (2.33-eV) laser. The characteristics of the ODMR spin
readout are in line with those found in a previous report [12] giving
evidences for their connection to carbon-related defects although using a
longer excitation wavelength (at least 594 nm or 2.09 eV) insufficient for
exciting all of these emitters. Hence the singly positive
C${}_{\text{B}}$C${}_{\text{N}}$-DAPs at larger separations (at least
$\sqrt{7}$, cf. Tab. 2) may explain the observations if the defect ground
state spin is $S=1/2$.

Fig. 10 compares the PSBs of optical emitter D1, which was experimentally
isolated by Chejanovsky et al. [12], to those of the candidates for its
identification. D1 has $E_{\text{ZPL}}$=1.71 eV with a very similar PSB to the
calculated ones and is paramagnetic with half-integer spin. We note that the
other reported single photon emitters with longer wavelength ZPLs than that of
1.71-eV D1 center in Ref. 12 have different PL lineshape, so the
characteristic PSBs are distinct from that of D1. On the other hand, all the
considered carbon substitutional defects have very similar PSB features. Thus
our paper now focuses on the identification of D1 center and it is beyond the
scope to identify the other single photon emitters.

An error of 0.2 eV was allowed to take any kind of inaccuracy into account
when the D1 ZPL energy is compared to the calculated ZPL energies. Hence the
present study provides four candidates: C${}_{\text{B}}$(0),
C${}_{\text{B}}$C${}_{\text{N}}$-DAP-2(+),
C${}_{\text{B}}$C${}_{\text{N}}$-DAP-4(-), C2C${}_{\text{N}}$(0) (cf. Tab. 2).
We note that we often anticipate an accuracy of about 0.1 eV for the
calculated ZPLs. Here we expect a larger error for some of the considered
defects as excitation to band edges is involved in the excited state which is
a subject of finite size error for the neutral defects. Tests on
C${}_{\text{B}}$ defect shows that the calculated $E_{\text{ZPL}}$ increases
by about 0.12 eV going from $8\times 8\times 1$ to $11\times 11\times 1$
supercell. Nevertheless, we provide the calculated ZPL energies within meV
throughout the paper which is technically convergent in terms of the given
supercell size ($8\times 8\times 1$) and DFT functional.

Table 2:  Zero-phonon-line energy ($E_{\text{ZPL}}$) of substitutional carbon defects in bulk hBN obtained from screened hybrid DFT (cf. text and Figs. 1-4) and optical emitter D1 (see Ref. 12). For C${}_{\text{B}}$C${}_{\text{N}}$-DAPs we provide ZPL energies at various charge states. Defect |  $E_{\text{ZPL}}$ (eV)
---|---
C${}_{\text{B}}$($0$) | 1.695
C${}_{\text{N}}$($0$) | 2.468
C${}_{\text{B}}$C${}_{\text{N}}$-DAPs | ($+$) | ($0$) | ($-$)
4 | 2.229 | 1.547 | 1.525
$\sqrt{13}$ | 2.194 | 1.655 | 1.446
$\sqrt{7}$ | 2.018 | 1.984 | 1.243
2 | 1.814 | 2.432 | 1.133
1 (dimer) | 1.055 | 4.131 | 0.442
Trimers |  |
C2C${}_{\text{B}}$($0$) | 1.360
C2C${}_{\text{N}}$($0$) | 1.623
Emitter D1 | 1.710
Figure 6:  Calculated photoluminescence (PL) spectrum for the neutral point
defects (C${}_{\text{B}}$ and C${}_{\text{N}}$, cf. Fig. 1) and basic trimers
(C2C${}_{\text{B}}$ and C2C${}_{\text{N}}$, cf. Fig. 4) in bulk hBN Figure 7:
Calculated photoluminescence (PL) spectrum for the five closest neutral DAPs
(cf. Figs. 2 and 3) in bulk hBN Figure 8:  Calculated photoluminescence (PL)
spectrum for the five closest singly positive DAPs (cf. Figs. 2 and 3) in bulk
hBN Figure 9:  Calculated photoluminescence (PL) spectrum for the five closest
singly negative DAPs (cf. Fig. 2 and 3) in bulk hBN Figure 10:  Comparing
calculated phonon sidebands (PSBs) of candidates (cf. Figs. 1-4) for
identification of optical emitter D1 experimentally isolated by Chejanovsky et
al. [12] in hBN

###  III.2  Charge transition levels

To modify the charge state of a specific defect $X$ it has to be understood
how it depends on the Fermi level $E_{\text{F}}$. The transition between two
charge states happens at the CTL position $E_{\text{F}}^{CT}$ where they are
equally favorable. The defect charge $q$ can be modified accordingly if this
transition is between the two most stable states and $E_{\text{F}}^{CT}$ lies
within the band gap.

The CTLs between neutral and singly charged ground states with respect to
valence band maximum ($\epsilon_{\text{VBM}}$) are calculated as

$\begin{gathered}E_{\text{F}(X)}^{\text{CT}(+1\mid
0)}+\epsilon_{\text{VBM}}=E_{0(X)}^{0}-E_{0(X)}^{+1}-\Delta_{+1}(X)\\\
E_{\text{F}(X)}^{\text{CT}(0\mid-1)}+\epsilon_{\text{VBM}}=E_{0(X)}^{-1}-E_{0(X)}^{0}+\Delta_{-1}(X)\end{gathered}$
(2)

where $E_{0(X)}^{q}$ is the total energy of the defect $X$ with charge $q$ and
$\Delta_{q}(X)$ is the charge correction. We note that we believe that the
threedimensional model corresponds to the experimental conditions where the
defects were found in multilayer hBN on a substrate which introduce a charge
screening.

Tab. 3 summarizes the CTL positions. All uncorrected CTLs lie within the band
gap as displayed in Fig. 11. C${}_{\text{N}}$(+1$\mid$0) resp.
C${}_{\text{B}}$(0$\mid$-1) is closest to the VBM resp. CBM. Since the charge
structure of the neutral DAPs is C${}_{\text{B}}^{+}$-C${}_{\text{N}}^{-}$ it
is possible to either add an electron to C${}_{\text{B}}^{+}$ or to remove an
electron from C${}_{\text{N}}^{-}$ in the presence of the other charged point
defect. Hence C${}_{\text{B}}$C${}_{\text{N}}$(+1$\mid$0) resp.
C${}_{\text{B}}$C${}_{\text{N}}$(0$\mid$-1) converges to
C${}_{\text{N}}$(0$\mid$-1) resp. C${}_{\text{B}}$(+1$\mid$0) with increasing
separation. Meanwhile C${}_{\text{B}}$C${}_{\text{N}}$(+1$\mid$0) is raised
whereas C${}_{\text{B}}$C${}_{\text{N}}$(0$\mid$-1) is lowered. This is not
only due to the diminishing, and finally vanishing, stabilizing effect of the
Coulomb interaction between the charged point defects on the neutral ground
state (cf. Eq. (2)). It is also due to the decreasing energy difference
between the orbitals involved in the charge transitions (highest occupied at
C${}_{\text{N}}^{-}$ and lowest empty at C${}_{\text{B}}^{+}$). As a
consequence of this a “forbidden region” in which the Fermi level
$E_{\text{F}}$ can be modified without affecting the charge of the point
defects and DAPs exists between C${}_{\text{N}}$(0$\mid$-1) and
C${}_{\text{B}}$(+1$\mid$0). The CTLs have also been calculated for the basic
trimers C2C${}_{\text{B}}$ and C2C${}_{\text{N}}$.
C2C${}_{\text{B}}$(+1$\mid$0) and C2C${}_{\text{N}}$(0$\mid$-1) are similar to
those of the isolated point defects without the “additional dimer” C2, whereas
C2C${}_{\text{B}}$(0$\mid$-1) and C2C${}_{\text{N}}$(+1$\mid$0) resemble those
of C2.

The charge corrections stabilize the charged systems by several hundreds of
meV (cf. Tab. 3) but the expected trend of convergence is not obvious for the
corrected CTLs (cf. Fig. 11).

Note that similar values for the corrected CTL positions of point defects and
dimer were found by Mackoit-Sinkevičienė et al. [9], while Weston et al. [2,
3] also give a similar value for C${}_{\text{B}}$(+1$\mid$0) but a somewhat
lower one for C${}_{\text{N}}$(0$\mid$-1). Just like in this work both
employed a bulk model, hybrid DFT and charge corrections according to
Freysoldt et al. [39, 40].

Table 3:  Charge transition levels (CTLs) of substitutional carbon defects in bulk hBN relative to the valence band maximum (VBM) obtained from screened hybrid DFT with (and without) a posteriori charge corrections [39, 40] (cf. text and Fig. 11) 222green: outside the indirect band gap (above the conduction band minimum (CBM)), blue: close to the CBM, orange: most likely above the CBM (too large oscillations in potential alignment) Defect | CTL(+1$\mid$0) (eV) | CTL(0$\mid$-1) (eV)
---|---|---
C${}_{\text{B}}$ | 3.81 (4.11) |  6.39 (5.79)
C${}_{\text{N}}$ | 0.21 (0.47) | 3.27 (2.72)
C${}_{\text{B}}$C${}_{\text{N}}$-DAPs |  |
1 | 0.85 (1.13) |  5.99 (5.42)
2 | 1.73 (2.01) | 5.25 (4.68)
$\sqrt{7}$ | 2.02 (2.31) | 5.10 (4.53)
$\sqrt{13}$ | 2.29 (2.57) | 4.81 (4.23)
4 | 2.37 (2.65) | 4.77 (4.20)
Trimers |  |
C2C${}_{\text{B}}$ | 3.66 (3.94) |  unclear (5.40)
C2C${}_{\text{N}}$ | 0.71 (0.98) | 3.29 (2.72)
Figure 11:  Charge transition levels (CTLs) of substitutional carbon defects
in bulk hBN relative to the valence band maximum (VBM) obtained directly from
screened hybrid DFT (light green (point defects) and dark green (complexes))
and with a posteriori charge corrections (light blue (point defects) and dark
blue (complexes)) [39, 40] (cf. text and Tab. 3)

###  III.3  cw ODMR spectra

The cw ODMR spectrum, in particular the line shape resp. broadening, serves as
a further criterion for identifying experimentally isolated emitters. It has
been analyzed for all ground state $S=1/2$ systems (doublets), i.e. the
neutral point defects and basic trimers as well as the singly charged DAPs
(see also [17]).

The theoretical ODMR spectra are mostly obtained based on a simplified spin
Hamiltonian $\hat{H}=\hat{H}_{\text{EZI}}+\sum_{k}\hat{H}_{\text{HFI}}(k)$
taking into account the interaction of the unpaired electron with the external
magnetic field B (electron Zeeman interaction $\hat{H}_{\text{EZI}}$) and the
spins of the nuclei $k$ (hyperfine interaction
$\sum_{k}\hat{H}_{\text{HFI}}(k)$). In the chosen approximation they are
determined by the external magnetic field B (assuming a gyromagnetic tensor
$g=2*\mathbb{1}_{3}$ for the electron), the isotopic composition and the
associated hyperfine matrices $A(k)$ from hybrid DFT. The hyperfine parameters
are calculated for the most abundant nuclear spin active isotopes (11B, 14N,
13C). Since nuclear spin active carbon is rare (13C of 1.1%), it is assumed
that only boron (11B of 80.1%) and nitrogen (14N of 99.6%) contribute to the
hyperfine splitting. For the neutral point defects and singly charged DAPs up
to second neighbors (d=$\sqrt{3}$) of the paramagnetic carbon atom carrying
the unpaired electron are taken into account. For the neutral basic trimers
those of each carbon atom were chosen instead. The results and the
justification of the underlying assumptions are discussed in the following.

We also considered a full spin Hamiltonian for test cases, e.g.,
C${}_{\text{B}}$ defect, which reads as

$\hat{H}=\hat{H}_{\text{EZI}}+\sum_{k}(\hat{H}_{\text{HFI}}+\hat{H}_{\text{NZI}}+\hat{H}_{\text{NQI}})(k)\text{,}$
(3)

which includes $\hat{H}_{\text{NZI}}$ nuclear Zeeman-term and
$\hat{H}_{\text{NQI}}$ nuclear quadrupole term too. The terms can be in Eq.
(3) expressed as

$\begin{gathered}\hat{H}_{\text{EZI}}=\mu_{\text{B}}\mathbf{B}^{\text{T}}g\hat{\mathbf{S}}\\\
\approx\mu_{\text{B}}\begin{pmatrix}B_{x}&B_{y}&B_{z}\end{pmatrix}\begin{pmatrix}2&&\\\
&2&\\\ &&2\end{pmatrix}\begin{pmatrix}\hat{S}_{x}\\\ \hat{S}_{y}\\\
\hat{S}_{z}\end{pmatrix}\\\
=2\mu_{\text{B}}\langle\mathbf{B}.\hat{\mathbf{S}}\rangle\\\
\hat{H}_{\text{HFI}}(k)=h\hat{\mathbf{S}}^{\text{T}}A(k)\hat{\mathbf{I}}_{k}\\\
\hat{H}_{\text{NZI}}(k)=-\mu_{k}g_{k}\langle\mathbf{B}.\hat{\mathbf{I}}_{k}\rangle\\\
\hat{H}_{\text{NQI}}(k)=h\hat{\mathbf{I}}_{k}^{\text{T}}Q(k)\hat{\mathbf{I}}_{k}\text{,}\end{gathered}$
(4)

where $\mu_{\text{B}}$ and $\mu_{k}$ are the Bohr-magneton of the electron and
nuclei $k$, $g_{k}$ and $Q(k)$ are the gyromagnetic factor and quadrupole
moment of nuclei $k$. In terms of broadening of the ESR signal, the most
important parameter is the hyperfine tensor $A(k)$ for nuclei $k$. We will
show that $\hat{H}_{\text{NZI}}$ and $\hat{H}_{\text{NQI}}$ can be typically
neglected at non-zero magnetic fields, although, $\hat{H}_{\text{NQI}}$ can
contribute to the final shape of the ESR or ODMR spectrum coming from the 14N
nuclei spins.

Tab. 4 summarizes the calculated line broadening and peak positions at an
external magnetic field of 42 Gauss as used for the ODMR measurement of D1
spin center (cf. Fig. 12 and Ref. 12). The line broadening
FWHM${}_{\text{(Gauss)}}$=2$\sqrt{2\ln(2)}\sigma_{\nu}$ is defined as the FWHM
of a Gaussian normal distribution with the same standard deviation
$\sigma_{\nu}$. The original peak position due to the electron Zeeman
interaction is at 118 MHz (cf. Tab. 9). The asymmetric influence of the
hyperfine interaction induces a defect-dependent blueshift of up to 15 MHz.
Neutral C${}_{\text{N}}$ exhibits a much larger line broadening (74 MHz) than
neutral C${}_{\text{B}}$ (43 MHz) as expected in view of the hyperfine
constants (cf. Tabs. 5 and 6 and [17]). The underlying reason is that the
gyromagnetic ratio of the boron isotope ($\gamma_{I}/2\pi$ = 1366 Hz/Gauss for
11B), which is a first neighbor of the paramagnetic carbon for
C${}_{\text{N}}$ but only a second neighbor for C${}_{\text{B}}$, is more than
four times larger than that of the nitrogen isotope ($\gamma_{I}/2\pi$ = 308
Hz/Gauss for 14N). In view of the charge structure the line broadening of the
singly positive resp. negative DAP has to converge to that of neutral
C${}_{\text{N}}$ resp. C${}_{\text{B}}$ with increasing separation. While no
unambiguous trend is found for the former, the presence of the charged point
defect increases it for the latter. The line broadening of the basic trimers
is smaller than that of the corresponding point defect without the adjacent
dimer. We note that the C2C${}_{\text{N}}$ defect create a dense eigenvalue
spectrum of the spin Hamiltonian with contribution of 14 nuclear spins which
will result in a relatively smooth ODMR spectrum in the simulation.

Table 4:  Calculated ODMR line broadening FWHM${}_{\text{(Gauss)}}$ (and peak center $\overline{\nu}$) for $S=1/2$ substitutional carbon defects in bulk hBN at an external magnetic field of 42 Gauss: FWHM${}_{\text{(Gauss)}}$=2$\sqrt{2\ln(2)}\sigma_{\nu}$ is the FWHM of a Gaussian normal distribution with the same standard deviation $\sigma_{\nu}$. $\overline{\nu}$ is the expectancy value of the microwave frequency. See also [17]. Defect type | FWHM${}_{\text{(Gauss)}}$ (MHz) ($\overline{\nu}$ (MHz))
---|---
Point defects | C${}_{\text{N}}$(0) | C${}_{\text{B}}$(0)
| 74 (132) | 43 (119)
C${}_{\text{B}}$C${}_{\text{N}}$-DAPs | positive (+) | negative (-)
4 | 74 (132) | 42 (119)
$\sqrt{13}$ | 73 (132) | 44 (119)
$\sqrt{7}$ | 77 (133) | 49 (119)
2 | 71 (131) | 54 (119)
1 (dimer) | 77 (125) | 70 (118)
Trimers | C2C${}_{\text{N}}$(0) | C2C${}_{\text{B}}$(0)
| 51 (127) | 27 (119)

The experimental ODMR spectrum of D1 is compared to those calculated for the
defect identification candidates C${}_{\text{B}}$(0),
C${}_{\text{B}}$C${}_{\text{N}}$-DAP-2(+),
C${}_{\text{B}}$C${}_{\text{N}}$-DAP-4(-), C2C${}_{\text{N}}$(0) (cf. Fig.
12). It is concluded that D1 has to be either C${}_{\text{B}}$(0) or at least
a very similar structure where the spin density is localized around
C${}_{\text{B}}$. As previously discussed the theoretical model fulfills the
conditions found by Chejanovsky et al. [12] and all defect identification
candidates exhibit PL in good agreement with the measured spectrum. However
C${}_{\text{B}}$(0) is the simplest and best explanation of the observed ODMR,
while non-similar structures provide only inferior spectra. The observed ODMR
peak is centered at 112-113 MHz. A constant shift of every single calculated
spectrum has been introduced to match this value, since the magnetic field
measurement (42 Gauss) may have an error of a few Gauss. For
C${}_{\text{B}}$(0) and C${}_{\text{B}}$C${}_{\text{N}}$-DAP-4($-$) ($-7$ MHz)
it is smaller than for C2C${}_{\text{N}}$(0) ($-14$ MHz) and
C${}_{\text{B}}$C${}_{\text{N}}$-DAP-2(+) ($-19$ MHz). While
C${}_{\text{B}}$(0) with 43 MHz line broadening yields a very good
reproduction of the observed peak, the spectra of C2C${}_{\text{N}}$(0) (51
MHz) and especially C${}_{\text{B}}$C${}_{\text{N}}$-DAP-2(+) (71 MHz) are
clearly broader. In view of the charge structure the results for
C${}_{\text{B}}$(0) and C${}_{\text{B}}$C${}_{\text{N}}$-DAP-4($-$) are close
to each other. In both defects, the spin density is localized around
C${}_{\text{B}}$, so one can conclude that C${}_{\text{B}}$ defect should be
involved in the D1 spin center.

Although the experimental background might reach up to 40% of the peak
maximum, it can not lead to another interpretation of the results or a
substantial reduction of their quality. The measured intensities within, and
only within, the interval 92-134 MHz are $>$40% while they come close to this
value at some other frequencies. It is hence an upper bound for the
background. Since even after subtraction of the full 40% (as background
correction) the normalized peak within the relevant interval is only slightly
less broad, so the conclusions remain unchanged.

The hyperfine constants of C${}_{\text{B}}$(0) are given in Tab. 5. Note that
the involvement of nuclear spin active carbon ($\gamma_{I}/2\pi$ = 1071
Hz/Gauss for 13C) would substantially split the ODMR peak of
C${}_{\text{B}}$(0) to a doublet because of the huge spin density on the
carbon atom attributed to the unpaired electron localized on its $p_{z}$
orbital. This behavior is also found for the other $S=1/2$ defects [17].
Measuring it would provide further insight into the structure of D1 and/or
other experimentally isolated emitters.

The hyperfine parameters include the contribution of the spin polarization of
the core electrons to the Fermi contact interaction. Tab. 6 takes neutral
C${}_{\text{N}}$ as an example to show that this contribution can be
significant, in particular if the carbon is nuclear spin active. Without it
Sajid et al. [4] found hyperfine constants in reasonable agreement employing a
2D model and the regular HSE06 functional. Our calculations imply that the
core polarization should be included in the calculations of accurate hyperfine
tensors.

Figure 12:  Measured ODMR spectrum of isolated optical emitter D1 [12] compared to line broadening calculated for the defect identification candidates C${}_{\text{B}}$(0), C${}_{\text{B}}$C${}_{\text{N}}$-DAP-2(+), C${}_{\text{B}}$C${}_{\text{N}}$-DAP-4(-), C2C${}_{\text{N}}$(0) in bulk hBN, all at an external magnetic field of 42 Gauss, theoretical spectra shifted to match the experimental peak position (cf. text) Table 5:  Hyperfine constants ($A_{xx}$, $A_{yy}$, $A_{zz}$) of the neutral C${}_{\text{B}}$ point defect in bulk hBN up to second neighbors obtained from screened hybrid DFT (# $\coloneqq$ numbers of atoms, d $\coloneqq$ distance within ideal lattice, cf. text and Fig. 1) total | # | d(C${}_{\text{B}}$) |  $A_{xx}$ (MHz) |  $A_{yy}$ (MHz) |  $A_{zz}$ (MHz)
---|---|---|---|---|---
C | 1 | 0 | 12.1 | 12.1 | 231.3
N | 3 | 1 | -9.0 | -5.1 | -9.0
B | 6 | $\sqrt{3}$ | 1.4 | -0.9 | 6.1
Table 6:  Hyperfine constants ($A_{xx}$, $A_{yy}$, $A_{zz}$) of the neutral C${}_{\text{N}}$ point defect in bulk hBN up to second neighbors obtained from screened hybrid DFT ($A_{1c}$ $\coloneqq$ core contribution, () $\coloneqq$ value without $A_{1c}$, # $\coloneqq$ numbers of atoms, d $\coloneqq$ distance within ideal lattice, cf. text and Fig. 1) total | # | d(C${}_{\text{N}}$) |  $A_{xx}$ (MHz) |  $A_{yy}$ (MHz) |  $A_{zz}$ (MHz) |  $A_{1c}$ (MHz)
---|---|---|---|---|---|---
C | 1 | 0 | -19.2 (44.8) | -19.2 (44.8) | 156.5 (220.5) | -64.0
B | 3 | 1 | -17.9 (-16.7) | -16.0 (-14.8) | -24.5 (-23.3) | -1.2
N | 6 | $\sqrt{3}$ | -0.2 (0.8) | -0.8 (0.1) | 2.4 (3.4) | -1.0

The influence of the isotopic composition on the ODMR spectra of the point
defects (see Tabs. 7 and 8) was analyzed, because 10B is also significantly
abundant (19.9%). Replacing 11B with 10B gradually reduces the line
broadening, because the hyperfine matrices scale linearly with the nuclear
gyromagnetic factor (g${}_{\text{n}}$ = $0.600$ instead of $1.792$) although
the nuclear spin is increased ($I=3$ instead of $3/2$) [42]. This effect is
more pronounced for neutral C${}_{\text{N}}$, since the boron isotopes are
first neighbors of the paramagnetic carbon. In this case it also affects the
peak center which is slightly redshifted. Considering only 11B is usually a
good approximation, since large changes of the ODMR spectra (dominant
contribution of 10B) are unlikely to be observed at the natural distribution
of isotopes.

Table 7:  Calculated effect of the isotopic composition on ODMR line broadening FWHM (and peak center $\overline{\nu}$) for the neutral C${}_{\text{B}}$ point defect in bulk hBN at an external magnetic field of 42 Gauss (boron isotopes are 2nd neighbors (d = $\sqrt{3}$) of C${}_{\text{B}}$, p $\coloneqq$ abundance of isotopic composition) #(11B) | #(10B) | p(%) |  $\overline{\nu}$ (MHz) | FWHM (MHz)
---|---|---|---|---
6 | 0 | 26.41 | 119 | 43
5 | 1 | 39.37 | 119 | 41
4 | 2 | 24.45 | 119 | 39
3 | 3 | 8.10 | 119 | 37
2 | 4 | 1.51 | 119 | 34
1 | 5 | 0.15 | 119 | 32
0 | 6 | 0.01 | 119 | 29
Table 8:  Calculated effect of the isotopic composition on ODMR line broadening FWHM (and peak center $\overline{\nu}$) for the neutral C${}_{\text{N}}$ point defect in bulk hBN at an external magnetic field of 42 Gauss (boron isotopes are 1st neighbours ($d=1$) of C${}_{\text{N}}$, p $\coloneqq$ abundance of isotopic composition) #(11B) | #(10B) | p (%) |  $\overline{\nu}$ (MHz) | FWHM (MHz)
---|---|---|---|---
3 | 0 | 51.39 | 132 | 74
2 | 1 | 38.30 | 129 | 66
1 | 2 | 9.52 | 126 | 57
0 | 3 | 0.79 | 123 | 45

The peak broadening is essentially determined by the hyperfine interaction and
the simplified spin Hamiltonian (evaluated with second-order perturbation
theory) yields already a good approximation of the ODMR spectra. This is
demonstrated using the most basic example of the neutral point defects
C${}_{\text{B}}$ and C${}_{\text{N}}$. However it has to be kept is mind that
a comprehensive discussion of the full spin Hamiltonian is beyond the scope of
the present study. Tab. 9 summarizes the results for line broadening and peak
center position obtained from different approaches. Further simplifying the
spin Hamiltonian via using only the hyperfine constants, instead of the full
hyperfine tensor $A$ is not justified. It leads to a significant
overestimation of the line broadening, in particular for C${}_{\text{N}}$ (38
MHz) whose peak center is additionally redshifted (5 MHz). Since second-order
perturbation theory cannot be combined with exact diagonalization for specific
nuclei in the code, first-order perturbation theory is also considered. Since
the C${}_{\text{N}}$ peak center is redshifted by 14 MHz, this seems to be a
non-accurate approach. Exact diagonalization allows to include the nuclear
Zeeman interaction and additionally the nuclear quadrupole interaction, which
have an indirect effect on the electronic transitions as the HFI couples the
electronic spin to the nuclear spins. Due to computational expense which
scales exponentially with the number of spins in the Hilbert space it was
performed only for the first (1st) neighbors of the respective point defect.
Including the NZI reverses the 14 MHz redshift of the C${}_{\text{N}}$ peak
center. However this might also be related to the underlying change of the
algorithm to solve the spin Hamiltonian. Note that the NZI leads to core
transition peaks in a low frequency range ($<$30MHz) that was not included in
the calculation of peak center and line broadening. Judging from Fig. 13 the
NQI can still have a non-negligible effect on the spectral shape although peak
center and line broadening are hardly affected (cf. Tab. 9). Indeed the
entries of the $Q$ matrices (cf. Eqs. (1) and (3) and [17]) are $\leq$1MHz.
Hence the line splittings induced by the NQI can only be correspondingly
small. However it can also induce significant intensity redistribution
explaining its importance in the given context. The results for
C${}_{\text{B}}$ are further evidence for its near-identification as
experimentally isolated optical emitter D1.

Table 9:  ODMR line broadening FWHM (and peak center $\overline{\nu}$) for the neutral point defects (C${}_{\text{B}}$, C${}_{\text{N}}$) in bulk hBN at an external magnetic field of 42 Gauss calculated with different approaches (see text for explanation). Labels “perturb2” and “perturb1” are second-order and first-order perturbation theories whereas “hybrid” means an exact diagonalization for NZI and NQI terms but the other terms are considered within second-order perturbation theory. Approach | FWHM (MHz) ($\overline{\nu}$ (MHz))
---|---
| C${}_{\text{N}}$(0) | C${}_{\text{B}}$(0)
EZI only |  0 (118) |  0 (118)
perturb2 ($A$ constants) |  112 (127) | 49 (119)
perturb2 (full $A$) |  74 (132) | 43 (119)
perturb1 (full $A$) |  73 (118) | 43 (118)
hybrid (1st: NZI) |  72 (132) | 42 (119)
hybrid (1st: NZI, NQI) |  72 (132) | 42 (119)
Figure 13:  Influence of the nuclear quadrupole interaction (NQI) on the
theoretical cw EPR/ODMR spectrum of the neutral point defects in bulk hBN at
an external magnetic field of 42 Gauss (a/c) C${}_{\text{B}}$/C${}_{\text{N}}$
without NQI (b/d) C${}_{\text{B}}$/C${}_{\text{N}}$ with NQI for 1st
neighbours (see text for explanation, cf. Tab. 9)

###  III.4  Summary

Paramagnetic substitutional carbon defects in hexagonal boron nitride have
been characterized via modeling their photoluminescence spectra, charge
transition levels and cw ODMR spectra. Clear trends in all of these were found
by analyzing the separation between C${}_{\text{B}}$ and C${}_{\text{N}}$.
Furthermore a near-identification of an experimentally isolated optical
emitter D1 as neutral C${}_{\text{B}}$ is given.

$E_{\text{ZPL}}$=1.695 eV (2.468 eV) was obtained for neutral C${}_{\text{B}}$
(C${}_{\text{N}}$) which are approximate values due to the finite size effects
and the accuracy in the applied DFT functional. $E_{\text{ZPL}}$ of the singly
negative (positive) C${}_{\text{B}}$C${}_{\text{N}}$-DAP converges to this, or
at least a similar, value from below with increasing separation. This behavior
is explained based on the relevant electronic transition and the charge
structure. The result for the neutral dimer (C${}_{\text{B}}$C${}_{\text{N}}$,
$E_{\text{ZPL}}$=4.131 eV), which allegedly causes the ubiquitous 4.1-eV
photoluminescence [9], is consistent with previous hybrid DFT results [9, 1].
$E_{\text{ZPL}}$ of the neutral C${}_{\text{B}}$C${}_{\text{N}}$-DAPs has to
decrease with separation due to the position of the states inside the band gap
(see also [8]). Each basic trimer has a lower $E_{\text{ZPL}}$ than its
“building blocks” point defect and dimer as expected in view of the overlap of
their states inside the band gap.

An intense feature in PSB next to ZPL at 0.15-0.17 eV (except for the neutral
dimer at 0.20 eV) was found for all defects. Similar results were obtained in
theoretical [1, 8] as well as experimental [12, 7, 41, 6] studies. In
particular ODMR experiments yielded evidences for the connection of these
features to carbon-related defects [7, 12]. Four candidates for the
identification of D1 spin center were provided based on the agreement with the
ZPL energies within 0.2 eV accuracy.

The cw ODMR spectra for all ground state $S=1/2$ systems have been calculated
at 42 Gauss external magnetic field based on a simplified spin Hamiltonian.
The asymmetric influence of the hyperfine interaction induces a defect-
dependent blueshift of the peak positions from 118 MHz up to 133 MHz.
C${}_{\text{N}}$ exhibits a much larger line broadening than C${}_{\text{B}}$
due to the hyperfine parameters. Replacing 11B with the less abundant 10B
reduces the peak width and can also shift the peak center. However large
changes are unlikely to be observed at the natural distribution of isotopes.
The charge structure explains the convergence of the line broadening of the
singly charged C${}_{\text{B}}$C${}_{\text{N}}$-DAPs.

The calculated ZPL energy and PSB as well as the ESR spectrum strongly implies
that C${}_{\text{B}}$ is involved as spin-active part of the D1 spin center,
and it could be the isolated neutral C${}_{\text{B}}$ defect. Measuring the
ODMR splitting due to 13C isotope would provide further insight about the
chemical composition of D1 or other emitters.

An exemplary analysis of the point defects demonstrates that the ODMR line
broadening is essentially due to hyperfine interaction. The simplified spin
Hamiltonian together with second-order perturbation theory results in a good
approximation of the spectra. Using only the hyperfine constants, i.e.
neglecting the Euler angles, is not appropriate as it leads to a large
overestimation of the line broadening. A hybrid approach allowed to include
the nuclear Zeeman and additionally the quadrupole interaction for the first
neighbors. The line broadening remained essentially unaffected. Nuclear
quadrupole interaction is important for an exact description of the spectral
shape as it can lead to non-negligible modifications of the intensities but it
does not effect the overall broadening of the ESR signal.

##  IV  Conclusions and outlook

The modeling of PL and ODMR spectra has yielded helpful criteria for the
identification of paramagnetic substitutional carbon defects in hBN. An
example for a near-identification as neutral C${}_{\text{B}}$ has been given.
The separation between C${}_{\text{B}}$ and C${}_{\text{N}}$ emerged as key to
control the ZPL wavelength and charge compensation of defects. All findings
are further steps towards the identification of these structures and later
their target-oriented engineering as quantum bits.

It is demonstrated that first principles theory is able to predict the complex
electron spin resonance spectrum of defects in hBN. It has been found that the
full hyperfine tensor of the nuclei is required to accurately determine the
linewidth and the central position of the ODMR signal at a given external
magnetic fields. As direct measurement of the full hyperfine tensors is
challenging our results demonstrate that the accurate interpretation of ESR
and ODMR signals in hBN requires a tight cooperation between experimental and
ab initio theoretical spectroscopy.

## Acknowledgements

We are grateful to Prof. Jörg Wrachtrup and Dr. Durga Dasari for sharing and
explaining their experimental results, Christoph Freysoldt for discussions
about his charge correction scheme and EasySpin creator Stefan Stoll for
helping us to use his program. This research was supported by the Ministry of
Innovation and Technology and the National Research, Development and
Innovation Office (NKFIH) within the Quantum Information National Laboratory
of Hungary, the National Quantum Technology Program (NKFIH Grant No.
2017-1.2.1-NKP-2017-00001), the National Excellence Program (NKFIH Grant No.
KKP129866), as well as the support of European Commission within the Quantum
Technology Flagship Project ASTERIQS (Grant No. 820394). We acknowledge that
the results of this research have been achieved using the DECI resource Eagle
HPC based in Poland at Poznan with support from the PRACE aisbl and resources
provided by the Hungarian Governmental Information Technology Development
Agency, including the project “gallium”.

## References

  * 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).
  * Weston  _et al._ [2018] L. Weston, D. Wickramaratne, M. Mackoit, A. Alkauskas, and C. G. V. de Walle, Native point defects and impurities in hexagonal boron nitride, Phys. Rev. B 97, 214104 (13) (2018).
  * Weston  _et al._ [2020] L. Weston, D. Wickramaratne, M. Mackoit, A. Alkauskas, and C. G. V. de Walle, Erratum: Native point defects and impurities in hexagonal boron nitride [phys. rev. b 97, 214104 (2018)], Phys. Rev. B 102, 099903 (2) (2020).
  * 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 (9) (2018).
  * 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).
  * 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).
  * 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. Mater. 20, 321 (2021).
  * Linderälv  _et al._ [2021] C. Linderälv, W. Wieczorek, and P. Erhart, Vibrational signatures for the identification of single-photon emitters in hexagonal boron nitride, Phys. Rev. B 103, 115421 (13) (2021).
  * Mackoit-Sinkevičienė  _et al._ [2019] M. Mackoit-Sinkevičienė, M. Maciaszek, C. G. V. 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 (4) (2019).
  * Korneychuk  _et al._ [2018] S. Korneychuk, G. Guzzinati, and J. Verbeeck, Measurement of the indirect band gap of diamond with eels in stem, Phys. Status Solidi A 215, 1800318 (7) (2018).
  * Cassabois  _et al._ [2016] G. Cassabois, P. Valvin, and B. Gil, Hexagonal boron nitride is an indirect bandgap semiconductor, Nat. Photonics 10, 262 (2016).
  * Chejanovsky  _et al._ [2021] N. Chejanovsky, A. Mukherjee, J. Geng, Y.-C. Chen, Y. Kim, A. Denisenko, A. Finkler, T. Taniguchi, K. Watanabe, D. B. R. Dasari, P. Auburger, A. Gali, J. H. Smet, and J. Wrachtrup, Single-spin resonance in a van der waals embedded paramagnetic defect, Nat. Mater. 10.1038/s41563-021-00979-4 (2021).
  * 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).
  * Ivády  _et al._ [2020] V. Ivády, G. Barcza, G. Thiering, S. Li, H. Hamdi, J.-P. Chou, O. Legeza, and A. Gali, Ab initio theory of the negatively charged boron vacancy qubit in hexagonal boron nitride, npj Comp. Mat. 6, 1 (2020).
  * Sajid  _et al._ [2020] A. Sajid, K. S. Thygesen, J. R. Reimers, and M. J. Ford, Edge effects on optically detected magnetic resonance of vacancy defects in hexagonal boron nitride, Communications Physics 3, 1 (2020).
  * Abdi  _et al._ [2018] M. Abdi, J.-P. Chou, A. Gali, and M. B. Plenio, Color Centers in Hexagonal Boron Nitride Monolayers: A Group Theory and Ab Initio Analysis, ACS Photonics 5, 1967 (2018).
  * [17] See the Supplemental Material at … to be provided by the editor … for a collection of data regarding hyperfine coupling and electronic electric field gradients (EFGs) at the atomic nuclei. The data for the full matrices can be made available upon request.
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B. 54, 11169 (18) (1996).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Heyd  _et al._ [2003] J. Heyd, G. E. Scuseria, and M. Ernzerhof, Hybrid functionals based on a screened coulomb potential, J. Chem. Phys. 118, 8207 (2003).
  * Heyd  _et al._ [2006] J. Heyd, G. E. Scuseria, and M. Ernzerhof, Erratum: Hybrid functionals based on a screened coulomb potential [j. chem. phys. 118, 8207 (2003)], J. Chem. Phys. 124, 219906 (1) (2006).
  * Doan  _et al._ [2016] T. C. Doan, J. Li, J. Y. Lin, and H. X. Jiang, Bandgap and exciton binding energies of hexagonal boron nitride probed by photocurrent excitation spectroscopy, Appl. Phys. Lett. 109, 122101 (5) (2016).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Perdew  _et al._ [1997] J. P. Perdew, K. Burke, and M. Ernzerhof, Erratum: Generalized gradient approximation made simple [phys. rev. lett. 77, 3865 (1996)], Phys. Rev. Lett. 78, 1396 (1997).
  * Paszkowicz  _et al._ [2002] W. Paszkowicz, J. B. Pelka, M. Knapp, T. Szyszko, and S. Podsiadlo, Lattice parameters and anisotropic thermal expansion of hexagonal boron nitride in the 10–297.5k temperature range, Appl. Phys. A 75, 431 (2002).
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (dft-d) for the 94 elements h-pu, J. Chem. Phys. 132, 154104 (19) (2010).
  * Smith  _et al._ [2016] D. G. A. Smith, L. A. Burns, K. Patkowski, and C. D. Sherrill, Revised damping parameters for the d3 dispersion correction to density functional theory, J. Phys. Chem. Lett. 7, 2197 (2016).
  * Grimme  _et al._ [2011] S. Grimme, S. Ehrlich, and L. Goerigk, Effect of the damping function in dispersion corrected density functional theory, J. Comput. Chem. 32, 1456 (2011).
  * Alkauskas  _et al._ [2014] A. Alkauskas, B. B. Buckley, D. D. Awschalom, and C. G. V. de Walle, First-principles theory of the luminescence lineshape for the triplet transition in diamond nv centres, New J. Phys. 16, 073026 (23) (2014).
  * Gali  _et al._ [2009] A. Gali, E. Janzén, P. Deák, G. Kresse, and E. Kaxiras, Theory of Spin-Conserving Excitation of the $N\ensuremath{-}{V}^{\ensuremath{-}}$ Center in Diamond, Physical Review Letters 103, 186404 (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).
  * Blöchl [2000] P. E. Blöchl, First-principles calculations of defects in oxygen-deficient silica exposed to hydrogen, Phys. Rev. B 62, 6158 (2000).
  * Szász  _et al._ [2013] K. Szász, T. Hornos, M. Marsman, and A. Gali, Hyperfine coupling of point defects in semiconductors by hybrid density functional calculations: The role of core spin polarization, Phys. Rev. B 88, 075202 (7) (2013).
  * Stoll [2003] S. Stoll,  _Spectral simulations in solid-state electron paramagnetic resonance_ , (phd) dissertation, ETH Zürich (2003).
  * Stoll and Schweiger [2006] S. Stoll and A. Schweiger, Easyspin, a comprehensive software package for spectral simulation and analysis in epr, J. Magn. Reson. 178, 42 (2006).
  * Stoll [2015] S. Stoll, Chapter six - cw-epr spectral simulations: Solid state, in  _Electron Paramagnetic Resonance Investigations of Biological Systems by Using Spin Labels, Spin Probes, and Intrinsic Metal Ions, Part A_, Methods in Enzymology, Vol. 563, edited by P. Z. Qin and K. Warncke (Academic Press, 2015) pp. 121–142.
  * Petrilli  _et al._ [1998] H. M. Petrilli, P. E. Blöchl, P. Blaha, and K. Schwarz, Electric-field-gradient calculations using the projector augmented wave method, Phys. Rev. B 57, 14690 (1998).
  * Stoll and Goldfarb [2017] S. Stoll and D. Goldfarb, Epr interactions – nuclear quadrupole couplings, in  _eMagRes_, Vol. 6 (American Cancer Society, 2017) pp. 495–510, https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470034590.emrstm1504 .
  * Freysoldt  _et al._ [2011] C. Freysoldt, J. Neugebauer, and C. G. V. de Walle, Electrostatic interactions between charged defects in supercells, Phys. Stat. Sol. B 248, 1067 (2011).
  * Freysoldt  _et al._ [2009] C. Freysoldt, J. Neugebauer, and C. G. V. de Walle, Fully Ab Initio finite-size corrections for charged-defect supercell calculations, Phys. Rev. Lett. 102, 016402 (4) (2009).
  * Museur  _et al._ [2008] L. Museur, E. Feldbach, and A. Kanaev, Defect-related photoluminescence of hexagonal boron nitride, Phys. Rev. B 78, 155204 (8) (2008).
  * [42] At this point the order of replacement is irrelevant in view of the nuclear symmetry and the neglect of core-core interaction in the spin Hamiltonian.


results to the individual spectra (“by observation”), and to the source-
representative spectra. See §5 for details. The number of degrees of freedom
for each fit (d.o.f.) is shown in parentheses after each best-fitting
$\chi^{2}$ value.


# Co-optimization of a piezoelectric energy harvesting system for broadband
operation

S Zhao1 U Radhakrishna2 S Hanly3 J Ma4 J H Lang2 and D Buss5 1 School of
Microelectronics, Tianjin University, Tianjin 300072, China  2 Department of
Electrical Engineering & Computer Science, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  3 Mide Technology, Medford, MA 02155,
USA  4 School of Computers, Guangdong University of Technology, Guangzhou
510006, China  5 Texas Instruments, Dallas, Texas 75243, USA
<EMAIL_ADDRESS>
###### Abstract

The goal of this research is to increase the bandwidth (BW) over which
substantial energy can be harvested using a piezoelectric energy harvester
(PEH). The key innovation is the use of bias-flip (BF) electronics at the
output of a PEH having a large electromechanical coupling coefficient
$\kappa_{e}^{2}$. For a PEH with large $\kappa_{e}^{2}$, the open-circuit
resonance frequency $f_{oc}$ is substantially larger than the short-circuit
resonance frequency $f_{sc}$. Over the intervening range, the reactive part of
the conjugate matched load impedance is small, and can be approximated using
BF electronics in which the BF voltage is sufficiently small and the BF losses
are small. This results in a large BW over which substantial energy can be
harvested. Experimental results using a commercially available PEH are
presented to demonstrate this concept. Design guidelines are provided for
achieving PEHs having increased $\kappa_{e}^{2}$.

##  1 Introduction

Kinetic energy harvesters have become a viable source of electric power for
low-power wireless sensor networks that might be connected to the Internet of
Things. Of the many harvesters through which kinetic energy can be converted
to electric energy, PEHs are often favored for their simple structure, their
high power density, their high output voltage, and their ease of self
starting. However, for the common low-loss resonant PEHs, the harvested power
drops off dramatically when the ambient vibration frequency deviates from the
resonance frequency of the PEH; only at the resonant frequency can the PEH
output maximum power. This is problematic when PEH manufacturing tolerances
and/or uncertainties in the vibration spectrum result in a frequency mismatch.
In such cases it is desirable to increase the energy harvesting BW.

Several nonlinear interface circuits such as SSHI, BF, SECE and SSDCI have
been proposed to increase the harvested power [1-5]. These techniques were
demonstrated to be efficient at the resonance frequency of the PEH. At
frequencies away from resonance, it was subsequently shown that BF employing
universal phase can effectively optimize output power [6]. The BF technique
shares features that are common with the independently developed P-SSHI-$\phi$
technique [7]. Even so, for a PEH with a small coupling coefficient
($\kappa_{e}^{2}$), the BF losses limit output power away from resonance
because the BF voltage for frequencies away from resonance becomes large.
Therefore, to compensate for the limited output power, this paper studies the
resonance splitting behaviour of PEHs with large $\kappa_{e}^{2}$, which was
initially explored in [8] with only resistive load. Resonance splitting with
large $\kappa_{e}^{2}$ gives two well-separated output power peaks between
which (and slightly beyond) BF electronics can harvest near peak power,
thereby extending the harvesting BW. With a combination of BF employing
universal phase, and a PEH having large $\kappa_{e}^{2}$, the energy
harvesting BW can be further increased. Experimental results with a
commercially-available PEH having $\kappa_{e}^{2}$=0.069 are presented to
demonstrate this concept. With the help of BF, a peak power of 112 $\mu$W, and
a 3-dB BW of 50 Hz over which substantial energy is harvested, are achieved.
Additionally, this paper provides design guidelines for achieving PEHs having
larger $\kappa_{e}^{2}$, and a new PEH is suggested to improve the BW
performance of the existing commercial PEH.

Figure 1: Schematic of a cantilever PEH.

Table 1: Four4 frequencies that are important to understanding a cantilever PEH. Value (Hz) | Definition
---|---
$f_{sc}$=673.00 | short-circuit resonance frequency
$f_{ZR1}$=673.43 | lower zero Thevenin reactance frequency
$f_{ZR2}$=695.72 | upper zero Thevenin reactance frequency
$f_{oc}$=696.00 | open circuit resonance frequency
|

##  2 Characterization of a Commercial PEH

The energy harvesting approach proposed here is studied experimentally with a
commercial PEH. Mide PPA2014 PEH was selected because of its large
$\kappa_{e}^{2}$. Figure 1 shows a cross section of this cantilever device.
Figure 2a shows the equivalent circuit for the PEH that defines the compact
model (CM) parameters. Extracted CM parameters and output voltage measurements
made using 1-g acceleration and a 5 mm clamping overlap are shown in Figure
2b. Four frequencies are important to understand power output from a
cantilevered PEH. They are explained in Table 1. The matched-load (ML)
resistance $R_{ML}$ gives maximum output power at the lower zero-Thevenin-
reactance frequency $f_{ZR1}$.

At the two frequencies $f_{ZR1}$ and $f_{ZR2}$, the reactive part of the
impedance looking into the output terminals of the circuit model in Figure 2a
is zero. At these frequencies the source force (voltage $V_{F}$) and velocity
(current $I_{S}$) are in phase, and maximum power is delivered to a resistive
load that is “matched” to the internal resistance of the device.

##  3 Bias-Flip Circuit

Figure 2c shows the DC rectification and storage (DCRS) circuit used here. The
BF inductor, shown in the shaded area, is used to implement the tunable
inductive or capacitive reactive impedance that is required to deliver maximum
output power to the storage cell at a given vibration frequency. Operation of
the BF or “switched” inductor has been described in detail elsewhere [1-3],
and used in a variety of ways to increase output power from a PEH [4, 5].
Here, the bias, that is $V_{Out}$ across $C_{P}$, is flipped so as to produce
a waveform having a phase $\phi(V_{Opt})$ in Figure 2d that would result if a
matched reactive load were used in place of the BF inductor.

The BF timing is determined as follows. Begin with the AC Matched Load (ACML)
circuit shown in Figure 2d. At each vibration frequency, the conjugate-matched
impedance that delivers maximum real power to the matched resistive load is
determined. With this impedance the phase $\phi(f)$ of the voltage $V_{Out}$
relative to the source voltage $V_{F}$ is determined to be

$\phi(f)=tan^{-1}\left[({1}/{(2\pi fC_{m})}-{2\pi fL_{m}})/{R_{m}}\right]$ (1)

The bias-flipped voltage $V_{Out}$ is then flipped at the phase $\phi(f)$
relative to $V_{F}$. Despite the BF timing, to achieve maximum output power,
in ACML case the matched resistive load $R_{L}$ is required, while for DCRS
case an optimum rectification voltage $V_{Rect}$ is required.

It is well-known that the maximum output power from a PEH is given by

$P_{Opt}={m^{2}a^{2}}/{(8\eta)}={V_{F}^{2}}/{(8R_{m})}$ (2)

For the ACML circuit, this maximum output power is achieved at all
frequencies. The black dashed curve of Figure 4a shows that, if the conjugate
matched reactive load is replaced by a BF circuit having 100% BF efficiency
(no inductor loss), and the resistive load is replaced by an ideal rectifier
having the proper $V_{Rect}$, the power delivered to the storage cell is
within 5% of $P_{Opt}$ over a wide frequency range. As will be seen in Figure
3, the $V_{Out}$ waveform in the BF/DCRS circuit looks very different from the
sinusoidal voltage in the ACML circuit. In spite of this, ACML output power is
a good predictor for BF/DCRS output power.

Figure 2: (a) Equivalent circuit of a PEH, (b) CM parameters, simulated and
measured voltage curves of the Mide PPA2014 PEH for various resistive loads
with a 5 mm clamping overlap, (c) DCRS circuit and (d) ACML circuit.

A limitation of the BF technique is that its efficiency is typically between
80% and 90% [6]; the losses occur in the inductor and switch. When a PEH with
low $\kappa_{e}^{2}$ is used, the BF voltage for frequencies away from
resonance becomes large, and BF losses limit output power [3, 7]. However, for
PEHs having large $\kappa_{e}^{2}$, $f_{sc}$ and $f_{oc}$ become well
separated, the BF voltage in this range becomes low, and BF losses become
small, thereby increasing the BW over which substantial energy can be
harvested. The low BF voltage in this range is illustrated by the simulations
in Figure 3. This behavior can be understood by analogy to the ACML circuit.
For each vibration frequency $f$, there is a matched load resistor that gives
optimum output power, and this resistor defines a resonance frequency
$f_{res}$. Between BF pulses, the voltage oscillates at $f_{res}$. For the
region between $f_{sc}$ and $f_{oc}$, $f_{res}$ is close to the vibration
frequency $f$. Because of this, a very small tweak is required by the BF
circuit to achieve the correct phase of $V_{Out}$. Far above $f_{oc}$ and far
below $f_{sc}$, there is a large difference between $f_{res}$ and $f$, and a
large BF voltage is required. In summary, the combination of BF and large
$\kappa_{e}^{2}$ can be used to significantly extend the range over which a
PEH can produce (near) maximum power output without the need for high-voltage
loading electronics.

##  4 Experimental and Simulation Results

To demonstrate widening of the frequency range over which (near) maximum power
can be harvested, experiments and simulations are performed with the Mide
PPA2014 PEH described in Section 2. A 1-g acceleration and a 5 mm clamping
overlap are used. The results shown in Figure 4 are based on the DCRS circuit
shown in Figure 2c in which the rectifier is a diode bridge. Output power is
optimized by varying $V_{Rect}$ at each frequency. The blue curve shows that,
even with 82% BF efficiency, the output power is very near the theoretical
maximum (green curve) between $f_{sc}$ and $f_{oc}$, given ideal diodes. Diode
voltage drop degrades output power (red curve) but the resulting 3-dB BW of 50
Hz (7.4% of $f_{sc}$) represents a substantial increase over a PEH with
negligible coupling, where the 3-dB BW is $f_{sc}/Q_{ML}$=17Hz (2.5% of
$f_{sc}$).

Figure 3: Waveforms of $V_{Out}$ for selected vibration frequencies. For $f\ll
f_{sc}$ and $f\gg f_{oc}$, the BF voltage $V_{BF}$ is large, and BF loss is
significant. However, for $f_{sc}<f<f_{oc}$ $V_{BF}$ is small. At the two
zero-Thevenin reactance frequencies $f_{ZR1}\approx f_{sc}$ and at
$f_{ZR2}\approx f_{oc}$, $V_{BF}=0V$. Figure 4: (a) Maximum output power and
(b) corresponding value of $V_{Rect}$ that gives maximum output power for a
variety of load conditions.

Figure 5 shows the effect of varying $\kappa_{e}^{2}$. These simulations show
output power for PEHs having the CM parameters shown in Figure 2b, but with
the parameter A varied to provide different $\kappa_{e}^{2}$, where A
represents the coupling between mechanical and electrical energy in a PEH. The
simulations are made assuming a realistic 82% BF efficiency. They further
assume the use of a “smart” switching rectifier in the DCRS/BF circuit of
Figure 2c to eliminate the diode-voltage loss around $f_{sc}$. When
$\kappa_{e}^{2}$ is small, maximum output power occurs at
$f_{ZR1}$$\approx$$f_{sc}$. As $\kappa_{e}^{2}$ increases, the 3-dB BW
increases. Use of the “smart” switching rectifier increases the BW with the
present device ($\kappa_{e}^{2}$=0.069) to 78 Hz (11.5% of $f_{sc}$); if
$\kappa_{e}^{2}$ can be increased by a factor of 2, the 3-dB BW increases to
104 Hz (15.4% of $f_{sc}$).

##  5 Design for Large $\kappa_{e}^{2}$

The coupling coefficient $\kappa_{e}^{2}$ is given by [9]

$\kappa_{e}^{2}={A^{2}C_{m}}/{C_{P}}={A^{2}}/{(kC_{P})}={A^{2}}/{((k_{PE}+k_{non-
PE})C_{P})}={(f_{oc}^{2}-f_{sc}^{2})}/{f_{sc}^{2}}$ (3)

where $C_{m}$ represents the beam stiffness, $C_{P}$ represents the plate
capacitance of the piezoelectric material, $k$ is the total spring constant of
a PEH, $k_{non-PE}$ and $k_{PE}$ are spring constants of non-PE and PE
component respectively. For the bimorph cantilever of Figure 1, A is given by
[10]

$A={3WY_{PE}b_{PE}d_{31}}/{L}$ (4)

where W and L are width and length of a PEH, $b_{PE}$ is the distance from the
center-line of the cantilever to the center-line of the piezoelectric (PE)
layer, $Y_{PE}$ is the Young’s modulus of PE material, and $d_{31}$ is the
piezoelectric coefficient.

Figure 5: Output power for PEHs having different $\kappa_{e}^{2}$.

Figure 6: (a) Cross section of baseline device and (b) proposed redesign.

A and $C_{P}$ depend only on the PE material. It is apparent that $k_{non-PE}$
should be made as small as possible. Thus, a composite cantilever designed for
large $\kappa_{e}^{2}$ should be designed such that the electrodes and non-PE
outer protection layers have low Young’s modulus, and are as thin as possible.
But the Young’s modulus of non-PE material shouldn’t be too small such that
shearing occurs and two PE layers vibrate independently. Additionally,
increasing the coupling term A also leads to a large $\kappa_{e}^{2}$. If the
PE material is predetermined in the design, increasing $b_{PE}$ is a good
choice to increase A. While changing of W or L is not suggested since not only
A but also $C_{P}$ and $k$ are related to W and L. To increase $b_{PE}$, thick
non-PE material should be inserted in between two PE layers, and the device
center-line has minimal impact on spring constant.

Figure 6a shows the cross section of the commercial PPA2014. The non-PE
material in this device contributes 41.6% of the total spring constant, and
the $\kappa_{e}^{2}$ is calculated to be 0.053 (0.069 measured) with 5 mm
clamping overlap. A redesigned device is shown in Figure 6b. The redesigned
device retains the original L and W. To reduce $k_{non-PE}$, the Cu electrodes
are removed, and the FR4 outer layer is replaced by polyimide which has a
lower Young’s modulus. Thicker PE capacitors are adopted to increase $k_{PE}$.
In addition, a thicker center FR4 layer is used to increase $b_{PE}$ so that A
is increased. In the redesigned device, with a 5 mm clamping overlap, the non-
PE material contributes 5.6% of the spring constant, and $\kappa_{e}^{2}$ is
predicted to be increased by 2x to $\kappa_{e}^{2}$=0.121.

##  6 Conclusions

One factor that has prevented the widespread commercial use of PEHs is the
narrow range of vibration frequencies over which (near) optimum power can be
harvested. This paper shows that, with a combination of bias-flip electronics,
smart rectification, and a PEH optimized for large $\kappa_{e}^{2}$, one can
achieve a fractional 3-dB BW in excess of 15%.

## References

## References

  * [1] Guyomar D, Badel A, Lefeuvre E and Richard C 2005 IEEE T. Ultrason. Ferroelectr. Freq. Control. 52 584
  * [2] Ramadass Y and Chandrakasan A 2010 IEEE J. Solid—St. Circ. 45 189
  * [3] Zhao S, Ramadass Y, Lang J, Ma J and Buss D 2014 Proc. IEEE S3S Conf. 1
  * [4] Lefeuvre E, Badel A, Richard C and Guyomar D 2005 J. Intel. Mat. Syst. Str. 16 865
  * [5] Wu W, Wickenheiser A, Reissman T and Garcia E 2009 Smart Mater. Struct. 18 055012
  * [6] Zhao S, Paidimarri A, Ickes N, Araghchini M, Lang J, Ma J, Ramadass Y and Buss D 2015 Proc. IEEE S3S Conf. 1
  * [7] Hsieh P, Chen C and Chen H 2015 IEEE. T. Power. Electr. 30 3142
  * [8] Dutoit N, Wardle B and Kim S 2005 Integr. Ferroelectr. 71 121
  * [9] Shu Y and Lien I 2006 Smart Mater. Struct. 15 1499
  * [10] Erturk A and Inman D J 2009 Smart Mater. Struct. 18 025009