# Calculating magnetic interactions in organic electrides

Taek Jung Kim  Department of Physics, KAIST, 291 Daehak-ro, Yuseong-gu,
Daejeon 34141, Republic of Korea  Hongkee Yoon  Department of Physics, KAIST,
291 Daehak-ro, Yuseong-gu, Daejeon 34141, Republic of Korea  Myung Joon Han
<EMAIL_ADDRESS>Department of Physics, KAIST, 291 Daehak-ro, Yuseong-gu,
Daejeon 34141, Republic of Korea  KAIST Institute for the NanoCentury, Korea
Advanced Institute of Science and Technology, Daejeon 305-701, Korea

###### Abstract

We present our calculation results for organic magnetic electrides. In order
to identify the ‘cavity’ electrons, we use maximally-localized Wannier
functions and ‘empty atom’ technique. The estimation of magnetic coupling is
then performed based on magnetic force linear response theory. Both short- and
long-range magnetic interactions are calculated with a single self-consistent
calculation of a primitive cell. With this scheme we investigate four
different organic electrides whose magnetic properties have been partly
unknown or under debate. Our calculation results unveil the nature of magnetic
moment and their interactions, and justify or defy the validity of preassumed
spin models. Our work not only provides useful insight to understand magnetic
electrides but also suggests a new paradigm to study the related materials.

##  I Introduction

Electrides are a special type of ionic crystals. In these fascinating
materials the electrons, trapped in cavity, are the anions. Due to this
characteristic feature, many possibly useful properties are realized such as
high magnetic susceptibility, low work-function, and strong reducing
character, highly variable conductivity, low temperature thermionic emission,
high hyperpolarizability Dye (2009); Yanagi  _et al._ (2012); Buchammagari
_et al._ (2007); Xu  _et al._ (2007). Recently room-temperature-stable organic
and inorganic electride have been synthesized Satoru Matsuishi (2003); Redko
_et al._ (2005). Most of organic electrides are known to be antiferromagnetic
(AFM) from the field response Redko  _et al._ (2005); Dawes  _et al._ (1991);
Huang  _et al._ (1997); Xie  _et al._ (2000); Wagner  _et al._ (2000);
Ichimura  _et al._ (2002); Wagner and Dye (1995). In the sense that their
magnetic properties are presumably originated from the electrons in the
cavity, the magnetism of electrides is of unique interest. Although some
features of their magnetic interactions have been modeled, such as that the
electrons are interacting via a vacant aisle, and classified accordingly Dye
_et al._ (1996); Dye (1997); Ryabinkin and Staroverov (2010a, b, 2011), a
large part of their fundamental nature still remains elusive. It is largely
due to the limitation of conventional ab initio calculation method.

The conventional way of investigating magnetic property from first-principles
is to calculate the interaction parameter by comparing multiple total energies
corresponding to the ground state and meta-stable magnetic orders. In this
way, not only the ground state spin configuration but the magnetic interaction
strength are also calculated as shown recently by Dale and Johnson for
electrides Dale and Johnson (2016). However, this conventional approach is
severely limited when the system size is large which is indeed the case for
many organic electrides. For large systems, it is difficult to calculate the
long-range interactions as the supercell contains too many atoms. While the
magnetic interaction in solid is typically classified into the long-range
(e.g., Ruderman–Kittel–Kasuya–Yosida (RKKY) interactions) Ruderman and Kittel
(1954); Yosida (1957); Kasuya (1956) or short-range (e.g., superexchange
interactions) Kanamori (1960); Goodenough (1963) nature, the identification of
even such basic character has been hampered by the large unitcell size for
organic electrides. On top of their intriguing features of interacting path
presumably through some cavity aisle, this practical issue limits the ab
initio study.

In order to meet this challenge, here we introduce a new approach. First we
employ so-called ‘magnetic force response theory (MFT)’ Liechtenstein  _et
al._ (1987); Antropov  _et al._ (1996, 1997); Katsnelson and Lichtenstein
(2000); Han  _et al._ (2004); Yoon  _et al._ (2018) for calculating magnetic
interactions. MFT enables us to calculate all the magnetic interactions
residing in a given material within a primitive unitcell and at one time.
Thus, without a supercell, one can estimate the magnetic coupling parameter
$J$ as a function of distance for both short and long range. In order to
understand the magnetic properties of organic electrides and to demonstrate
the capability of our computation scheme, we take four different materials;
namely, Rb+(cryptand[2.2.2])e- (Fig. 2(a)), Li+(cryptand[2.1.1])e- (Fig.
2(b)), [Cs+(15C5)(18C6)e-]6(18C6) (Fig. 3(a)), and K+(cryptand[2.2.2])e- (Fig.
4(a)). Our calculations clearly show that the magnetic interactions in these
electrides indeed come from the localized electrons as anions. Further we
unveil their short-range versus long-range nature of the interactions. In
fact, for some electrides, there is an indication of oscillating $J$ which is
a signature of RKKY type magnetic couplings. In order to make MFT feasible,
one has to identify the trapped electron states properly. For this purpose, we
further employ maximally-localized Wannier function (MLWF) technique Marzari
and Vanderbilt (1997); Souza  _et al._ (2001). Another difficulty in dealing
with electrides within first-principles framework is about controlling
magnetic order for the cavity electrons. With a special constraint DFT scheme,
we successfully stabilized the magnetic solution of K+(cryptand[2.2.2])e- for
the first time. Our current work provides useful information to understand the
magnetism of organic electrides.

##  II  Computational methods

###  II.1 Magnetic force response theory

MFT is a method to calculate magnetic interactions at a given electronic
structure. In this method the exchange coupling is estimated as a response to
small spin tiltings as a perturbation from the given converged solution
Liechtenstein  _et al._ (1987); Yoon  _et al._ (2018):

$J_{ij}({\bf{q}})=\frac{1}{\pi}{\rm{Im}}\int\int^{\epsilon_{\rm{F}}}d{\rm{k}}\,d\epsilon\rm{\,Tr}[V_{{\bf{k}},i}^{\downarrow\uparrow}{\mathbf{G}}_{{\bf
k},ij}^{\uparrow\uparrow}(\epsilon)V_{{\bf{k+q}},j}^{\uparrow\downarrow}{\mathbf{G}}_{{\bf
k+q},ji}^{\downarrow\downarrow}(\epsilon)].$ (1)

Here i and j are the site indices, and up and down arrow indicate the spin
direction. Green’s function $\mathbf{G}$ is represented as

$\mathbf{G}^{\uparrow\uparrow}_{{\bf{k}},ij}(\epsilon)=\sum_{n}\frac{\ket{\psi_{{\bf{k}},i}^{\uparrow}}\bra{\psi^{\uparrow}_{{\bf{k}},j}}}{\epsilon-\epsilon_{\uparrow
n,{\bf{k}}}+i\eta}$ (2)

where $\epsilon_{n,{\bf{k}}}$ and $\ket{\psi_{{\bf{k}},i}}$ refers to the n-th
eigenvalue and eigenstate, respectively. $\mathbf{V}$ is given by

$V_{{\bf{k}},i}^{\downarrow\uparrow}={\frac{1}{2}}({\bf{H}}_{{\bf
k},i}^{\downarrow\downarrow}-{\bf{H}}_{{\bf k},i}^{\uparrow\uparrow})$ (3)

where ${\bf{H}}_{{\bf k},i}^{\uparrow\uparrow(\downarrow\downarrow)}$ is Kohn-
Sham Hamiltonian corresponding to the collinear up(down) spin.

The exchange interaction between two sites is calculated by Fourier
transformation from $\bf{k}$ to real space. Thus one can just take the minimal
size of unitcell with no need to consider large supercells. Once the localized
magnetic sites are well defined, MFT provides the exchange constants, $J$’s,
as a function of distance. Note that in MFT we extract all the information
from a single self-consistently converged electronic structure. The further
details of our implementation and the results of some classical example can be
found in our previous studies Han  _et al._ (2004); Yoon  _et al._ (2018).

###  II.2 Calculation details

We perform density functional theory (DFT) calculations within generalized
gradient approximation proposed by Perdew et al. (GGA-PBE) Perdew  _et al._
(1996) by employing LCPAO (linear combination of pseudo-atomic orbitals)
method Ozaki (2003); Ozaki and Kino (2004) as implemented in our ‘OpenMX’
software package ope . It should be noted that the limitation of GGA in
describing the correlation effect can cause the overestimation of magnetic
couplings or the underestimation of the charge occupation Han  _et al._
(2004); Solovyev and Terakura (1998); Johnson  _et al._ (2013). 3×3×3 k-points
and 500 Ry energy cutoff are used for numerical integration. Poisson equations
are solved by using fast Fourier transformations, and the projector expansion
method is used to accurately calculate three-center integrals associated with
the deep neutral atom potential Ozaki and Kino (2005).

Figure 1:  (Color online) (a), (b) The calculated band dispersion of (a)
Rb+(cryptand[2.2.2])e- and (b) Li+(cryptand[2.1.1])e- (blue line). The
calculated MLWF band dispersion is expressed in red lines.

In order to describe the characteristic feature of electride, namely the
electrons trapped in the cavity space, we employ the ‘empty atom’ technique
for all of our systems. The empty atom technique, also called as ‘ghost atom’
or ‘empty sphere’, has been used, for example, to correctly estimate the basis
set superposition error (BSSE) Latajka and Scheiner (1989); Simon  _et al._
(1996), and to treat a large void space Becker  _et al._ (1996); Kalpana  _et
al._ (1996); Ossicini and Bisi (1989); Sque  _et al._ (2007) within local
basis schemes. In this method, an empty atom is represented as an atom with a
nuclear charge of zero, which acts as a basis for describing the wave function
of the empty space. We used two $s$, one $p$, and one $d$ orbitals with a
cutoff radius of 11 a.u. as the basis of empty atom. We determine the number
of empty atoms by comparing the band structure with the plane-wave result Dale
_et al._ (2014). MLWF Marzari and Vanderbilt (1997); Souza  _et al._ (2001) is
also used to identify the magnetic bands as a localized state. We found that
the position of Wannier functions was well compared with the distances known
from experiments Huang  _et al._ (1997); Xie  _et al._ (2000); Wagner and Dye
(1995); Dye (1990). This combination of empty atom and MLWF techniques enables
us to perform the MFT calculation to estimate the magnetic couplings. The
s-wave symmetry for Wannier functions is considered to describe the cavity-
electron states since they are clearly of s-orbital character from spin
density plots seen in Fig. 2(a) and 2(b) for example. Figure 1(a) and 1(b)
shows the calculated band dispersion of Rb+(cryptand[2.2.2])e- and
Li+(cryptand[2.1.1])e-, respectively (blue line). The MLWF band is expressed
by red lines. An excellent overlap of the two bands (blue and red) shows that
MLWF well identifies the electronic states in the cavity.

MFT calculations are conducted based on the previously-known magnetic phase
Dale and Johnson (2016), namely, G-type AFM order (in which all of the
directly connected neighbors have different spins) except for
[Cs+(15C5)(18C6)e-]6(18C6). For [Cs+(15C5)(18C6)e-]6(18C6), we consider the
intra-ring AFM order (see Fig. 3(a)). In order to get the AFM phase for all
other electrides except K+(cryptand[2.2.2])e-, the initial spin polarization
is applied to the surrounding hydrogen atoms as in the previous calculation
Dale and Johnson (2016).

For K+(cryptand[2.2.2])e-, however, even such an ad hoc technique does not
work as reported in the previous theoretical study Dale and Johnson (2016). As
a result, the magnetic property of this electride has never been theoretically
addressed. Here, we used a kind of special constraint DFT technique to
stabilize the magnetic order. In this scheme, the initial spin density is
assigned onto the ‘empty atom’ sites which means that the constraint is
directly imposed on the cavity. We further confirmed that this magnetic
solution is not just consistent with experiment but also robust through the
further self-consistent steps. It is eventually converged into the magnetic
solution without any constraint. From this process, we successfully obtained
the well-stabilized self-consistent magnetic solution of K+(cryptand[2.2.2])e-
for the first time, on top of which MFT can be conducted.

##  III RESULTS AND DISCUSSION

###  III.1 Rb+(cryptand[2.2.2])e- and Li+(cryptand[2.1.1])e-

Rb+(cryptand[2.2.2])e- and Li+(cryptand[2.1.1])e- are the electrides with
‘ladder-like’ channel Huang  _et al._ (1997); Xie  _et al._ (2000).
Experimentally, the magnetic property has been studied with electron
paramagnetic resonance (EPR) spectroscopy as a function of temperature, and
fitting to a certain exchange model. For these two materials, the EPR data is
well fit to ‘FN (first-neighbor)-1D Heisenberg’ model. However, a recent
calculation study Dale and Johnson (2016) raises a question against this
conclusion com . The calculated $J$ from the total energy difference mapped
onto the same FN-1D model is found to be quite different from the experimental
data Huang  _et al._ (1997); Xie  _et al._ (2000). It can be argued that this
difference is attributed to the next and longer-range magnetic couplings Dale
and Johnson (2016). Importantly, however, a solid conclusion could not be made
because the next neighbor interactions were not accessible within the
conventional computation scheme.

Figure 2:  (Color online) (a), (b) The calculated spin density of (a)
Rb+(cryptand[2.2.2])e- and (b) Li+(cryptand[2.1.1])e- where red and green
spheres represent the up and down spin density, respectively. We used the
isosurface value of 0.00075 in atomic unit (a.u.). $J_{1}$, $J_{2}$, and
$J_{3}$ refer to the first, second, and third neighbor interactions in (a) and
(b). (c), (d) The calculated magnetic coupling parameters for (c)
Rb+(cryptand[2.2.2])e- and (d) Li+(cryptand[2.1.1])e-. Our calculation results
by MFT (dark blue circles) are compared with the previous calculation by total
energy difference (green triangles; Ref. Dale and Johnson, 2016) and
experiment (magenta squares; Ref. Huang  _et al._ , 1997; Xie  _et al._ ,
2000). Note that both short and long range interactions are calculated from
MFT while only nearest neighbor values can be obtained from experiments and
total-energy-based computation scheme.

Hereby using MFT combined with MLWF, we calculated the exchange coupling
constants as a function of distance; from the nearest neighbor to the long-
range interactions. The results are presented in Figs. 2(a)-(d). For
Rb+(cryptand[2.2.2])e- (Fig. 2(a) and 2(c)), MFT calculation shows that the
first neighbor interaction is dominant, $J_{1}/k_{B}$=56.3 K, while the second
and third neighbor interaction are much smaller, $J_{2}/k_{B}$=0.17 K and
$J_{3}/k_{B}$= $-$4.56 K (see Table 1). Thus our results confirm that
Rb+(cryptand[2.2.2])e- is well classified as a FN-1D Heisenberg system.

On the other hand, Li+(cryptand[2.1.1])e- is not the case (Fig. 2(b) and
2(d)). The calculated first neighbor interaction is $J_{1}/k_{B}$=51.7 K, and
the second and third neighbor value is quite significant; $J_{2}/k_{B}$=6.69 K
and $J_{3}/k_{B}$ =14.0 K (see Table 1). Note that this corresponds to 12.9%
and 27.1% of $J_{1}$, respectively. Thus it is hardly regarded as a FN-1D spin
system contrary to the previous study.

Our results demonstrate the usefulness of the current computation scheme for
magnetic electrides. The MFT in combination with MLWF provides the long-range
interaction parameters and therefore one can have a reliable picture for the
magnetism even if the system size is too large to be calculated by
conventional total energy method. Further, our calculation confirms that the
magnetic properties measured in the experiment indeed come from the localized
electron state in the cavity sites. It is because MFT estimates the magnetic
force in between two sites, and in the current case, these ‘sites’ are defined
as the cavity electrons by means of empty atoms and MLWF technique. We also
emphasize that, in our scheme, one does not need to build any a priori model
considering only a few neighbor interactions. The exchange parameters are
calculated as the response to spin tilting. This feature is advantageous when
one needs to consider any type of model building or to justify a certain
model.

###  III.2 [Cs+(15C5)(18C6)e-]6(18C6)

Figure 3:  (Color online) (a) The calculated spin density of
[Cs+(15C5)(18C6)e-]6(18C6) where red and green spheres represent the up and
down spin density, respectively. We used the isosurface value of 0.00075 a.u..
The spin order in this material can be referred as ‘intra-ring
antiferromagnetic (intra-ring AFM)’ configuration: Black lines show the
connections of the intra-vacant spaces of ‘six-memebered rings’ while yellow
lines show the connections between the different ‘six-memebered rings’. Yellow
points represent the connections which are invisible in the current view of
Figure. Note that each site is connected by two black lines and two yellow
lines or points. In terms of distance, the black lines correspond to the first
neighbor interactions ($J_{1}$) and the yellow dots/lines to the second
neighbors ($J_{2}$). The stacking sequence of the six-membered rings is A-B-
C-A-… from the side view. (b) The calculated exchange interaction for
[Cs+(15C5)(18C6)e-]6(18C6). We estimated $J$ values based on two different
spin density configurations (see main text for more details), and the results
are presented with dark-blue circles and pink diamonds. For comparison, the
previous calculation (Ref. Dale and Johnson, 2016) and experimental (Ref.
Wagner and Dye, 1995) data are also presented. The two inset figures show the
schematic spin configuration based on which we performed the MFT calculation
for $J$; namely, intra-ring AFM and ferromagnetic (FM) order.

[Cs+(15C5)(18C6)e-]6(18C6) is an electride having ‘six-membered ring’ cluster
structure as shown in Fig. 3(a), and its magnetic susceptibility was reported
by Wagner and Dye Wagner and Dye (1995). We note that, in this experimental
study, the measured data has been analyzed based on a simple model assumption
that the intra-ring couplings are the only magnetic interactions Wagner and
Dye (1995). From this analysis, it was concluded that this electride has
strong AFM interactions.

The theoretical study was also conducted under the same assumption that only
the first-neighbor interaction is important Dale and Johnson (2016). It is
presumably because the computational cost is too large to calculate multiple
total energies within the enlarged supercell geometry. In fact, in order to
simulate the theoretically-estimated G-type spin ground state, one needs to
calculate the 4080-atom supercell. In Ref. Dale and Johnson, 2016, the cell
size of 510 atoms has been used within which only the nearest neighboring $J$
is accessible. This result is plotted in Fig. 3(b) (denoted by green triangle)
showing that the difference between the previous calculation and experiment is
about 395 K. The origin of the relatively large difference is unclear. It was
indeed concluded in Ref. Dale and Johnson, 2016 that this difference is
attributed to the supercell size effect.

In this context, MFT can give a useful insight since it does not require any
supercell calculation but it does provide the longer range interactions. Fig.
3(b) shows our MFT calculation results of exchange couplings as a function of
distance (see dark-blue circles and pink diamonds). Since the magnetic ground
state cannot be represented within the structural primitive cell, we
considered two different magnetic solutions which can be realized within this
primitive cell; namely, FM (pink diamonds) and AFM (dark-blue circles). Here
AFM order refers to the AFM order in between the intra-rings (see the inset of
Fig. 3(b)). We first note that the two calculation results based on FM and AFM
spin density are quite similar; the difference is less than 6.5 K. It not only
indicates that the spin-polarized electron states in cavity are basically well
localized as previously discussed Yoon  _et al._ (2018), but also implies that
we do not need the real spin ground state density in order to estimate $J$
values Yoon  _et al._ (2018).

$J_{n}/k_{B}$ Name | Rb+(cryptand[2.2.2])e- | Li+(cryptand[2.1.1])e- |  [Cs+(15C5)(18C6)e-]6(18C6) | K+(cryptand[2.2.2])e-
---|---|---|---|---
|  MFT |  $\Delta$E | Exp |  MFT |  $\Delta$E | Exp |  MFT |  $\Delta$E | Exp |  MFT |  $\Delta$E | Exp
$J_{1}$ | 56.3 | 78.2 | 30 | 51.7 | 177 | 54 | 19.4 | 15.2 | 410 | 9.6 | - | 440
$J_{2}$ | 0.2 | - | - | 6.7 | - | - | 24.9 | - | - | 5.7 | - | -
$J_{3}$ |  $-$4.6 | - | - | 14.0 | - | - | 6.6 | - | - | 148 | - | -
$J_{7}$ | 0.1 | - | - | 0.02 | - | - |  $-$1.3 | - | - |  $-$33.5 | - | -
Table 1:  Exchange parameters of 4 different organic electrides in the unit of
K. Our calculation results by MFT are compared with the previous calculations
($\Delta$E; Ref. Dale and Johnson (2016)) and experiments (Exp). The
experimental data can be found in Ref. Xie  _et al._ (2000), Ref. Huang  _et
al._ (1997), Ref. Wagner and Dye (1995) and Ref. Ichimura  _et al._ (2002) for
Rb+(cryptand[2.2.2])e-, Li+(cryptand[2.1.1])e-, [Cs+(15C5)(18C6)e-]6(18C6) and
K+(cryptand[2.2.2])e-, respectively.

The sum of all magnetic interactions is $J_{tot}/k_{B}$ = ($J_{1}$+$J_{2}$+…
)$/k_{B}$= 53 K which is notably smaller than the experimental value of 410 K
(magenta square) Wagner and Dye (1995) (see Table 1). It means that the
difference between the previous calculation and experiment is not originated
from the longer-range interactions or supercell-size effects. Comparing our
MFT result with the previous total energy-based estimation, the first neighbor
$J_{1}$ is in good agreement with each other within 4.2 K. Importantly, the
second neighbor $J_{2}$ is comparable with and slightly larger than $J_{1}$.
Namely, the inter-ring interaction is sizable and the nearest neighbor spin
model is not relevant to this material. Our calculations clearly show that the
inter-ring coupling needs to be taken into account.

Further, our calculation confirms that the spin ground state of this material
is indeed G-type. As mentioned above, this material was speculated to have
G-type AFM ground state Dale and Johnson (2016). However, there is no
experimental or theoretical evidence for that. Theoretical confirmation has
been hampered by the large supercell size of 4080 atoms. Our MFT results of
Fig. 3(b) clearly shows that $J_{1}$ and $J_{2}$ are the two dominant
couplings and both of them are AFM.

###  III.3 K+(cryptand[2.2.2])e-

Figure 4:  (Color online) (a) The calculated spin density of
K+(cryptand[2.2.2])e- where the red and green spheres represent the up and
down spins, respectively. We used the isosurface value of 0.0003 a.u.. (b) The
calculation results of magnetic interaction as a function of distance. The
magenta square shows the experimental result (Ref. Ichimura  _et al._ , 2002)
while the dark-blue circles are our calculation results. For this material,
there is no previous calculation result because of the difficulty in
stabilizing the magnetic solution (see the main text for more details).

Our final example, K+(cryptand[2.2.2])e-, is also known to be an AFM ordered
organic electride Ichimura  _et al._ (2002). For this material, however, there
has been no successful calculation in obtaining the magnetic solution. The
conventional DFT calculation gives the paramagnetic spin ground state as the
converged solution even if starting from the spin polarized initial condition.
Even with ad hoc treatment in which initial spins are assigned to the hydrogen
atoms around the vacant space, the magnetic solution is hardly achieved Dale
and Johnson (2016).

First of all, we successfully obtained the magnetic solution by applying
initial spins directly on the vacant space. We developed a constrained DFT
scheme for assigning the initial spin moment to ‘empty sphere’. This process
is not straightforward in the sense that electron spins need to be polarized
within the empty spheres. We applied the magnetic constraint during the
initial $~{}20$ self-consistent steps, through which empty spheres are
occupied by polarized electrons. After that, the usual DFT self-consistent
calculations are performed with the constraint turned off. With this scheme,
we successfully generated the well-stabilized AFM solution as shown in Fig.
4(a).

On top of this spin density we performed MFT calculation and successfully
estimated magnetic interactions for the first time. Our results are summarized
in Fig. 4(b). First of all, our calculation confirms that this magnetic phase
is indeed G-type AFM ground state as speculated in the previous study Dale and
Johnson (2016). It is also consistent with an experiment Ichimura  _et al._
(2002). The largest interaction is the third neighbor $J_{3}$ which is notably
larger than $J_{1}$, $J_{2}$ and others. Interestingly the second largest
interaction is $J_{7}$ which is about 25% of the $J_{3}$. Our results is
consistent with the experiment Ichimura  _et al._ (2002) in the sense that the
magnetic property of this material can be described with two parameters in two
dimension. The total sum of all our $J/k_{B}$’s is about 110 K (see Table 1).

##  IV Summary

A new theoretical approach is applied to study magnetic electrides. Spin-
polarized electrons trapped in the cavity are identified by empty atom and
MLWF method, and their interactions calculated within MFT. The usefulness of
this scheme is shown by calculating four different organic electrides for
which the validity of pre-assumed models have remained unclear. The long range
magnetic interaction profile as a function of distance is calculated and
compared, which has not been accessed by the conventional total energy
calculations. For K+(cryptand[2.2.2])e-, we apply a constraint DFT method to
stabilize the magnetic solution and calculate the magnetic interaction for the
first time. Our study provides useful insights to understand magnetic
electrides and related materials.

##  V Acknowledgments

This work was supported by Basic Science Research Program through the National
Research Foundation of Korea (NRF) funded by the Ministry of Education
(2018R1A2B2005204).

## References

  * Dye (2009) J. L. Dye, Accounts of Chemical Research 42, 1564 (2009).
  * Yanagi  _et al._ (2012) H. Yanagi, T. Kuroda, K.-B. Kim, Y. Toda, T. Kamiya, and H. Hosono, Journal of Materials Chemistry 22, 4278 (2012).
  * Buchammagari  _et al._ (2007) H. Buchammagari, Y. Toda, M. Hirano, H. Hosono, D. Takeuchi, and K. Osakada, Org. Lett. 9, 4287 (2007).
  * Xu  _et al._ (2007) H.-L. Xu, Z.-R. Li, D. Wu, B.-Q. Wang, Y. Li, F. L. Gu, and Y. Aoki, Journal of the American Chemical Society 129, 2967 (2007).
  * Satoru Matsuishi (2003) H. H. Satoru Matsuishi, Science 301, 626 (2003).
  * Redko  _et al._ (2005) M. Y. Redko, J. E. Jackson, R. H. Huang, and J. L. Dye, Journal of the American Chemical Society 127, 12416 (2005).
  * Dawes  _et al._ (1991) S. B. Dawes, J. L. Eglin, K. J. Moeggenborg, J. Kim, and J. L. Dye, Journal of the American Chemical Society 113, 1605 (1991).
  * Huang  _et al._ (1997) R. H. Huang, M. J. Wagner, D. J. Gilbert, K. A. Reidy-Cedergren, D. L. Ward, M. K. Faber, and J. L. Dye, J. Am. Chem. Soc. 119, 3765 (1997).
  * Xie  _et al._ (2000) Q. Xie, R. H. Huang, A. S. Ichimura, R. C. Phillips, W. P. Pratt, and J. L. Dye, J. Am. Chem. Soc. 122, 6971 (2000).
  * Wagner  _et al._ (2000) M. J. Wagner, A. S. Ichimura, R. H. Huang, R. C. Phillips, and J. L. Dye, The Journal of Physical Chemistry B 104, 1078 (2000).
  * Ichimura  _et al._ (2002) A. S. Ichimura, M. J. Wagner, and J. L. Dye, The Journal of Physical Chemistry B 106, 11196 (2002).
  * Wagner and Dye (1995) M. J. Wagner and J. L. Dye, Journal of Solid State Chemistry 117, 309 (1995).
  * Dye  _et al._ (1996) J. L. Dye, M. J. Wagner, G. Overney, R. H. Huang, T. F. Nagy, and D. Tománek, J. Am. Chem. Soc. 118, 7329 (1996).
  * Dye (1997) J. L. Dye, Inorganic Chemistry 36, 3816 (1997).
  * Ryabinkin and Staroverov (2010a) I. G. Ryabinkin and V. N. Staroverov, Phys. Rev. A 82, 022505 (2010a).
  * Ryabinkin and Staroverov (2010b) I. G. Ryabinkin and V. N. Staroverov, Phys. Rev. A 81, 032509 (2010b).
  * Ryabinkin and Staroverov (2011) I. G. Ryabinkin and V. N. Staroverov, Physical Chemistry Chemical Physics 13, 21615 (2011).
  * Dale and Johnson (2016) S. G. Dale and E. R. Johnson, Phys. Chem. Chem. Phys. 18, 27326 (2016).
  * Ruderman and Kittel (1954) M. A. Ruderman and C. Kittel, Physical Review 96, 99 (1954).
  * Yosida (1957) K. Yosida, Physical Review 106, 893 (1957).
  * Kasuya (1956) T. Kasuya, Progress of Theoretical Physics 16, 45 (1956).
  * Kanamori (1960) J. Kanamori, Journal of Applied Physics 31, S14 (1960).
  * Goodenough (1963) J. B. Goodenough, “Magnetism and the Chemical Bond by Goodenough, J B: Interscience Publishers Hardcover - HPB-Dallas,” (1963).
  * Liechtenstein  _et al._ (1987) A. I. Liechtenstein, M. I. Katsnelson, V. P. Antropov, and V. A. Gubanov, Journal of Magnetism and Magnetic Materials 67, 65 (1987).
  * Antropov  _et al._ (1996) V. P. Antropov, M. I. Katsnelson, B. N. Harmon, M. Van Schilfgaarde, and D. Kusnezov, Physical Review B 54, 1019 (1996).
  * Antropov  _et al._ (1997) V. P. Antropov, M. I. Katsnelson, and A. I. Liechtenstein, Physica B: Condensed Matter Proceedings of the Yamada Conference XLV, the International Conference on the Physics of Transition Metals, 237, 336 (1997).
  * Katsnelson and Lichtenstein (2000) M. I. Katsnelson and A. I. Lichtenstein, Physical Review B 61, 8906 (2000).
  * Han  _et al._ (2004) M. J. Han, T. Ozaki, and J. Yu, Physical Review B 70 (2004), 10.1103/PhysRevB.70.184421.
  * Yoon  _et al._ (2018) H. Yoon, T. J. Kim, J.-H. Sim, S. W. Jang, T. Ozaki, and M. J. Han, Physical Review B 97 (2018), 10.1103/PhysRevB.97.125132.
  * Marzari and Vanderbilt (1997) N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
  * Souza  _et al._ (2001) I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Ozaki (2003) T. Ozaki, Phys. Rev. B 67, 155108 (2003).
  * Ozaki and Kino (2004) T. Ozaki and H. Kino, Phys. Rev. B 69, 195113 (2004).
  * (35) http://www.openmx-square.org/ .
  * Solovyev and Terakura (1998) I. V. Solovyev and K. Terakura, Physical Review B 58, 15496 (1998).
  * Johnson  _et al._ (2013) E. R. Johnson, A. Otero-de-la Roza, and S. G. Dale, The Journal of Chemical Physics 139, 184116 (2013).
  * Ozaki and Kino (2005) T. Ozaki and H. Kino, Physical Review B 72 (2005), 10.1103/PhysRevB.72.045121.
  * Latajka and Scheiner (1989) Z. Latajka and S. Scheiner, Journal of Molecular Structure: THEOCHEM 199, 9 (1989).
  * Simon  _et al._ (1996) S. Simon, M. Duran, and J. J. Dannenberg, The Journal of Chemical Physics 105, 11024 (1996).
  * Becker  _et al._ (1996) U. Becker, M. F. Hochella, and E. Apra, American Mineralogist 81, 1301 (1996).
  * Kalpana  _et al._ (1996) G. Kalpana, B. Palanivel, R. M. Thomas, and M. Rajagopalan, Physica B: Condensed Matter 222, 223 (1996).
  * Ossicini and Bisi (1989) S. Ossicini and O. Bisi, Surface Science 211-212, 572 (1989).
  * Sque  _et al._ (2007) S. J. Sque, R. Jones, and P. R. Briddon, physica status solidi (a) 204, 3078 (2007).
  * Dale  _et al._ (2014) S. G. Dale, A. Otero-de-la Roza, and E. R. Johnson, Phys. Chem. Chem. Phys. 16, 14584 (2014).
  * Dye (1990) J. L. Dye, Science 247, 663 (1990).
  * (47) “Ref. Dale and Johnson (2016) also mentions that there are materials in which their calculations and experiments are well matched. Among them, we proceed further calculation for Cs+(15C5)2e- to check reliability of our methodology intensively. As a result, Ref. Dale and Johnson (2016) ($j_{1}$ = 12.2K) and our results ($j_{1}$ = 10.8K) have agreed well. In addition, we confirm that our calculation result agree well with experiment in the sense that it can be represented as 1D-chain model as $j_{1}$ is 10 times larger than other $j$ values.” .


# Signatures of Non-Markovianity in Cavity-QED with Color Centers in 2D
Materials

Mark Kamper Svendsen  CAMD, Department of Physics, Technical University of
Denmark, 2800 Kgs. Lyngby, Denmark  Sajid Ali  CAMD, Department of Physics,
Technical University of Denmark, 2800 Kgs. Lyngby, Denmark  Nicolas Stenger
Department of electrical and photonics engineering, Technical University of
Denmark, 2800 Kgs. Lyngby, Denmark  Center for Nanostructured Graphene,
Technical University of Denmark, 2800 Kgs. Lyngby, Denmark  NanoPhoton –
Center for Nanophotonics, Technical University of Denmark, 2800 Kgs. Lyngby,
Denmark  Kristian Sommer Thygesen  CAMD, Department of Physics, Technical
University of Denmark, 2800 Kgs. Lyngby, Denmark  Center for Nanostructured
Graphene, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark  Jake
Iles-Smith<EMAIL_ADDRESS>Department of Physics and
Astronomy, The University of Manchester, Oxford Road, Manchester M13 9PL,
United Kingdom  Department of Electrical and Electronic Engineering, The
University of Manchester, Oxford Road, Manchester M13 9PL, United Kingdom

(January 2022)

###### Abstract

Light-matter interactions of defects in two dimensional materials are expected
to be profoundly impacted by strong coupling to phonons. In this work, we
combine _ab initio_ calculations of a defect in hBN, with a fully quantum
mechanical and numerically exact description of a cavity-defect system to
elucidate this impact. We show that even at weak light-matter coupling, the
dynamical evolution of the cavity-defect system has clear signatures of non-
markovian phonon effects, and that linear absorption spectra show the
emergence of hybridised light-matter-phonon states in regimes of strong light-
matter coupling. We emphasise that our methodology is general, and can be
applied to a wide variety of material/defect systems.

††preprint: APS/123-QED††preprint: APS/123-QED

Single-photon emission has been observed from a broad range of two dimensional
materials (2DM) [1, 2, 3, 4, 5, 6, 7]. Of particular interest to quantum
technologies is the emission from defect complexes in hexagonal Boron Nitride
(hBN) [8, 9, 10, 11, 12, 13], owing to the materials wide band-gap [14] and
stability of the emitters [10, 15]. Not only does the precise nature of these
defect complexes remain under scrutiny [10, 11, 16], but also their quantum
optical properties remain largely unexplored, particularly in regimes of
strong light-matter coupling [17, 18]. Crucially, the optical properties of
condensed matter systems cannot be understood through the standard theory of
quantum optics, since strong system-environment interactions lead to a
breakdown of the underpinning peturbative methods [19, 20]. Of particular
relevance to the optical properties of defect complexes in 2DM are strong
electron-phonon interactions. Specifically, defect emission spectra typically
show sharp and well resolved phonon sidebands [21, 22, 23, 24], a hallmark of
vibronic state formation and long lived correlations between electronic and
vibrational degrees of freedom. It is therefore crucial that any theory
describing the dynamical or optical behaviour of a 2DM defect complex,
accurately accounts for electron-phonon interactions.

To this end, we study the cavity quantum electrodynamics (cQED) of defects in
2DM [17, 18], focusing on two colour centres in hBN proposed as single-photon
emitters, the CBVN [25, 10] and C2CN defect complexes [26, 27]. We accurately
account for electron-phonon interactions by combining _ab initio_ calculations
of the defects, with a fully quantum mechanical and numerically exact
description of the dynamics of the defect states interacting with a single
mode optical cavity. The former provides an atomistic description of the
phonon modes of the material system. The latter employs the time evolved
matrix product operator (TEMPO) algorithm [28], where non-Markovian influences
of the phonon environment are encoded within a tensor network. When combined
with tensor compression methods, TEMPO provides an efficient approach for
calculating the reduced dynamics of the system [29, 30, 31], which we extend
to extract the linear response spectrum of the defect-cavity system. This
hybrid approach is inspired by [32] and allows us to maintain the accuracy and
insight of first-principles calculations, while enabling coherent dynamics to
be captured in previously inaccessible regimes of light-matter coupling,
information that is typically inaccessible to first-principle methods.

A central finding of this work is that electron-phonon interactions in the
studied defect complexes lead to highly non-Markovian dynamics, even in
regimes of moderate light-matter coupling strengths, which are accessible to
state-of-the-art experiments. We attribute this behaviour to coherent coupling
between the electronic states of the defect and high-quality ($Q$) factor
phonon modes, leading to long-lived correlations between the system and its
phonon environment. This provides a mechanism to optically manipulate
mechanical modes of a condensed matter system that would otherwise not
directly couple to light. Furthermore, from linear response spectra, we see
evidence of hybridisation between cavity-defect polaritons and high-$Q$ phonon
modes when entering regimes of strong light-matter coupling. These states bare
resemblance to recently predicted exciton-photon-phonon hybridisation in hBN
[33], however in this instance they occur in a continuum of modes coupled to a
defect state, and therefore will inherit non-linear features from the
localised electronic states.

Figure 1: (a) A schematic figure of the CBVNdefect in hBN interacting with a
single quantised field mode. This mode could be plasmonic or photonic in
nature. (b) A plot of the _ab initio_ (solid) and effective width (dashed)
spectral densities. Using the atomistic simulation, we ascribe the dominant
peak at $\sim 125$ meV to a phonon breathing mode. (c) shows the bath
correlation function for the _ab intio_ (solid) and effective width (dashed)
spectral densities. The sharp peaks of the _ab initio_ spectral density leads
to long lived oscillations in the correlation function, which are washed out
when an effective width approximation is used.

For simplicity we focus on the CBVNdefect complex in the manuscript, leaving
the discussion of C2CN to the supplementary information (SI). The electronic
structure of the CBVNdefect consists of two separate manifolds [34]: the first
consists of the $(2)^{1}$A1 and $(1)^{1}B_{1}$ excited states, and a
$(1)^{1}A_{1}$ singlet ground state; the second contains the triplet states
$(1)^{3}\mathrm{B}_{1}$ and $(2)^{3}\mathrm{B}_{1}$ [16], separated in energy
by $\hbar\omega_{e}=2$ eV [34]. By assuming the inter-system crossing due to
spin–orbit coupling occurs on a much longer timescale ($10-100$ ns) than
cavity-enhanced optical transitions ($<1$ ns), we restrict our attention to
the triplet state manifold, which reduces to a two level system with ground
and exited states, $\ket{g}$ and $\ket{e}$ respectively 111This analysis is
repeated in the supplementary information for the C2CN defect.. The electronic
transition interacts with a single mode optical cavity as shown schematically
in Fig. 1a, with a cavity resonance $\hbar\Omega_{c}$ and described by a
Jaynes-Cummings type interaction:

$H_{\mathrm{S}}=\hbar\omega_{e}\sigma^{\dagger}\sigma+\hbar
g(\sigma^{\dagger}a+\sigma a^{\dagger})+\hbar\Omega_{c}a^{\dagger}a,$ (1)

where $\sigma=\ket{g}\\!\bra{e}$ is the transition dipole operator, $a$
($a^{\dagger}$) is the annihilation (creation) operator of the cavity mode,
and $g$ is the light-matter coupling strength. Here we limit ourselves to the
single excitation subspace, spanned by the basis
$\\{\ket{g,0},\ket{e,0},\ket{g,1}\\}$. The composite cavity-emitter system
interacts with both vibrational and electromagnetic environments, such that
the global Hamiltonian is written as
$H=H_{\mathrm{S}}+H_{I}^{\mathrm{EM}}+H_{I}^{\mathrm{Ph}}+H_{\mathrm{B}}^{\mathrm{EM}}+H_{\mathrm{B}}^{\mathrm{Ph}}$.

The emitter-cavity system couples to an external electromagnetic environment
either through direct emission from the defect into free space, or leakage
from the cavity mode. In the rotating wave approximation the interaction to
the electromagnetic field takes the form
$H_{I}^{\mathrm{EM}}=\sum_{j}\sum_{\mathbf{l}}(\hbar
f_{j,\mathbf{l}}\hat{A}_{j}^{\dagger}c_{j,\mathbf{l}}+\mathrm{h.c.}),$ where
$\hat{A}_{1}=\sigma$ and $\hat{A}_{2}=a$, and $c_{j,\mathbf{l}}$ is the
annihilation operator for the $\mathbf{l}^{\mathrm{th}}$ mode of the
$j^{\mathrm{th}}$ electromagnetic environment. In the limit of weak coupling
to these fields, we can make a flat spectrum approximation, such that
$f_{1,\mathbf{1}}=\sqrt{2\pi\Gamma}$ and $f_{2,\mathbf{1}}=\sqrt{2\pi\kappa}$
[36, 37]. Here $\Gamma$ is the spontaneous emission rate from the TLE into
free space, and $\kappa$ is the cavity leakage rate. Using these
approximations and tracing over the external EM fields, we can describe the
optical contribution to the system evolution through the superoperator
$\mathcal{V}_{t}=\exp(t~{}\mathcal{L}_{\mathrm{S}})$, where
$\mathcal{L}_{\mathrm{S}}\rho=-i[H_{\mathrm{S}},\rho]+\Gamma
L_{\sigma}\rho+\kappa L_{a}\rho$, and we have introduced the Lindblad
dissipators $L_{o}\rho=o\rho o^{\dagger}-\\{o^{\dagger}o,\rho\\}/2$. For a
full derivation and discussion of these approximations please refer to the SI.

We take a linear electron-phonon interaction [38] of the form
$H_{\mathrm{I}}^{\mathrm{Ph}}=\sigma^{\dagger}\sigma\sum_{\mathbf{k}}\hbar
g_{\mathbf{k}}(b_{\mathbf{k}}^{\dagger}+b_{-\mathbf{k}}),$ where
$b_{\mathbf{k}}$ is the annihilation operator for a phonon mode with wave
vector $\mathbf{k}$, with the corresponding free field evolution
$H_{\mathrm{B}}^{\mathrm{Ph}}=\sum_{\mathbf{k}}\hbar\omega_{\mathbf{k}}b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}$,
and $g_{\mathbf{k}}$ is the electron-phonon coupling strength of the
$\mathbf{k}^{\mathrm{th}}$ phonon mode. The electron-phonon coupling can be
fully characterised in terms of the spectral density
$\mathcal{J}_{\mathrm{Ph}}(\omega)=\sum_{\mathbf{k}}S_{\mathbf{k}}\delta(\omega-\omega_{\mathbf{k}})$,
where $S_{\mathbf{k}}=\omega_{\mathbf{k}}^{-2}|{}g_{\mathbf{k}}|^{2}$ is the
partial Huang-Rhys factor of mode $\mathbf{k}$. The total Huang-Rhys parameter
can then be recovered by integrating the spectral density through
$S_{\mathrm{Tot}}=\int_{0}^{\infty}d\omega\mathcal{J}_{\mathrm{Ph}}(\omega)$.

To determine the partial Huang-Rhys factors, and therefore the spectral
density, we treat the electron-phonon coupling from first principles, using
density functional theory (DFT) to calculate the normal modes of the ground-
and excited states of the hBN lattice with the considered defect complex.
Importantly, the normal modes of the ground- and excited states of the system
can be very different and can even be of different types. To properly account
for this, we employ the method outlined by Duschinsky [39] to calculate the
partial Huang-Rhys factors. For completeness we also checked the generating
function approach [40] and found almost identical results, meaning that the
modes in the ground and excited states for the system studied in the present
work are in fact very similar. The DFT calculations of the phonons were
performed for periodically repeated defect complexes in hBN monolayers using
the Vienna Ab Initio Simulation Package (VASP) [41]. In order to avoid
interactions between the periodic copies of the defects, we used a 9x9x1
supercell for the calculations, which we relax on a 3x3x1 K-point grid, to a
force convergence of $10^{-3}$ eV Å-1 using a plane wave cut-off of 700 eV.
All calculations were performed with the HSE06 exchange correlation
functional, which is essential for an accurate description of the electron-
phonon coupling [42]. The normal modes and the dynamical matrix was calculated
at the $\Gamma$-point. Note that we relax the defect in the in-plane
constrained C${}_{2}V$ configuration.

Figure 2: Time dependent emitter population (top) and cavity occupation
(bottom) for different coupling strengths and cavity widths. In contrast to
the effective width spectral density (dashed), the _ab initio_ spectral
density (solid orange) predicts significant oscillations in the cavity
occupation over all parameter regimes. Parameters used are $T=4$ K, $\Gamma=4$
meV, and $\hbar\Omega_{c}=\hbar\omega_{e}-\hbar\lambda$ where
$\lambda=\int_{0}^{\infty}d\omega\omega^{-2}\mathcal{J}_{\mathrm{Ph}}(\omega)$
is the reorganisation energy. The step size used to obtain convergence was
$dt=3.2$ fs, with SVD cut-off $\epsilon_{\mathrm{C}}=10^{-6}$.

The resulting spectral density is shown in Fig. 1(b), where we see multiple
sharp peaks present. These peaks correspond to high-$Q$ phonon modes present
in hBN; notably from the atomistic simulations, we can assign the dominant
contributions to the peak $\sim 125~{}$meV from de-localized defect breathing
modes, in which the atoms surrounding the defect oscillate along the dipole
direction of the defect.

To understand the influence of these modes will have on the cavity-defect
system, we can calculate the bath correlation function
$C(\tau)=\int_{0}^{\infty}d\omega\mathcal{J}_{\mathrm{Ph}}(\omega)\left(\cos(\omega\tau)\coth(\omega/2k_{\mathrm{B}}T)-i\sin(\omega\tau)\right)$,
shown in Fig. 1(c), where $T$ is the temperature, $k_{\mathrm{B}}$ is the
Boltzmann constant, and we have assumed the phonon bath to be initially in
thermal equilibrium. This correlation function encodes the timescales over
which memory effects last between the system and environment [43], for a full
derivation we refer the reader to the SI. The sharp peaks in the _ab intio_
spectral density lead to long lived oscillation in the correlation function,
which can be attributed to phonon modes with large $Q$-factors. To highlight
the importance of this structure we can replace the spectral density with a
broadened function with the same total Huang-Rhys parameter, as shown by the
dashed curve in Fig. 1(b). The corresponding correlation function, shown in
Fig. 1(c), has only small oscillations which are rapidly damped.

To study the time evolution of the composite cavity-emitter system, we employ
the time-evolving matrix product operator (TEMPO) algorithm initially
developed by Strathearn  _et al_ [28]. TEMPO is a powerful method for the
study of open quantum systems in strong coupling regimes [29, 30], and has
been applied to study quantum thermodynamics in the strong coupling regime
[44], as well as non-additive phenomena in quantum optics [31], and optimal
control [45]. The starting point for the TEMPO formalism, alongside other
numerically exact path integral methods [46, 47, 48], is the Trotter
decomposition [49], where for a sufficiently small time increment $\delta t$,
the propagator for the open system can be factorised such that
$\mathcal{U}_{\delta t}=e^{(\mathcal{L}_{S}+\mathcal{L}_{B})\delta
t}\approx\mathcal{V}_{\delta t}^{1/2}\mathcal{W}_{\delta t}\mathcal{V}_{\delta
t}^{1/2}+\mathcal{O}(\delta t^{3})$. The superoperator $\mathcal{V}_{\delta
t}$ is as defined above, and captures the evolution of the system and
dissipation through the external electromagnetic fields. The interaction
between the system and phonon environment is captured through
$\mathcal{W}_{\delta t}=\exp(\delta t\mathcal{L}_{\mathrm{B}})$, and is given
by
$\mathcal{L}_{\mathrm{B}}\rho=-i\left[H_{\mathrm{I}}^{\mathrm{Ph}}+H_{\mathrm{B}}^{\mathrm{Ph}},\rho\right]$.

This partitioning allows one to construct a discrete-time influence functional
of Feynman-Vernon type [46, 47], which captures the influence of the
environment to all orders in the interaction strength. The reduced state after
$k$-time-steps can be expressed as:

$\rho^{\alpha_{k}}=\sum_{\vec{\alpha},\vec{\beta}}\rho^{\alpha_{0}}\prod_{j=1}^{k}\mathcal{V}_{\beta_{j}}^{\alpha_{j}}\mathcal{V}^{\beta_{j}}_{\alpha_{j-1}}\mathcal{F}_{\beta_{k}\cdots\beta_{1}},$
(2)

where we have introduced the compound indices $\alpha_{k}=(s_{k},r_{k})$ and
$\beta_{k}=(t_{k},u_{k})$, yielding the density matrix elements
$\rho^{\alpha_{k}}=\bra{r_{k}}\rho\ket{s_{k}}$, and the system superoperator
is given by
$\mathcal{V}^{\alpha_{k}}_{\beta_{k}}=\bra{r_{k}}\mathcal{V}_{\delta
t}[\ket{t_{k}}\bra{u_{k}}]\ket{s_{k}}$. The influence of the phonon
environment is captured by the influence tensor
$\mathcal{F}_{\beta_{k}\cdots\beta_{1}}=\tr_{E}(\mathcal{W}^{\beta_{k}}\dots\mathcal{W}^{\beta_{1}}[\tau_{B}])$,
where $\mathcal{W}^{\beta_{k}}=\bra{t_{k}}\mathcal{W}_{\delta
t}[\ket{t_{k}}\bra{u_{k}}]\ket{u_{k}}$. By taking the initial state of the
environment to be in thermal equilibrium
$\tau_{B}=\exp(-\sum_{k}\nu_{k}b_{k}^{\dagger}b_{k}/k_{B}T)/\tr[\exp(-\sum_{k}\nu_{k}b_{k}^{\dagger}b_{k}/k_{B}T)]$,
the trace over the environmental degrees of freedom can be done analytically
[46, 47].

Crucially, the influence tensor scales exponentially in the number of time-
steps taken [46, 47]; while applying a finite time memory approximation can
reduce the computational cost of propagating the reduced state of the system
out to long times [46, 47, 50], it limits one to scenarios where key dynamics
occur on short timescales. A key insight of Strathearn  _et al_ [28] was that
the influence tensor may be represented in matrix product operator (MPO) form.
This allows one to encode the exponentially growing tensor as a tensor
network, and apply tensor compression [51] to reduce the rank of the elements
in the network, thereby circumventing exponential scaling. The reduced state
of the electronic system is then calculated by contracting the network down
after each time-step. The convergence of the TEMPO algorithm is sensitive to
taking a sufficiently small timestep $\delta t$, and the degree to which the
tensors are compressed. Further details of TEMPO and its convergence
properties are discussed in the SI.

We now consider the dynamics of the cavity-defect system initialised in the
state $\rho(0)=\ket{e,0}\\!\bra{e,0}$. Fig. 2 shows the time evolution of the
excited state population and the cavity mode occupation for different cavity
coupling parameters. We find that across all parameter regimes TEMPO predicts
complex oscillations in the cavity occupation. We attribute these oscillations
to high-$Q$ phonon modes in the environment, which lead to long-lived
oscillations in the bath correlation function shown in Fig. 1c, indeed these
oscillations follow closely those observed in the bath correlation function.
We confirm this by repeating the TEMPO calculations with the effective width
spectral density shown in Fig. 1b. Since this spectral density has the same
total Huang-Rhys parameter, naïvely we might expect similar dynamics to
emerge. However, the resulting dynamics in Fig. 2 show no oscillations at weak
light-matter coupling, suggesting it is indeed the structure in the spectral
density that has the dominant contribution to non-Markovian behaviour.
Interestingly, phonon induced oscillations are present in the cavity dynamics
across all parameter regimes studied, but only emerge in the emitter dynamics
in the strong light-matter coupling regime. This is a consequence of phonon
processes inducing fluctuations to the defects’ excited state energy; the
cavity is sensitive to these fluctuations through the dipole coupling with the
emitter, leading to the observed oscillations in the occupation. However,
these effects do become increasingly visible in the emitter dynamics with
increasing coupling strength, with the onset of vacuum Rabi oscillations.

We now turn our attention to the linear spectrum of the cavity-defect system.
The linear response of an electronic system coupled to a quantized mode can be
formulated in two ways, by an external field coupling to the electronic
degrees of freedom, or by an external current pumping the cavity mode [52,
53]. In both cases the driving can be included as a perturbation at the
Hamiltonian level, such that $H^{\prime}(t)=H+E(t)\cdot\mu$. Here the system
transition operator can take on two values
$\mu=\\{\sigma^{\dagger}+\sigma,a^{\dagger}+a\\}$. The first is the dipole
operator of the defect, describing the scattering of light directly off the
two level transition, and forms the basis of standard linear response theory
[54]. The second is the quadrature operator of the cavity field, which induces
a static polarisation of the cavity mode and excites real photons into the
field [53]. $E(t)$ then corresponds to the external field or the external
current respectively.

Figure 3: (a-c) the linear absorption spectra when probing the cavity mode,
for various values of the scaling parameter $\alpha_{\mathrm{HRF}}$. Here
structure in the absorption spectra becomes increasingly important for large
coupling, and broad cavity widths. (d) shows the change in peak position as a
function of $\alpha_{\mathrm{HRF}}$ (points) for the cavity parameters studied
in (a-c). All other parameters are the same as Fig. 2

Treating the driving as weak, we can extract the linear response spectrum
using density matrix perturbation theory [54]. Taking the semi-impulsive limit
such that $E(t)\propto e^{i\omega_{D}t}\delta(t)$, where $\omega_{D}$ is the
driving frequency, we obtain the absorption spectrum:

$A(\omega)=2\mathrm{Re}\left[\int_{0}^{\infty}dte^{i\omega
t}S^{(1)}(t)\right],$ (3)

where we have introduced the first-order response function
$S^{(1)}(t)=\tr(\mu(t)[\mu(0),\chi(-\infty)])$. The global equilibrium density
operator is given by
$\chi(-\infty)=\ket{g,0}\\!\bra{g,0}\otimes\tau_{B}\otimes\ket{\\{0\\}}\\!\bra{\\{0\\}}$,
where $\ket{\\{0\\}}$ denotes the vacuum state of the electromagnetic field.
The above expressions allow us to extend TEMPO to calculate the first-order
response function by propagating the initial system state
$\rho(0)=\mu(0)\ket{g,0}\\!\bra{g,0}$.

Figs. 3(a-c) show the absorption spectra of the cavity mode when driven by an
external current for various cavity parameters. To understand the role of
phonons in these spectra, we artificially include a scaling parameter to
spectral density
$\mathcal{J}_{\mathrm{Ph}}(\omega)\rightarrow\alpha_{\mathrm{HRF}}\mathcal{J}_{\mathrm{Ph}}(\omega)$,
where $\alpha_{\mathrm{HRF}}\in[0,1]$, equivalent to scaling the total Huang-
Rhys parameter. Fig. 3a shows absorption spectrum at the onset of the strong
coupling regime, where $g=\kappa$. In the absence of phonons
($\alpha_{\mathrm{HRF}}=0$) we see a clear Rabi splitting, with the separation
of the peaks determined by the the light-matter coupling $g$. As
$\alpha_{\mathrm{HRF}}$ increases, the separation of these peaks is reduced.
We can understand this reduction by appealing to the polaron formalism
commonly used to study the behaviour semi-conductor quantum dots [20, 55]:
here, the light-matter coupling strength is reduced by the Frank-Condon factor
$\mathscr{F}=\exp(-\alpha_{\mathrm{HRF}}S_{\mathrm{Tot}}/2)$, which accounts
for the geometric difference of the phonon modes associated to the emitter
ground and excited states. Fig. 3d compares the peak separation as a function
of the scaling parameter $\alpha_{\mathrm{HRF}}$ for the three parameter
regimes. For the parameters associated with Fig. 3a, the peak separation
follows the Frank-Condon factor $\mathscr{F}$ (solid curve).

In regimes of stronger light-matter coupling show in Fig. 3b and c, a more
complex picture emerges: significant structure becomes apparent in the spectra
as $\alpha_{\mathrm{HRF}}$ increases. Of particular interest is the departure
from the well understood polaronic physics seen in Fig. 3a, which is
highlighted in Fig. 3d. Here we see that at stronger coupling regimes ($g=50$
meV), the renormalisation no longer follows the Frank-Condon factor, and in
regime of strong light-matter coupling ($g=100$ meV), we in-fact see the
splitting increase with $\alpha_{\mathrm{HRF}}$. We attribute this behaviour
to a hybridisation of the light-matter polariton and high-$Q$ phonon modes.
The resultant state is a tri-partite quasi-particle with characteristics of
light, matter, and vibrations. This can be seen most clearly in the upper-
polariton of Fig. 3c, where an additional splitting emerges when
$\alpha_{\mathrm{HRF}}\rightarrow 1$. This interpretation is supported by the
C2CN calculations shown in the SI; this defect complex has a spectral density
with little low energy structure, such that at $g=50$ meV no hybridisation
occurs, and the renormalisation of the Rabi frequency follows closely the
Frank-Condon factor $\mathscr{F}$. It is only at higher light-matter coupling
strengths ($g=100$ meV) for C2CN, when the polariton splitting approaches
resonance with a high-Q phonon mode that we observe a departure from Frank
Condon physics, heralding hybridised polariton-polaron states.

_Conclusion—_ In this letter we have combined atomistic simulations of a
defect complex in hBN with TEMPO, a numerically exact and fully quantum
mechanical simulation method. Our hybrid approach allows us to study realistic
emitters beyond phenomenological and approximate treatments [22], providing a
complete description of electron-phonon interactions in condensed matter
single photon emitters in optical or plasmonic cavities, with direct relevance
to on-going experiments. Furthermore, by considering the cavity quantum
electrodynamics of a defect across the weak and strong light-matter coupling
regimes, we have shown that strong coupling to high-Q vibrational modes play a
significant role in determining the dynamics of the cavity-defect system, even
in the weak light-matter coupling regime. At strong light-matter coupling, the
absorption spectrum show clear signatures of hybridisation between the light-
matter polaritons and phonon modes inherent to hBN.

The method we present here is general, and not restricted to the material
system or specific defect complexes studied here, with potential application
to organic polaritons [56]. It is worth noting that we have restricted
ourselves to low temperatures in the above discussion; at elevated
temperatures phonon processes beyond linear electron-phonon coupling become
important through mechanisms such as the Jahn-Teller effect [57].

_Acknowledgments—_ We thank Mortiz Fischer and Emil V. Denning for insightful
conversations. The work is supported by the BIOMAG project funded by the Novo
Nordisk Foundation (grant no. NNF21OC0066526) We also acknowledge funding from
the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program Grant No. 773122 (LIMA). K. S. T. is a Villum
Investigator supported by Villum Foundation (grant no. 37789). This work was
supported by the Danish National Research Foundation through NanoPhoton -
Center for Nanophotonics, Grant No. DNRF147 and Center for Nanostructured
Graphene, Grant No. DNRF103. NS acknowledges support from the Villum
Foundation through grant no. 00028233 and from the Independent Research Fund
Denmark - Natural Sciences (Project no. 0135-00403B). J.I.-S. acknowledges
support from the Royal Commission for the Exhibition of 1851.

## References

  * Chakraborty  _et al._ [2015] C. Chakraborty, L. Kinnischtzke, K. M. Goodfellow, R. Beams, and A. N. Vamivakas, Voltage-controlled quantum light from an atomically thin semiconductor, Nature Nanotechnology 10, 507 (2015).
  * He  _et al._ [2015] Y.-M. He, G. Clark, J. R. Schaibley, Y. He, M.-C. Chen, Y.-J. Wei, X. Ding, Q. Zhang, W. Yao, X. Xu, C.-Y. Lu, and J.-W. Pan, Single quantum emitters in monolayer semiconductors, Nature Nanotechnology 10, 497 (2015).
  * Koperski  _et al._ [2015] M. Koperski, K. Nogajewski, A. Arora, V. Cherkez, P. Mallet, J.-Y. Veuillen, J. Marcus, P. Kossacki, and M. Potemski, Single photon emitters in exfoliated wse2 structures, Nature Nanotechnology 10, 503 (2015).
  * Srivastava  _et al._ [2015] A. Srivastava, M. Sidler, A. V. Allain, D. S. Lembke, A. Kis, and A. Imamoğlu, Optically active quantum dots in monolayer $\mathrm{WSe}_{2}$, Nature Nanotechnology 10, 491 (2015).
  * Klein  _et al._ [2019] J. Klein, M. Lorke, M. Florian, F. Sigger, L. Sigl, S. Rey, J. Wierzbowski, J. Cerne, K. Müller, E. Mitterreiter, _et al._ , Site-selectively generated photon emitters in monolayer mos2 via local helium ion irradiation, Nature communications 10, 1 (2019).
  * Michaelis de Vasconcellos  _et al._ [2022] S. Michaelis de Vasconcellos, D. Wigger, U. Wurstbauer, A. W. Holleitner, R. Bratschitsch, and T. Kuhn, Single-photon emitters in layered van der waals materials, physica status solidi (b) , 2100566 (2022).
  * Kianinia  _et al._ [2022] M. Kianinia, Z.-Q. Xu, M. Toth, and I. Aharonovich, Quantum emitters in 2d materials: Emitter engineering, photophysics, and integration in photonic nanostructures, Applied Physics Reviews 9, 011306 (2022).
  * Tran  _et al._ [2016] T. T. Tran, K. Bray, M. J. Ford, M. Toth, and I. Aharonovich, Quantum emission from hexagonal boron nitride monolayers, Nature Nanotechnology 11, 37 (2016).
  * Aharonovich and Toth [2017] I. Aharonovich and M. Toth, Quantum emitters in two dimensions, Science 358, 170 (2017).
  * Fischer  _et al._ [2021] M. Fischer, J. Caridad, A. Sajid, S. Ghaderzadeh, M. Ghorbani-Asl, L. Gammelgaard, P. Bøggild, K. S. Thygesen, A. Krasheninnikov, S. Xiao, _et al._ , Controlled generation of luminescent centers in hexagonal boron nitride by irradiation engineering, Science Advances 7, eabe7138 (2021).
  * 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, _et al._ , Identifying carbon as the source of visible single-photon emission from hexagonal boron nitride, Nature materials 20, 321 (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, _et al._ , Initialization and read-out of intrinsic spin defects in a van der waals crystal at room temperature, Nature materials 19, 540 (2020).
  * Hoese  _et al._ [2020] M. Hoese, P. Reddy, A. Dietrich, M. K. Koch, K. G. Fehler, M. W. Doherty, and A. Kubanek, Mechanical decoupling of quantum emitters in hexagonal boron nitride from low-energy phonon modes, Science advances 6, eaba6038 (2020).
  * Cassabois  _et al._ [2016] G. Cassabois, P. Valvin, and B. Gil, Hexagonal boron nitride is an indirect bandgap semiconductor, Nature photonics 10, 262 (2016).
  * Tran  _et al._ [2019] T. T. Tran, C. Bradac, A. S. Solntsev, M. Toth, and I. Aharonovich, Suppression of spectral diffusion by anti-stokes excitation of quantum emitters in hexagonal boron nitride, Applied Physics Letters 115, 071102 (2019).
  * Li  _et al._ [2021] K. Li, T. Smart, and Y. Ping, $\mathrm{C}_{2}\mathrm{C}_{\mathrm{N}}$ as a 2 ev single-photon emitter candidate in hexagonal boron nitride, arXiv preprint arXiv:2110.01787 (2021).
  * Vogl  _et al._ [2019] T. Vogl, R. Lecamwasam, B. C. Buchler, Y. Lu, and P. K. Lam, Compact cavity-enhanced single-photon generation with hexagonal boron nitride, ACS Photonics 6, 1955 (2019).
  * Häußler  _et al._ [2021] S. Häußler, G. Bayer, R. Waltrich, N. Mendelson, C. Li, D. Hunger, I. Aharonovich, and A. Kubanek, Tunable fiber-cavity enhanced photon emission from defect centers in hbn, Advanced Optical Materials 9, 2002218 (2021).
  * Denning  _et al._ [2020] E. V. Denning, J. Iles-Smith, N. Gregersen, and J. Mork, Phonon effects in quantum dot single-photon sources, Opt. Mater. Express 10, 222 (2020).
  * Nazir and McCutcheon [2016] A. Nazir and D. P. S. McCutcheon, Modelling exciton–phonon interactions in optically driven quantum dots, Journal of Physics: Condensed Matter 28, 103002 (2016).
  * Vuong  _et al._ [2016] T. Vuong, G. Cassabois, P. Valvin, A. Ouerghi, Y. Chassagneux, C. Voisin, and B. Gil, Phonon-photon mapping in a color center in hexagonal boron nitride, Physical review letters 117, 097402 (2016).
  * Feldman  _et al._ [2019] M. A. Feldman, A. Puretzky, L. Lindsay, E. Tucker, D. P. Briggs, P. G. Evans, R. F. Haglund, and B. J. Lawrie, Phonon-induced multicolor correlations in hbn single-photon emitters, Phys. Rev. B 99, 020101 (2019).
  * Khatri  _et al._ [2019] P. Khatri, I. J. Luxmoore, and A. J. Ramsay, Phonon sidebands of color centers in hexagonal boron nitride, Phys. Rev. B 100, 125305 (2019).
  * Grosso  _et al._ [2020] G. Grosso, H. Moon, C. J. Ciccarino, J. Flick, N. Mendelson, L. Mennel, M. Toth, I. Aharonovich, P. Narang, and D. R. Englund, Low-temperature electron–phonon interaction of quantum emitters in hexagonal boron nitride, ACS Photonics 7, 1410 (2020).
  * Sajid and Thygesen [2020] A. Sajid and K. S. Thygesen, Vncb defect as source of single photon emission from hexagonal boron nitride, 2D Materials 7, 031007 (2020).
  * Jara  _et al._ [2021] C. Jara, T. Rauch, S. Botti, M. A. Marques, A. Norambuena, R. Coto, J. Castellanos-Águila, J. R. Maze, and F. Munoz, First-principles identification of single photon emitters based on carbon clusters in hexagonal boron nitride, The Journal of Physical Chemistry A 125, 1325 (2021).
  * Li  _et al._ [2022] K. Li, T. J. Smart, and Y. Ping, Carbon trimer as a 2 ev single-photon emitter candidate in hexagonal boron nitride: A first-principles study, Physical Review Materials 6, L042201 (2022).
  * Strathearn  _et al._ [2018] A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W. Lovett, Efficient non-markovian quantum dynamics using time-evolving matrix product operators, Nature communications 9, 1 (2018).
  * Jørgensen and Pollock [2019] M. R. Jørgensen and F. A. Pollock, Exploiting the causal tensor network structure of quantum processes to efficiently simulate non-markovian path integrals, Phys. Rev. Lett. 123, 240602 (2019).
  * Gribben  _et al._ [2020] D. Gribben, A. Strathearn, J. Iles-Smith, D. Kilda, A. Nazir, B. W. Lovett, and P. Kirton, Exact quantum dynamics in structured environments, Phys. Rev. Research 2, 013265 (2020).
  * Gribben  _et al._ [2022] D. Gribben, D. M. Rouse, J. Iles-Smith, A. Strathearn, H. Maguire, P. Kirton, A. Nazir, E. M. Gauger, and B. W. Lovett, Exact dynamics of nonadditive environments in non-markovian open quantum systems, PRX Quantum 3, 010321 (2022).
  * Svendsen  _et al._ [2021] M. K. Svendsen, Y. Kurman, _et al._ , Combining density functional theory with macroscopic qed for quantum light-matter interactions in 2$\mathrm{D}$ materials, Nature communications 12, 1 (2021).
  * Latini  _et al._ [2021] S. Latini, U. De Giovannini, E. J. Sie, N. Gedik, H. Hübener, and A. Rubio, Phonoritons as hybridized exciton-photon-phonon excitations in a monolayer $\mathrm{h-BN}$ optical cavity, Phys. Rev. Lett. 126, 227401 (2021).
  * Sajid  _et al._ [2018] A. Sajid, J. R. Reimers, and M. J. Ford, Defect states in hexagonal boron nitride: Assignments of observed properties and prediction of properties relevant to quantum computation, Phys. Rev. B 97, 064101 (2018).
  * Note [1] This analysis is repeated in the supplementary information for the C2CN defect.
  * Carmichael [1999] H. J. Carmichael,  _Statistical methods in quantum optics 1: master equations and Fokker-Planck equations_ , Vol. 1 (Springer Science & Business Media, 1999).
  * Maguire  _et al._ [2019] H. Maguire, J. Iles-Smith, and A. Nazir, Environmental nonadditivity and franck-condon physics in nonequilibrium quantum systems, Physical review letters 123, 093601 (2019).
  * Mahan [2013] G. D. Mahan,  _Many-particle physics_ (Springer Science & Business Media, 2013).
  * Borrelli  _et al._ [2013] R. Borrelli, A. Capobianco, and A. Peluso, Franck–condon factors—computational approaches and recent developments, Canadian Journal of Chemistry 91, 495 (2013).
  * Alkauskas  _et al._ [2014] A. Alkauskas, B. B. Buckley, D. D. Awschalom, and C. G. Van de Walle, First-principles theory of the luminescence lineshape for the triplet transition in diamond nv centres, New Journal of Physics 16, 073026 (2014).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Physical review B 47, 558 (1993).
  * Heyd and Scuseria [2004] J. Heyd and G. E. Scuseria, Efficient hybrid density functional calculations in solids: Assessment of the heyd–scuseria–ernzerhof screened coulomb hybrid functional, The Journal of chemical physics 121, 1187 (2004).
  * Breuer  _et al._ [2002] H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Popovic  _et al._ [2021] M. Popovic, M. T. Mitchison, A. Strathearn, B. W. Lovett, J. Goold, and P. R. Eastham, Quantum heat statistics with time-evolving matrix product operators, PRX Quantum 2, 020338 (2021).
  * Fux  _et al._ [2021] G. E. Fux, E. P. Butler, P. R. Eastham, B. W. Lovett, and J. Keeling, Efficient exploration of hamiltonian parameter space for optimal control of non-markovian open quantum systems, Phys. Rev. Lett. 126, 200401 (2021).
  * Makri and Makarov [1995a] N. Makri and D. E. Makarov, Tensor propagator for iterative quantum time evolution of reduced density matrices. i. theory, The Journal of chemical physics 102, 4600 (1995a).
  * Makri and Makarov [1995b] N. Makri and D. E. Makarov, Tensor propagator for iterative quantum time evolution of reduced density matrices. ii. numerical methodology, The Journal of chemical physics 102, 4611 (1995b).
  * Cygorek  _et al._ [2022] M. Cygorek, M. Cosacchi, A. Vagov, V. Axt, B. Lovett, J. Keeling, and E. Gauger, Simulation of open quantum systems by automated compression of arbitrary environments, Nature Physics (2022).
  * Trotter [1959] H. F. Trotter, On the product of semi-groups of operators, Proceedings of the American Mathematical Society 10, 545 (1959).
  * Strathearn  _et al._ [2017] A. Strathearn, B. W. Lovett, and P. Kirton, Efficient real-time path integrals for non-markovian spin-boson models, New Journal of Physics 19, 093009 (2017).
  * Orús [2014] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Annals of Physics 349, 117 (2014).
  * Ruggenthaler  _et al._ [2018] M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Appel, and A. Rubio, From a quantum-electrodynamical light–matter description to novel spectroscopies, Nature Reviews Chemistry 2, 1 (2018).
  * Flick  _et al._ [2019] J. Flick, D. M. Welakuh, M. Ruggenthaler, H. Appel, and A. Rubio, Light–matter response in nonrelativistic quantum electrodynamics, ACS photonics 6, 2757 (2019).
  * Mukamel [1999] S. Mukamel,  _Principles of nonlinear optical spectroscopy_ , 6 (Oxford University Press on Demand, 1999).
  * Iles-Smith  _et al._ [2017] J. Iles-Smith, D. P. McCutcheon, A. Nazir, and J. Mørk, Phonon scattering inhibits simultaneous near-unity efficiency and indistinguishability in semiconductor single-photon sources, Nature Photonics 11, 521 (2017).
  * Del Pino  _et al._ [2018] J. Del Pino, F. A. Schröder, A. W. Chin, J. Feist, and F. J. Garcia-Vidal, Tensor network simulation of non-markovian dynamics in organic polaritons, Physical review letters 121, 227401 (2018).
  * White  _et al._ [2021] S. White, C. Stewart, A. S. Solntsev, C. Li, M. Toth, M. Kianinia, and I. Aharonovich, Phonon dephasing and spectral diffusion of quantum emitters in hexagonal boron nitride, Optica 8, 1153 (2021).