#  $T$ centres in photonic silicon-on-insulator material

E. R. MacQuarrie  Department of Physics, Simon Fraser University, Burnaby,
British Columbia, Canada  Photonic Inc., Vancouver, British Columbia, Canada
C. Chartrand  Department of Physics, Simon Fraser University, Burnaby, British
Columbia, Canada  D. B. Higginbottom  Department of Physics, Simon Fraser
University, Burnaby, British Columbia, Canada  K. J. Morse  Department of
Physics, Simon Fraser University, Burnaby, British Columbia, Canada  Photonic
Inc., Vancouver, British Columbia, Canada  V. A. Karasyuk  Department of
Physics, Simon Fraser University, Burnaby, British Columbia, Canada  S. Roorda
Department of Physics, Université de Montréal, Montreal, Quebec, Canada  S.
Simmons  Department of Physics, Simon Fraser University, Burnaby, British
Columbia, Canada  Photonic Inc., Vancouver, British Columbia, Canada
<EMAIL_ADDRESS>
###### Abstract

Global quantum networks will benefit from the reliable creation and control of
high-performance solid-state telecom photon-spin interfaces. T radiation
damage centres in silicon provide a promising photon-spin interface due to
their narrow $O$-band optical transition near $1326$ nm and long-lived
electron and nuclear spin lifetimes. To date, these defect centres have only
been studied as ensembles in bulk silicon. Here, we demonstrate the reliable
creation of high concentration T centre ensembles in the 220 nm device layer
of silicon-on-insulator (SOI) wafers by ion implantation and subsequent
annealing. We then develop a method that uses spin-dependent optical
transitions to benchmark the characteristic optical spectral diffusion within
these T centre ensembles. Using this new technique, we show that with minimal
optimization to the fabrication process high densities of implanted T centres
localized $\lesssim 100$ nm from an interface display $\sim 1$ GHz
characteristic levels of total spectral diffusion.

##  I Introduction

Point defect colour centres in solid state systems provide a promising
platform for emerging quantum technologies. A number of colour centres are
currently in development Awschalom et al. (2018); Zhang et al. (2020), but few
of these leading candidates are both natively integrated into silicon and
intrinsically operate at telecom wavelengths. Radiation damage centres in
silicon have recently emerged as promising light-matter interfaces that
simultaneously meet both of these criteria Buckley et al. (2017); Chartrand et
al. (2018); Beaufils et al. (2018); Redjem et al. (2020). Within this family
of defects, the T centre notably possesses highly coherent electron and
nuclear spin degrees of freedom and narrow, spin-dependent ensemble optical
transitions near $1326$ nm in the telecommunications $O$-band Bergeron et al.
(2020). These properties make the T centre a competitive candidate for
integration into future commercial-scale quantum networks with long-lived
quantum memory and computing capabilities.

Incorporating T centres into such quantum networks requires a reliable means
of making the centres in device-ready material such as silicon-on-insulator
(SOI). Other photon-spin interfaces have been produced by ion implantation and
thermal treatment Meijer et al. (2005); Chu et al. (2014); Falk et al. (2013);
Beaufils et al. (2018); Phenicie et al. (2019); Buckley et al. (2020); Redjem
et al. (2020). Here, we use similar techniques to generate high concentrations
of T centres in the 220nm device layer of commercial SOI wafers. We achieve T
centre concentrations $>30$ times higher than in previously measured bulk
samples and, in follow-on work kur , present evidence that these SOI
concentrations are $>1.7\times 10^{13}$ cm-3. Moreover, the inhomogeneous
broadening of the zero phonon line (ZPL) within implanted ensembles is only
$\sim 35\%$ larger than in T centre ensembles generated by electron
irradiating and annealing comparable bulk material.

The fabrication process used to generate T centres should preserve the
excellent spin and optical properties that T centre ensembles exhibit in bulk
silicon where the majority of centres are located far from material
interfaces. The increase in electric and magnetic field noise near those
interfaces can dramatically degrade a centre’s spin and optical properties.
Encouragingly, donor spin defects in silicon have already demonstrated
competitive electron and nuclear spin coherence times ($>0.5$ s, $>30$ s) at
distances of just $10$ nm from ${}^{28}\text{Si}$/SiO2 interfaces Muhonen et
al. (2014), and T centres can be expected to demonstrate comparable or better
spin properties at our target depth of $110$ nm. In the optical domain,
spectral diffusion — the time-varying changes in a centre’s optical transition
energies — often dominates the optical properties of emitters located near
interfaces. Spectral diffusion can reduce the indistinguishability of emitted
photons, limiting an emitter’s potential in quantum technologies. Measuring
the characteristic spectral diffusion of the generated T centre ensembles thus
provides an important metric for benchmarking our fabrication process. To this
end, we develop a method that uses spin-dependent optical transitions to
efficiently measure the characteristic long-term (total) spectral diffusion
within an ensemble of emitters.

Applying this new technique to our high-concentration T centre ensembles in
SOI, we measure characteristic total spectral diffusion of $\sim 1$ GHz.
Although this value is $\sim 10^{4}$ times larger than the $169$ kHz lifetime-
limited linewidth, excellent photon indistinguishability could still be
achieved by employing fabrication or measurement techniques designed to
stabilize the electromagnetic environment de Oliveira et al. (2017); Anderson
et al. (2019); Sangtawesin et al. (2019); Kasperczyk et al. (2020), by
filtering the emitted photons Bernien et al. (2013); Guo et al. (2019), or by
Purcell-enhancing the radiative decay rate Grange et al. (2015); Giesz et al.
(2015) as the T centre is expected to have a high radiative efficiency
Bergeron et al. (2020).

##  II T Centre Implantation

As shown in Fig. 1a, the T centre is thought to consist of two carbon atoms
sharing the substitutional site of one silicon atom; an interstitial hydrogen
terminates one of these carbon atoms while the other remains stable with an
unpaired electron Safonov et al. (1996). The current formation model for this
defect consists of an interstitial carbon capturing a hydrogen atom and then
migrating to a substitutional carbon site during heat treatment at
temperatures between $350^{\circ}$ and $600^{\circ}$C Safonov et al. (1996);
Bergeron et al. (2020). Starting from SOI with a $220$ nm thick Czochralski
(CZ) natural Si (${}^{\text{nat}}$Si) device layer and $3$ µm of buried SiO2,
we realize this formation process by implanting carbon and then hydrogen in a
fixed [C:H]=1:1 dose ratio, annealing between the two implants and then again
after the hydrogen implant.

We performed two fabrication runs: one with a 12C implant and doses of
$7\times 10^{12}$ cm-2 and a second with 13C and doses of $7\times 10^{13}$
cm-2. These doses, which were chosen from measurements detailed in the
Supplementary Materials (SM) SI , place the initial carbon concentrations well
above the carbon solubility limit for silicon at $1000^{\circ}$C Pichler
(2004). Carbon implants were performed at an energy of $38$ keV, while
hydrogen was implanted at $9$ keV to target overlapping implantation profiles
with a mean depth of $110$ nm.

In between the two implants, we perform a rapid thermal anneal at
$1000^{\circ}$C for $20$ seconds in an argon background to repair lattice
damage generated during the carbon implant and to incorporate carbon onto
substitutional lattice sites Berhanuddin et al. (2012); Beaufils et al.
(2018). The subsequent hydrogen implant then generates lattice damage and
interstitial carbon to promote T centre formation, but the smaller mass of the
hydrogen ions generates less total damage than the carbon implant SI . After
the hydrogen implant, we boil the samples for $1$ hour in deionized water to
further increase the hydrogen concentration before performing a final rapid
thermal anneal for $3$ min at a temperature ranging from $345^{\circ}$ to
$480^{\circ}$C in a nitrogen atmosphere.

Figure 1: T centre implantation development. (a) Atomic structure of the T
centre in Si as proposed in Ref. Safonov et al. (1996). (b) Integrated T
centre zero phonon line (ZPL) area in SOI as a function of the final anneal
temperature. (c,d) T centre photoluminescence (PL) at a temperature of $4.2$ K
for bulk 28Si and SOI implanted at fluences of $7\times 10^{12}$ cm-2 and
$7\times 10^{13}$ cm-2. The measured signal has been normalized to the
intensity of the free exciton (FE) recombination line. Panel (d) highlights
the collection window for our photoluminescence excitation (PLE) spectroscopy
measurements. The sharp dips within this window are water vapor absorption
lines. (e) ZPL PL for T centres in bulk float-zone ${}^{\text{nat}}$Si that
were created either by implantation or electron irradiation. The signal has
been normalized to the fitted T centre ZPL amplitude.

Fig. 1b shows the integrated area of the T centre ZPL as a function of the
final anneal temperature for a series of SOI samples. For these samples, we
find optimum temperatures of $400^{\circ}$C for a fluence of $7\times 10^{12}$
cm-2 and $420^{\circ}$C for $7\times 10^{13}$ cm-2. Unless otherwise noted, we
have normalized photoluminescence (PL) spectra to the free exciton (FE)
recombination line at $\sim 1097$ meV. Although the FE recombination line can
vary in intensity across samples as, for instance, the density of
recombination sites changes, it provides a means of roughly quantifying the T
centre concentration to enable comparison amongst samples.

In Fig. 1c,d, we plot the PL spectra under above-band ($532$ nm) excitation
for these particular SOI samples alongside that of an isotopically-purified
28Si bulk sample with an initial $1.5\times 10^{15}$ cm-3 carbon
concentration. T centres were created in this bulk sample by electron
irradiating and annealing (details in Ref. Bergeron et al. (2020)). Fig. 1d
provides a detailed view of the T centre sideband spectra, showing an isotopic
shift in the local vibrational modes (LVMs) as the $7\times 10^{13}$ cm-2 SOI
sample was implanted with 13C isotopes, the $7\times 10^{12}$ cm-2 SOI sample
was implanted with 12C, and the bulk 28Si sample contains a natural
distribution of carbon isotopes ($98.9\%$ 12C). In Fig. 1d, we have also
labeled both ZPLs within the T centre exciton doublet, where $\text{TX}_{1}$
sits $1.76$ meV above $\text{TX}_{0}$ Bergeron et al. (2020).

By normalizing the integrated area of the T centre ZPL to the nominal sample
thickness, we estimate that the implanted T centre densities for the $7\times
10^{12}$ cm-2 and $7\times 10^{13}$ cm-2 SOI samples are $>5$ and $>30$ times
larger than the T centre density in the bulk 28Si sample, respectively. Here,
we have used $220$ nm as the SOI thickness and the $2.5$ µm absorption depth
of $532$ nm light in silicon at cryogenic temperatures Dash and Newman (1955);
Jellison Jr. and Joshi (2018) as the nominal 28Si thickness. As the free
excitons generated by the $532$ nm excitation are not confined to the
absorption region in our bulk sample, this latter thickness serves as a lower
bound, making our relative concentration estimates lower bounds on the true
ratios. In the SM, we verify that the implanted T centres are localized to the
SOI device layer by selectively etching away the device layer and also by
using ultraviolet PL excitation, which has a much shallower absorption depth
than $532$ nm SI . We also note that this estimate of the T centre density
neglects differences in the defects’ local optical environment that could
modify their emission properties Barnes (1998).

In Fig. 1e, we compare the ZPL PL for T centres generated by implantation to
that of T centres created by electron irradiation. Both samples in this figure
are bulk, float-zone (FZ) ${}^{\text{nat}}$Si. One was implanted with 13C and
hydrogen at a fluence of $7\times 10^{12}$ cm-2, while the other was electron
irradiated and annealed to convert native carbon and hydrogen into T centres
Bergeron et al. (2020). We once again see the expected 13C isotopic shift in
the implanted ZPL, and at $4.2$ K, we report ZPL linewidths of $6.67\pm 0.09$
GHz in the irradiated FZ sample and $9.01\pm 0.08$ GHz for the implanted FZ
wafer. As will be shown below, these linewidths lie well above the linewidth
of our 28Si sample and below that of our CZ SOI samples. Among these
measurements, however, this comparison between FZ ${}^{\text{nat}}$Si samples
provides the best comparison between implanted and irradiated linewidths as
each used similar material outside the “semiconductor vacuum” of 28Si Saeedi
et al. (2013). We thus report only a $\sim 35\%$ increase in the inhomogeneous
broadening for our implanted ensembles compared to ensembles generated by
electron irradiation, and we expect that the broader lines we observe in CZ
SOI will be reduced for T centres created by implantation into a FZ SOI device
layer.

##  III Hyperpolarization

In other solid state emitters that lack inversion symmetry, spectral diffusion
has been shown to increase when centres are generated by ion implantation van
Dam et al. (2019); Kasperczyk et al. (2020) or as centres approach interfaces
due to a dramatic increase in electric and magnetic field noise in these
regions Ishikawa et al. (2012); Evans et al. (2016); Crook et al. (2020).
Typically, spectral diffusion is measured by tracking spectral lines of
individual emitters. Here, we develop a novel technique for measuring spectral
diffusion within an ensemble of emitters and use it to compare the optical
properties of the T centres in our SOI samples to those in bulk samples.

Under spin-selective resonant optical driving, the T centre ground state
electron spin will hyperpolarize through a process shown schematically in Fig.
2a. This unpaired electron provides a spin degree of freedom with an isotropic
$g$-factor of $g_{e}=2.005$ in the orbital ground state. Upon resonant
generation of a bound exciton (BE), the ground state electron pairs with the
BE electron to form a spin-$0$ singlet. The BE hole then provides an unpaired
spin with an anisotropic $g$-factor $g_{h}$ in the BE $\text{TX}_{0}$ state. A
magnetic field $B_{z}$ thus generates a spin splitting $g_{e}\mu_{B}B_{z}$
($g_{h}\mu_{B}B_{z}$) in the ground (excited, or $\text{TX}_{0}$) state where
$\mu_{B}$ is the Bohr magneton. Because these splittings can differ in
magnitude, for Zeeman splittings $\varepsilon_{B}(B_{z})=(\pm g_{e}\pm
g_{h})\mu_{B}B_{z}$ larger than the homogeneous optical linewidth, only one
spin sublevel will be excited by the resonant optical drive. Spin non-
conserving relaxation will then shelve the system in whichever state is not
addressed by the optical drive, hyperpolarizing the electron spin as the T
centre is optically cycled. For the T centre, this hyperpolarization mechanism
is efficient Bergeron et al. (2020).

Figure 2: Hyperpolarization in T centres (a) Hyperpolarization mechanism for T
centres. (b) Depiction of how spectral diffusion can limit T centre
hyperpolarization. (c) PLE measurements of T centres in 28Si performed as a
function of magnetic field at a temperature of 4.2 K. Detuning is measured
from the ZPL energy of $935.062$ meV. The dashed lines indicate the maximum
and minimum Zeeman splittings expected for this magnetic field and sample
orientation. The insets show a Lorentzian fit to the residual-field PLE
spectrum and a fit of the $\Delta=0$ µeV linecut to the hyperpolarization
model developed in the text. (d) Levels and transitions included in the rate
model for spin-dependent optical transitions. $L_{\text{B(C)}}$ represents
optical excitation in the spin-down (spin-up) manifold, and the green arrows
show fluorescent relaxation pathways. (e) Schematic depiction of the loss in
PLE signal as a magnetic field splits the T centre spin sublevels. The
resulting PL spectrum (green curve) is described by Eq. 1.

Using this process, we can reveal the characteristic homogeneous optical
linewidth of an ensemble by measuring how the resonantly-driven ZPL evolves in
a magnetic field. In the absence of any homogeneous broadening effects, an
individual T centre’s optical linewidth will approach its lifetime-limited
lower bound and observable hyperpolarization at low magnetic fields is
expected. This onset of hyperpolarization only weakly depends on the
inhomogeneously-broadened ensemble lineshape and will apply to ensembles in
28Si and ${}^{\text{nat}}$Si alike. Homogeneous effects such as thermal
broadening or spectral diffusion will delay this onset of hyperpolarization to
larger magnetic fields as shown schematically in Fig. 2b. If these homogeneous
mechanisms have a characteristic broadening on the same scale as the spin
splitting, they will ensure that each spin sublevel is intermittently
addressed by the optical field, reducing the level of hyperpolarization in the
system and increasing the measured PL.

In Fig. 2c, we perform PLE spectroscopy of the T centre ZPL in the bulk 28Si
sample by scanning a tunable laser over the optical resonance and measuring
the PL emission within the collection window highlighted in Fig. 1d. Repeating
this measurement as a function of $B_{z}$, we see that, rather than Zeeman
splitting along the dashed lines of Fig. 2c, the PLE trace remains stationary
and slowly decays as the interrogated T centres hyperpolarize. This
measurement was performed at a temperature of $4.2$ K where the ZPL linewidth
is thermally-broadened Bergeron et al. (2020) and with $B_{z}$ aligned along
the $\langle 100\rangle$ crystal axis. For this field orientation, we expect
the hole $g$-factors for the twelve possible T centre orientational subsets
relative to $B_{\langle 100\rangle}$ to take on the values
$g_{h}=\\{0.91,2.55\\}$ with multiplicities of $4$ and $8$, respectively
Safonov et al. (1995); SI .

The measured hyperpolarization dynamics can be modeled by noting that both
spin sublevels must be driven to prevent population shelving. With a magnetic
field applied, the spin sublevels split by $\varepsilon_{B}(B_{z})$ where each
sublevel is assumed to have the same optical homogeneous linewidth $\Gamma$.
As the spin states separate and the probability of driving both sublevels
decreases, the system starts to hyperpolarize. Assuming weak optical driving,
we model this process with a rate model using the levels and transitions shown
schematically in Fig. 2d and find that the amplitude of a partially
hyperpolarized PLE spectrum (Fig. 2e) is proportional to

$A(B_{z},\Delta)\propto\frac{L_{\text{B}}(B_{z},\Delta)L_{\text{C}}(B_{z},\Delta)}{L_{\text{B}}(B_{z},\Delta)+L_{\text{C}}(B_{z},\Delta)}$
(1)

where $L_{\text{B}(\text{C})}(B_{z},\Delta)$ is the Lorentzian amplitude of
the spin-down (spin-up) sublevel and $\Delta$ is the optical detuning from the
zero-field resonance. For simplicity, we have neglected optical transitions
between states with different spin orientations, which have smaller spectral
overlap with other transitions and have been shown for T centres to be much
weaker than optical transitions between spin states of the same orientation
Bergeron et al. (2020). The effect of these transitions is examined in the SM
SI .

Evaluating Eq. 1 and normalizing to the $B_{z}=0$, $\Delta=0$ value, we obtain
an expression for the PLE amplitude as a function of magnetic field for a
single T centre:

$A_{\text{h}}(B_{z},\Delta)=\frac{\Gamma^{2}}{4\Delta^{2}+\Gamma^{2}+\varepsilon_{B}(B_{z})^{2}}.$
(2)

Because we simultaneously interrogate a large number of T centres in twelve
possible orientational subsets, we model our measurement by averaging Eq. 2
across the twelve values of $g_{h}$. Here, we have assumed that each
orientation contributes equally to the measured signal.

Fitting the data in Fig. 2c to this model, we extract a homogeneous linewidth
of $\Gamma=$ $250\pm 100$ MHz, which agrees very well with the purely thermal
linewidth of $230$ MHz at $4.2$ K calculated from the thermal broadening data
in Ref. Bergeron et al. (2020) as well as the
$\Lambda_{\text{PLE}}=~{}$$329.8\pm 0.6$ MHz linewidth obtained by fitting the
PLE spectrum in the $B_{\langle 100\rangle}=36$ G residual magnetic field of
our setup. These fits are shown in the insets to Fig. 2c, and the SM provides
additional plots demonstrating good agreement between the modeled and measured
data SI . The large uncertainty in $\Gamma$ is dominated by uncertainties in
sample alignment and $g$-factor calculations as detailed in the SM SI .
Simultaneous magnetic resonance measurements in a magnetic field could address
these sources of uncertainty by providing in situ calibration of the
$g$-factors for an arbitrary sample orientation.

##  IV Spectral Diffusion in Bulk Si

When cooled to $1.4$ K, the thermal broadening contribution to the ZPL
linewidth is $33$ kHz, and the ZPL PLE linewidth in the 28Si sample narrows to
$\Lambda_{\text{PLE}}=$ $38\pm 6$ MHz SI . In the absence of significant
thermal broadening, we may assume that measuring $\Gamma$ amounts to measuring
the total spectral diffusion linewidth $\Gamma_{\text{sd}}$. To measure
$\Gamma_{\text{sd}}$ at $1.4$ K, we resonantly drive the T centre ZPL and ramp
$B_{\langle 100\rangle}$. Fitting the normalized data shown in Fig. 3a to Eq.
2 averaged over the twelve orientational subsets, we find
$\Gamma_{\text{sd}}=$ $27\pm 12$ MHz for our bulk 28Si sample which agrees
with the measured $\Lambda_{\text{PLE}}$. We stress this model only considers
homogeneous broadening and ignores any inhomogeneities that might be present.
For our measurements of $\Gamma_{\text{sd}}$, we integrate PL counts for $\geq
1$ s at each value of $B_{\langle 100\rangle}$, making this a long-term
average value of the characteristic spectral diffusion within the T centre
ensemble.

To demonstrate that this method can provide a means of measuring ensemble
homogeneous linewidths that are not resolved in an inhomogeneously broadened
PLE spectrum, we repeat this measurement on T centres in our bulk FZ
${}^{\text{nat}}$Si sample, which has an inhomogeneously-broadened optical
linewidth of $\Lambda_{\text{PLE}}=$ $6.0\pm 0.2$ GHz at $1.4$ K. To include
this inhomogeneous broadening in our hyperpolarization model, we replace
$\Delta\rightarrow\Delta+\delta E$ in Eq. 2 to account for inhomogeneously-
shifted resonances and convolve the result with a distribution of $\delta E$
described by a Gauss-Lorentz product with a full width at half maximum equal
to the measured optical linewidth $\Lambda_{\text{PLE}}$. For simplicity, we
restrict ourselves to $\Delta=0$, and once again, we model contributions from
different T centre orientational subsets by averaging the calculated response
over the expected distribution of $g_{h}$.

Figure 3: Spectral diffusion in bulk silicon (a,b) T centre zero phonon line
(ZPL) amplitude as a function of magnetic field at a temperature of 1.4 K for
(a) 28Si and (b) ${}^{\text{nat}}$Si. The data in (a) is fit to the
homogeneously-broadened model given by Eq. 2. The model curves in (b) assume
dynamics dictated either by $\Gamma_{\text{sd}}$ with inhomogeneous
contributions (orange line) or by the optically-measured inhomogeneous
linewidth (green line). The insets in (a,b) show the T centre ZPL measured at
$B_{\langle 100\rangle}\sim 36$ G (blue) and $\sim 450$ G (orange). In the
insets to (a) and (b), detuning is measured from the ZPL resonance at
$935.064$ meV and $935.167$ meV, respectively.

Because our single crystal ${}^{\text{nat}}$Si sample has been cut along an
unknown crystallographic axis, we can only extract an upper bound on
$\Gamma_{\text{sd}}$ by fitting the measured data shown in Fig. 3b as a
function of magnetic field alignment and selecting the orientation that
maximizes $\Gamma_{\text{sd}}$ as detailed in the SM SI . From this maximal
orientation, we find $\Gamma_{\text{sd}}\leq$ $16\pm 7$ MHz (the orange curve
in Fig. 3b). For comparison, we then assume the measured decay is entirely
dictated by inhomogeneous broadening and plot the behaviour predicted by Eq. 2
with $\Gamma$ equal to the inhomogeneous linewidth $\Lambda_{\text{PLE}}=$
$6.0$ GHz. The resulting model curve (shown in green in Fig. 3b) fails to
match the data. As these models are calculated for bounding cases, we can
comfortably conclude that this method extracts linewidths narrower than
$\Lambda_{\text{PLE}}$.

##  V Spectral Diffusion in SOI

Lastly, we use this method to benchmark the optical properties of implanted T
centres in CZ SOI. At $1.4$ K, the SOI samples shown in Fig. 1d display
inhomogeneously-broadened ZPL linewidths of $\Lambda_{\text{PLE}}=$ $20.6\pm
0.3$ GHz and $22.8\pm 0.5$ GHz for the samples implanted with fluences of
$7\times 10^{12}$ cm-2 and $7\times 10^{13}$ cm-2, respectively. We measure
the PLE amplitude as a function of $B_{\langle 100\rangle}$ over the range of
$B_{\langle 100\rangle}$ for which our optical setup remains stable, and we
find $\Gamma_{\text{sd}}=$ $1.0\pm 0.2$ GHz and $1.3\pm 0.3$ GHz,
respectively. Notably, the sample that received the lower implant fluence
shows a slightly reduced level of spectral diffusion. This suggests that at
least some portion of $\Gamma_{\text{sd}}$ could be attributable to our
fabrication process and might be improved with additional surface treatment
Sangtawesin et al. (2019), reduced implant fluence van Dam et al. (2019);
Wolfowicz et al. (2020), or further processing.

Figure 4: Spectral diffusion in silicon-on-insulator (SOI) (a,b) T centre zero
phonon line (ZPL) amplitude as a function of magnetic field at a temperature
of 1.4 K for SOI implanted at a fluence of (a) $7\times 10^{12}$ cm-2 and (b)
$7\times 10^{13}$ cm-2. The insets show the T centre ZPL measured at
$B_{\langle 100\rangle}\sim 36$ G (blue) and $\sim 2000$ G (orange). In the
insets to (a) and (b), detuning is measured from the ZPL resonance at
$935.176$ meV and $935.237$ meV, respectively.

It is important to note that the assumption of equally weighted contributions
from each T centre orientation is likely not true. To address this, we have
also fit our data under the assumption that the measured signal comes entirely
from the T centre orientation that minimizes (maximizes) $\Gamma_{\text{sd}}$
to obtain the bounding sets $\Gamma_{\text{sd}}\in[0.8\pm 0.3,1.7\pm 0.1]$ GHz
and $\Gamma_{\text{sd}}\in[1.0\pm 0.3,2.1\pm 0.1]$ GHz for the $7\times
10^{12}$ cm-2 and $7\times 10^{13}$ cm-2 samples, respectively. Once again,
simultaneous magnetic resonance measurements in a magnetic field could improve
the precision of this technique by calibrating the relative signal
contributions of different orientations.

Even with this uncertainty, however, this method provides efficient
comparative measurements of $\Gamma_{\text{sd}}$ and rapid feedback on
fabrication processes designed to reduce $\Gamma_{\text{sd}}$. Measuring PLE
decay due to spin-dependent optical transitions can provide the characteristic
$\Gamma_{\text{sd}}$ within an ensemble of emitters when complementary
techniques such as spectral hole burning Völker (1989) are not possible.
Moreover, because our measurements interrogate a large ensemble of T centres
simultaneously, the technique developed here enables quick access to the
characteristic value of $\Gamma_{\text{sd}}$ across a sample, which can be
time consuming to obtain by measuring individual centres.

##  VI Outlook

Although far from the lifetime-limited linewidth of $\Gamma_{1}=169$ kHz
Bergeron et al. (2020) and $\sim 100$ times larger than the values we report
in bulk samples, these values for $\Gamma_{\text{sd}}$ of implanted T centres
in CZ SOI are already competitive with other solid state emitters near
interfaces. It is common for spectral diffusion in defect centres that lack
inversion symmetry to be on the order of several GHz in the absence of careful
surface treatment or experimental protocols Ishikawa et al. (2012); Evans et
al. (2016); Anderson et al. (2019); van Dam et al. (2019); Crook et al.
(2020). With minimal optimization, implanted T centres display $\sim 1$ GHz
characteristic spectral diffusion, and this could be minimized further as
fabrication and measurement techniques continue to develop. Moreover, the
characteristic spectral diffusion reported here is a long-term, ensemble
average. Individual centres within this ensemble could have much lower
$\Gamma_{\text{sd}}$ Chu et al. (2014); Kasperczyk et al. (2020), and any
slowly varying effects that contribute to this long-term average might be
addressed by feedback or filtering techniques Bernien et al. (2013); Guo et
al. (2019).

To estimate the effect spectral diffusion will have on the performance of a T
centre-based photon-spin interface, we calculate the indistinguishability of
two photons emitted by a given T centre as
$I=\xi\Gamma_{1}/(\xi\Gamma_{1}+\Gamma_{\text{sd}})$ where $\xi=0.23$ is the T
centre Debye-Waller factor Bylander et al. (2003); Grange et al. (2015);
Bergeron et al. (2020). For our bulk and SOI samples, respectively, we find
$I\sim 2\times 10^{-3}$ and $I\sim 4\times 10^{-5}$. We stress that both of
these values can likely be improved by employing fabrication or measurement
techniques designed to stabilize the electromagnetic environment de Oliveira
et al. (2017); Anderson et al. (2019); Sangtawesin et al. (2019); Kasperczyk
et al. (2020) or by filtering the emitted photons Bernien et al. (2013); Guo
et al. (2019). Moreover, embedding T centres inside optical resonators can
Purcell enhance $\Gamma_{1}$ to increase the photon indistinguishabilty.
Optical resonators in SOI capable of demonstrating Purcell factors above
$10^{4}$ have already been demonstrated Zain et al. (2008). Because the T
centre radiative efficiency is expected to be near-unity Bergeron et al.
(2020), T centres incorporated into such an optical cavity could provide high
quality photon-spin interfaces in the telecommunications $O$-band.

In summary, we have demonstrated a reliable method for using ion implantation
and annealing treatments to incorporate T centres into industry-standard
silicon integrated photonics material. We then developed an efficient method,
which could be applied to a wide variety of colour centres, for using spin-
dependent optical transitions to benchmark spectral diffusion. Finally, we
used this new technique to show that implanted T centres within a $220$ nm SOI
device layer can display GHz-scale levels of spectral diffusion with minimal
optimization. These results represent an important step towards developing
quantum technologies in silicon.

##  VII Acknowledgments

We thank M. L. W. Thewalt and L. Childress for fruitful discussions as well as
N. V. Abrosimov for bulk sample growth and C. Clément for rapid thermal
annealing.

##  VIII Funding

This work was supported by the Natural Sciences and Engineering Research
Council of Canada (NSERC), the Canada Research Chairs program (CRC), the
Canada Foundation for Innovation (CFI), the B.C. Knowledge Development Fund
(BCKDF), the Canadian Institute for Advanced Research (CIFAR) Quantum
Information Science program, the CIFAR Catalyst Fund, and Le Fonds de
recherche du Québec – Nature et technologies (FRQNT). The 28Si samples used in
this study were prepared from the Avo28 crystal produced by the International
Avogadro Coordination (IAC) Project (2004–2011) in cooperation among the BIPM,
the INRIM (Italy), the IRMM (EU), the NMIA (Australia), the NMIJ (Japan), the
NPL (UK), and the PTB (Germany).

##  IX References

## References

  * Awschalom et al. (2018) D. D. Awschalom, R. Hanson, J. Wrachtrup, and B. B. Zhou, Nature Photonics 12, 516 (2018).
  * Zhang et al. (2020) G. Zhang, Y. Cheng, J.-P. Chou, and A. Gali, Applied Physics Reviews 7, 031308 (2020).
  * Buckley et al. (2017) S. Buckley, J. Chiles, A. N. McCaughan, G. Moody, K. L. Silverman, M. J. Stevens, R. P. Mirin, S. W. Nam, and J. M. Shainline, Applied Physics Letters 111, 141101 (2017).
  * Chartrand et al. (2018) C. Chartrand, L. Bergeron, K. J. Morse, H. Riemann, N. V. Abrosimov, P. Becker, H.-J. Pohl, S. Simmons, and M. L. W. Thewalt, Phys. Rev. B 98, 195201 (2018).
  * Beaufils et al. (2018) C. Beaufils, W. Redjem, E. Rousseau, V. Jacques, A. Y. Kuznetsov, C. Raynaud, C. Voisin, A. Benali, T. Herzig, S. Pezzagna, et al., Phys. Rev. B 97, 035303 (2018).
  * Redjem et al. (2020) W. Redjem, A. Durand, T. Herzig, A. Benali, S. Pezzagna, J. Meijer, A. Y. Kuznetsov, H. S. Nguyen, S. Cueff, J. M. Gérard, et al., Nat. Electron. 3, 738 (2020).
  * Bergeron et al. (2020) L. Bergeron, C. Chartrand, A. T. K. Kurkjian, K. J. Morse, H. Riemann, N. V. Abrosimov, P. Becker, H.-J. Pohl, M. L. W. Thewalt, and S. Simmons, PRX Quantum 1, 020301 (2020).
  * Meijer et al. (2005) J. Meijer, B. Burchard, M. Domhan, C. Wittmann, T. Gaebel, I. Popa, F. Jelezko, and J. Wrachtrup, Applied Physics Letters 87, 261909 (2005).
  * Chu et al. (2014) Y. Chu, N. de Leon, B. Shields, B. Hausmann, R. Evans, E. Togan, M. J. Burek, M. Markham, A. Stacey, A. Zibrov, et al., Nano Letters 14, 1982 (2014).
  * Falk et al. (2013) A. L. Falk, B. B. Buckley, G. Calusine, W. F. Koehl, V. V. Bobrovitski, A. Politi, C. A. Zorman, P. X.-L. Feng, and D. D. Awschalom, Nat. Commun. 4, 1819 (2013).
  * Phenicie et al. (2019) C. M. Phenicie, P. Stevenson, S. Welinski, B. C. Rose, A. T. Asfaw, R. J. Cava, S. A. Lyon, N. P. de Leon, and J. D. Thompson, Nano Letters 19, 8928 (2019).
  * Buckley et al. (2020) S. M. Buckley, A. N. Tait, G. Moody, B. Primavera, S. Olson, J. Herman, K. L. Silverman, S. P. Rao, S. W. Nam, R. P. Mirin, et al., Opt. Express 28, 16057 (2020).
  * (13) A. T. K. Kurkjian, D. B. Higginbottom, et al, in preparation.
  * Muhonen et al. (2014) J. T. Muhonen, J. P. Dehollain, A. L. ans Fay E. Hudson, R. Kalra, T. Sekiguchi, K. M. Itoh, D. N. Jamieson, J. C. McCallum, A. S. Dzurak, and A. Morello, Nature Nanotech. 9, 986 (2014).
  * de Oliveira et al. (2017) F. F. de Oliveira, D. Antonov, Y. Wang, P. Neumann, S. A. Momenzadeh, T. Häußermann, A. Pasquarelli, A. Denisenko, and J. Wrachtrup, Nat Commun 8, 15409 (2017).
  * Anderson et al. (2019) C. P. Anderson, A. Bourassa, K. C. Miao, G. Wolfowicz, P. J. Mintun, A. L. Crook, H. Abe, J. Ul Hassan, N. T. Son, T. Ohshima, et al., Science 366, 1225 (2019).
  * Sangtawesin et al. (2019) S. Sangtawesin, B. L. Dwyer, S. Srinivasan, J. J. Allred, L. V. H. Rodgers, K. De Greve, A. Stacey, N. Dontschuk, K. M. O’Donnell, D. Hu, et al., Phys. Rev. X 9, 031052 (2019).
  * Kasperczyk et al. (2020) M. Kasperczyk, J. A. Zuber, A. Barfuss, J. Kölbl, V. Yurgens, S. Flågan, T. Jakubczyk, B. Shields, R. J. Warburton, and P. Maletinsky, Phys. Rev. B 102, 075312 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.102.075312.
  * Bernien et al. (2013) H. Bernien, B. Hensen, W. Pfaff, G. Koolstra, M. S. Blok, L. Robledo, T. H. Taminiau, M. Markham, D. J. Twitchen, L. Childress, et al., Nature 497, 86 (2013).
  * Guo et al. (2019) K. Guo, H. Springbett, T. Zhu, R. A. Oliver, Y. Arakawa, and M. J. Holmes, Appl. Phys. Lett. 114, 112109 (2019).
  * Grange et al. (2015) T. Grange, G. Hornecker, D. Hunger, J.-P. Poizat, J.-M. Gérard, P. Senellart, and A. Auffèves, Phys. Rev. Lett. 114, 193601 (2015).
  * Giesz et al. (2015) V. Giesz, S. L. Portalupi, T. Grange, C. Antón, L. De Santis, J. Demory, N. Somaschi, I. Sagnes, A. Lemaître, L. Lanco, et al., Phys. Rev. B 92, 161302 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.92.161302.
  * Safonov et al. (1996) A. N. Safonov, E. C. Lightowlers, G. Davies, P. Leary, R. Jones, and S. Öberg, Phys. Rev. Lett. 77, 4812 (1996).
  * (24) See Supplementary Information.
  * Pichler (2004) P. Pichler, _Isovalent Impurities_ (Springer Vienna, Vienna, 2004), pp. 281–329, ISBN 978-3-7091-0597-9.
  * Berhanuddin et al. (2012) D. D. Berhanuddin, M. A. Lournco, R. M. Gwilliam, and K. P. Homewood, Adv. Funct. Mater. 22, 2709 (2012).
  * Dash and Newman (1955) W. C. Dash and R. Newman, Phys. Rev. 99, 1151 (1955).
  * Jellison Jr. and Joshi (2018) G. E. Jellison Jr. and P. C. Joshi, _Spectroscopic Ellipsometry for Photovoltaics_ (Springer, 2018), vol. 1, chap. 8.
  * Barnes (1998) W. L. Barnes, Journal of Modern Optics 45, 661 (1998).
  * Saeedi et al. (2013) K. Saeedi, S. Simmons, J. Z. Salvail, P. Dluhy, H. Riemann, N. V. Abrosimov, P. Becker, H.-J. Pohl, J. J. L. Morton, and M. L. W. Thewalt, Science 342, 830 (2013).
  * van Dam et al. (2019) S. B. van Dam, M. Walsh, M. J. Degen, E. Bersin, S. L. Mouradian, A. Galiullin, M. Ruf, M. IJspeert, T. H. Taminiau, R. Hanson, et al., Phys. Rev. B 99, 161203 (2019).
  * Ishikawa et al. (2012) T. Ishikawa, K.-M. C. Fu, C. Santori, V. M. Acosta, R. G. Beausoleil, K. Watanabe, S. Chikata, and K. M. Itoh, Nano Lett. 12, 2083 (2012).
  * Evans et al. (2016) R. E. Evans, A. Sipahigil, D. D. Sukachev, A. S. Zibrov, and M. D. Lukin, Phys. Rev. Applied 5, 044010 (2016).
  * Crook et al. (2020) A. L. Crook, C. P. Anderson, K. C. Miao, A. Bourassa, H. Lee, S. L. Bayliss, D. O. Bracher, X. Zhang, H. Abe, T. Ohshima, et al., Nano Letters 20, 3427 (2020).
  * Safonov et al. (1995) A. N. Safonov, E. C. Lightowlers, and G. Davies, Materials Science Forum 196-201, 909 (1995).
  * Wolfowicz et al. (2020) G. Wolfowicz, F. J. Heremans, C. P. Anderson, S. Kanai, H. Seo, A. Gali, G. Galli, and D. D. Awschalom, arXiv:2010.16395 (2020).
  * Völker (1989) S. Völker, Annu. Rev. Phys. Chem. 40, 499 (1989).
  * Bylander et al. (2003) J. Bylander, I. Robert-Philip, and I. Abram, Eur. Phys. J. D 22, 295 (2003).
  * Zain et al. (2008) A. R. M. Zain, N. P. Johnson, M. Sorel, and R. M. D. L. Rue, Opt. Express 16, 12084 (2008).

##  X Supplementary Information

###  X.1 Measurement Details

The samples are mounted in a liquid helium cryostat. The sample temperature is
set by pumping on the liquid helium bath. The PL and PLE spectroscopy
measurements were performed as reported in Ref. Bergeron et al. (2020). For
PLE spectroscopy measurements a Toptica DL100 tunable diode laser was used for
measurements of 28Si, and a tunable nanoplus DFB laser diode was used for the
broader SOI and bulk ${}^{\text{nat}}$Si samples.

###  X.2 Implantation Dose Optimization

To select an implant dose for our T centre implantation study, we first
generated T centres in float zone (FZ) silicon wafers by co-implanting carbon
and hydrogen. We swept the fluence of the implants while keeping fixed C:H
ratios of 1:1, 2:1, and 1:0. Carbon implants were again performed at an energy
of $38$ keV, and hydrogen was implanted at $9$ keV. The samples were then
annealed in air following a $(100,300,350,400,450)^{\circ}$C step-wise
temperature profile with dwell times of $(60,30,30,60,60)$ min. This
temperature profile was chosen as it previously generated high density T
centre concentrations in electron-irradiated, high-carbon bulk samples.

Fig. 5a shows the integrated area of the FE-normalized T centre ZPL for each
fluence we measured, and Fig. 5b presents the PL spectra for the implant
fluences that generated the brightest T centre ZPL for each C:H ratio. Each
spectrum shows a strong T centre ZPL, but the 1:0 C:H sample has a much higher
broadband radiation damage background. By comparison, the 1:1 and 2:1 C:H
implants show comparable T centre intensity with a much lower background for
carbon fluences of $7\times 10^{12}$ cm-2. From this data, we decided to fix
the C:H ratio at 1:1 and implant at fluences of $7\times 10^{12}$ cm-2 and
$7\times 10^{13}$ cm-2.

Figure 5: T centre implantation development. (a) Integrated T centre zero
phonon line (ZPL) area as a function of implant dose for fixed carbon-to-
hydrogen ratios. (b) T centre photoluminescence (PL) spectra for carbon and
hydrogen co-implantation into bulk float zone silicon.

###  X.3 T Centre Localization

To confirm the generated T centres are localized within the SOI device layer,
we etched away the device layer from one of our SOI samples and measured PL
from the exposed buried oxide and handle wafer material. As shown in Fig. 6a,
the resulting $532$ nm PL spectrum shows no significant T centre PL as
compared to the same sample in an unetched region, suggesting that the handle
material does not contain a significant concentration of T centres.

Figure 6: T centre localization. (a) Photoluminescence (PL) spectra of an
implanted sample that has had a portion of its device layer etched away. (b)
PL spectra of implanted SOI samples under $532$ nm and UV excitation.

As further confirmation that the generated T centres are localized within the
SOI device layer, we performed PL spectroscopy with UV excitation. Because the
absorption depth of $532$ nm light in Si is $2.5$ µm at cryogenic
temperatures, the SOI PL spectra presented in the main text sample both the
SOI device layer and the silicon handle wafer beneath the buried oxide. UV
wavelengths have a penetration depth of $<1$ µm in Si at cryogenic
temperatures Dash and Newman (1955); Jellison Jr. and Joshi (2018), making PL
with UV excitation more sensitive to the device layer than PL with $532$ nm
excitation. Fig. 6b shows a representative PL spectrum measured on the same
sample with both $532$ nm excitation and UV excitation. Because the UV
measurement shows a negligible FE recombination line, we have normalized both
spectra to their integrated signal. Both spectra show strong T centre lines,
confirming that there are T centres within the device layer.

Taken together, these two measurements localize our implanted T centres to
within the device layer of the SOI materal.

###  X.4 Calculating g-factors

Our calculation of the hole $g$-factors expands upon the work of Ref. Safonov
et al. (1995). We write the hole Hamiltonian as the sum of a strain and a
magnetic term $H_{h}=H_{s}(\epsilon)+H_{b}(B)$ where the strain Hamiltonian is

$H_{s}(\epsilon)=-b\sum_{i}(J_{i}^{2}-I)\epsilon_{ii}-\frac{d}{\sqrt{3}}\sum_{i\neq
j}[J_{i},J_{j}]\epsilon_{ij}.$ (3)

Here, $J_{i}$ are the angular momentum operators,
$[J_{i},J_{j}]=\frac{1}{2}(J_{i}J_{j}+J_{j}J_{i})$, $\epsilon_{ij}$ are strain
tensor components, and $b$, $d$ are deformation parameters. An external
magnetic field $\vec{B}$ enters the Hamiltonian as

$H_{b}(\vec{B})=\mu_{B}\left(g_{1}\sum_{i}B_{i}J_{i}+g_{2}\sum_{i}B_{i}J_{i}^{3}\right).$
(4)

where $\mu_{B}$ is the Bohr magneton and $g_{1}$, $g_{2}$ are the hole
$g$-factors.

We use the internal strain and deformation parameters obtained by Ref. Safonov
et al. (1995): $\epsilon_{YY}=-0.65\times 10^{-3}$, $\epsilon_{ZZ}=-0.26\times
10^{-3}$, $b=-0.8$ eV, and $d=-2.7$ eV. The internal strain parameters were
defined relative to the defect axes where $\hat{y}$ is parallel to $\langle
110\rangle$ and $\hat{z}$ is in the $(110)$ plane, tilted from the $\langle
001\rangle$ axis by $4^{\circ}$. In our treatment, we rotate the internal
strain from the defect coordinate system to that of the crystal so that
magnetic field alignments can be referenced to crystallographic axes.

Rather than using the $g_{1}=1.3$ and $g_{2}=-0.1$ values from Ref. Safonov et
al. (1995), we obtain $g_{1}=1.505\pm 0.022$ and $g_{2}=-0.138\pm 0.012$ by
fitting the model above to the $g$-factors obtained in Ref. Bergeron et al.
(2020), restricting ourselves to $g_{i}$ values near those reported in Ref.
Safonov et al. (1995). The highly-resolved 28Si data from Ref. Bergeron et al.
(2020) enables us to fit the hole $g$-factors using all twelve T centre
orientational subsets. For this fit, we also introduce a misalignment of the
magnetic field as an additional free parameter. Our fit returns an angular
inclination of $10.1^{\circ}\pm 0.1^{\circ}$ and rotation of $71.9^{\circ}\pm
0.6^{\circ}$ about the nominal $\langle 110\rangle$ field alignment in that
work. Table 1 provides comparisons between our fit results and the $g_{h}$
values reported in Ref. Bergeron et al. (2020) for a $\langle 110\rangle$
magnetic field alignment.

Subset | Measured $g_{h}^{\langle 110\rangle}$ Bergeron et al. (2020) | Fitted $g_{h}^{\langle 110\rangle}$ | Calculated $g_{h}^{\langle 100\rangle}$
---|---|---|---
1 | $3.457\pm 0.007$ | $3.460\pm 0.084$ | $0.91\pm 0.10$
2 | $3.457\pm 0.007$ | $3.459\pm 0.084$ | $0.91\pm 0.11$
3 | $2.233\pm 0.009$ | $2.269\pm 0.057$ | $0.91\pm 0.16$
4 | $2.165\pm 0.014$ | $2.252\pm 0.057$ | $0.91\pm 0.16$
5 | $1.970\pm 0.012$ | $2.036\pm 0.052$ | $2.55\pm 0.42$
6 | $1.871\pm 0.022$ | $2.015\pm 0.052$ | $2.55\pm 0.42$
7 | $1.851\pm 0.014$ | $1.779\pm 0.049$ | $2.55\pm 0.13$
8 | $1.770\pm 0.008$ | $1.758\pm 0.049$ | $2.55\pm 0.13$
9 | $1.596\pm 0.006$ | $1.566\pm 0.044$ | $2.55\pm 0.08$
10 | $1.497\pm 0.011$ | $1.537\pm 0.045$ | $2.55\pm 0.08$
11 | $1.082\pm 0.007$ | $1.029\pm 0.024$ | $2.55\pm 0.34$
12 | $1.069\pm 0.007$ | $1.023\pm 0.024$ | $2.55\pm 0.34$
Table 1: Comparing hole $g$-factors for a $\langle 110\rangle$ magnetic field
orientation as measured in Ref. Bergeron et al. (2020) to those fitted by the
model introduced in Ref. Safonov et al. (1995).

With model parameters in hand, we then calculate $g_{h}$ for each T centre
orientational subset by finding the $B_{i}$ components for each subset,
computing the eigenvalues of $H_{h}$, and dividing the resulting spin
splitting by $\mu_{B}B$. $H_{h}$ predicts a negligible zero-field splitting
within the $\text{TX}_{0}$ orbital branch and negligible interaction with the
$\text{TX}_{1}$ orbital at the moderate magnetic fields used in this work,
justifying this approach for calculating the $g_{h}$ values.

The $g$-factors for $\vec{B}\parallel\langle 100\rangle$ used in the main text
are listed with their calculated error bars in Table 1. For these
calculations, we have assigned inclination and azimuthal $\vec{B}$-field
alignment errors of $\pm 10^{\circ}$ each.

Figure 7: Bounding calculation of $\Gamma_{\text{sd}}$ in bulk
${}^{\text{nat}}$Si. (a,b) Fitted values of $\Gamma_{\text{sd}}$ as a function
of magnetic field orientation. The orientations in white were unable to fit
the data and were excluded from the plot. (b) A higher resolution sweep of
orientations in the red square of (a) that are near the maximal
$\Gamma_{\text{sd}}$.

For our bounding calculation of $\Gamma_{\text{sd}}$ in
${}^{\text{nat}}\text{Si}$, we fit the measured data as a function of the
magnetic field alignment. The fitted values for $\Gamma_{\text{sd}}$ are
plotted as a function of the magnetic field orientation in Fig. 7. Here, we
have omitted magnetic field orientations for which the model is unable to fit
the measured data. The model curve in Fig. 3b of the main text is plotted
using inclination and azimuthal angles of $0.30\pi$ and $0.25\pi$,
respectively.

###  X.5 Optical Transitions Between Different Spin Orientations

In the model presented in the main text, we have neglected optical transitions
between different spin orientations. Including these transitions into the rate
model used to derive Eq. 1 of the main text, we can rewrite the amplitude
proportionality as

$A\propto\frac{L_{B}L_{C}+L_{A}L_{C}+L_{B}L_{D}+L_{A}L_{D}}{L_{A}+L_{B}+L_{C}+L_{D}}$
(5)

where the optical transitions are labeled by subscripts as in Fig. 8a. The
amplitudes $L_{B}$, $L_{C}$ have been measured to be $\sim 10$ times larger
than $L_{A}$, $L_{D}$ Bergeron et al. (2020); kur . Moveover, since
$\left|\frac{\partial L_{A}}{\partial B_{z}}\right|,\left|\frac{\partial
L_{D}}{\partial B_{z}}\right|>\left|\frac{\partial L_{B}}{\partial
B_{z}}\right|,\left|\frac{\partial L_{C}}{\partial B_{z}}\right|$, the
approximation given by Eq. 1 of the main text becomes better as $B_{z}$
increases.

To quantify the effect of the $L_{A}$ and $L_{D}$ optical transitions on our
hyperpolarization model, we fit the thermally-broadened $T=4.2$ K data
displayed in Fig. 2c of the main text using the four-transition model of Eq.
5. Ignoring inhomogeneous broadening, we find $\Gamma=250\pm 100$ MHz with
$L_{A,D}=0$ and $\Gamma=260\pm 100$ MHz with $L_{A,D}=0.1L_{B,C}$. To further
explore this dependency, we fix $L_{A,D}+L_{B,C}=1$, fit the measured data as
a function of $L_{A,D}$, and compare the homogeneous linewidth extracted from
the model to the optically measured linewidth $\Lambda_{\text{PLE}}=$
$329.8\pm 0.6$ MHz. Plotting the results in Fig. 8b, we see that $\Gamma$ only
becomes strongly dependent on $L_{A,D}$ when $L_{A,D}>L_{B,C}$, as expected.
In this regime, the hyperpolarization model returns values of $\Gamma$ that
are much larger than the optically-measured value, supporting the assertion
that $L_{A,D}\ll L_{B,C}$ that has been experimentally verified elsewhere
Bergeron et al. (2020); kur .

Figure 8: Effect of optical transitions between different spin orientations.
(a) T centre optical transitions. (b) T centre homogeneous linewidth extracted
from the $4.2$ K hyperpolarization data presented in the text. Good agreement
between the optical and hyperpolarization measurements is only found when
transitions between parallel spin states dominate ($L_{B,C}>L_{A,D}$).

###  X.6 Fitting residual-field PLE spectra

The poles of the electromagnet used in these measurements generated a residual
field of $B_{0}\sim 36$ G. For our 28Si sample, the Zeeman splitting from this
residual field is comparable to the measured optical linewidth
$\Lambda_{\text{PLE}}$. To account for this in our analysis, we fit the
residual field spectra to an ensemble average of Eq. 2 with $B_{z}$ along the
$\langle 100\rangle$ axis and take $\Gamma=\Lambda_{\text{PLE}}$.

For the ${}^{\text{nat}}$Si and SOI samples where the lineshape is better
described by the Gauss-Lorentz product (GLP), we define the GLP linewidth to
be $\Gamma=\sqrt{\Lambda_{\text{PLE}}^{2}+\varepsilon_{B}(B_{z})^{2}}$.

###  X.7 Hyperpolarization model linecuts

In Fig. 2c of the main text, we measure T centre hyperpolarization in our 28Si
sample at $4.2$ K and compare the measured data to the behaviour predicted by
the hyperpolarization model developed in the text. As additional comparison,
we plot the simulated 2d data in Fig. 9a. We then fit both the measured and
the modeled PLE spectra to Lorentzian lineshapes and plot the fitted
linewidths against each other in Fig. 9b.

Figure 9: Hyperpolarization in T centres (a) Simulated version of Fig. 2c in
the main text. (b) Predicted and measured PLE linewidth as a function of
magnetic field for our 28Si sample at 4.2 K.


# Boosting the Search Performance of B+-tree for Non-volatile Memory with
Sentinels

Chongnan Ye  School of Information Science and Technology
ShanghaiTech University
<EMAIL_ADDRESS>Chundong Wang *Corresponding Author School of
Information Science and Technology
ShanghaiTech University
<EMAIL_ADDRESS>
###### Abstract

The next-generation non-volatile memory (NVM) is striding into computer
systems as a new tier as it incorporates both DRAM’s byte-addressability and
disk’s persistency. Researchers and practitioners have considered building
persistent memory by placing NVM on the memory bus for CPU to directly load
and store data. As a result, cache-friendly data structures have been
developed for NVM. One of them is the prevalent B+-tree. State-of-the-art in-
NVM B+-trees mainly focus on the optimization of write operations (insertion
and deletion). However, search is of vital importance for B+-tree. Not only
search-intensive workloads benefit from an optimized search, but insertion and
deletion also rely on a preceding search operation to proceed. In this paper,
we attentively study a sorted B+-tree node that spans over contiguous cache
lines. Such cache lines exhibit a monotonically increasing trend and searching
a target key across them can be accelerated by estimating a range the key
falls into. To do so, we construct a probing Sentinel Array in which a
sentinel stands for each cache line of B+-tree node. Checking the Sentinel
Array avoids scanning unnecessary cache lines and hence significantly reduces
cache misses for a search. A quantitative evaluation shows that using Sentinel
Arrays boosts the search performance of state-of-the-art in-NVM B+-trees by up
to 48.4% while the cost of maintaining of Sentinel Array is low.

###### Index Terms:

Non-volatile memory, B+-tree, Key-value Database, Cache

##  I Introduction

A conventional computer employs DRAM as main memory and hard disk as storage.
The disk provides persistency but demands a long time to load and store data
[1]. DRAM is volatile and cannot reach large capacity because of its density
limitation [2, 3, 4]. Recently, the next-generation non-volatile memory (NVM),
such as phase-change memory (PCM) [5], spin-transfer torque memory (STT-RAM)
[6], and resistive RAM (RRAM) [7], emerges as a promising medium that combines
both DRAM’s byte-addressability and disk’s persistency. Hence, researchers and
practitioners have proposed to place NVM on the memory bus to build persistent
memory, which acts as a new tier of computer architecture and blurs the
boundary between DRAM and disk.

Multiple classic data structures have been tuned for NVM. Among them B+-tree
is one that has been widely used for databases and file systems since the era
of disks [8]. B+-tree is a multi-level structure composed of internal and leaf
nodes. Internal nodes (INs) are used to index lower-level INs and LNs while
leaf nodes (LNs) hold actual key-value (KV) pairs. INs and LNs have the same
structure but the pointers of INs point to the address of LNs while the
pointers of LNs point to the value data that stored in the B+-tree. And INs
have one more pointer than LNs which indicates the leftmost LN of the internal
nodes. Keys are sorted in both INs and LNs. When tuning B+-tree for NVM,
researchers need to consider the fact that CPU or memory controller may
reorder writes from CPU cache to NVM [9, 10, 11]. One way to retain a writing
order is to employ a combination of cache line flush and memory fence [12],
which, however, incurs significant performance overheads. With further and
further understanding of B+-tree’s properties, state-of-the-art B+-trees
designed for NVM have evolved from being with sorted nodes (e.g., CDDS-Tree
[13]) to unsorted nodes (e.g., NV-Tree, wB+-tree, and FPTree [14]), and again
sorted nodes (e.g., FAST-FAIR and Circ-Tree). When developing such B+-tree
variants, researchers mainly placed their emphasis on reducing the use of
cache line flush and memory fence in order to achieve a high write performance
for insertion and deletion.

Nevertheless, the search operation with a target key for the key’s
corresponding value is of paramount importance for B+-tree. Search is not only
an independent operation, but insertion and deletion also rely on a preceding
search to proceed. As the latest FAST-FAIR and Circ-Tree are with sorted
nodes, we focus on optimizing the search operation over sorted KV pairs stored
in NVM. It is revealed that linear search is more efficient than binary search
over sorted keys although the latter has a theoretically lower time
complexity. The reason is that binary search incurs more cache misses while
linear search facilitates CPU’s branch prediction and prefetch. Our aim of
this paper is to further reduce cache misses in searching. We find that no
matter if we use binary search or linear search, we must access and check keys
that are not the target one. Although accessing those keys leads us to
continue the search, it is the cause of cache misses. This is the observation
of read amplification in searching a B+-tree node. We intend to reduce as much
read amplification as possible. In other words, we want to timely find the
exact cache line in which the target key resides. It is impossible except
using some auxiliary information. Previously, Oukid et al. [14] proposed to
use 1B hash values as fingerprints for keys. A key and its fingerprint have
the same offset. Searching a key transforms to finding the key’s fingerprint
that indicates the location of the key. However, due to the collisions of 1B
fingerprints and the calculation cost, fingerprinting is not efficient in
practice.

We closely study a sorted B+-tree node. We find that each cache line has a
range of keys and all cache lines exhibit a monotonically increasing trend. We
can choose a sentinel key in a uniform fashion, e.g., the smallest or the
largest key, from each cache line to reflect such a trend. When searching a
target key, instead of scanning cache lines, we go through the Sentinel Array
to determine which cache line the target key stays in. As a result, searching
a target key converts to searching the Sentinel Array and one actual cache
line of KV pairs. The Sentinel Array contains representatives of cache lines
of a B+-tree node and it thus occupies much fewer cache lines. Also, a linear
scan of it facilitates CPU’s prefetch and branch prediction. Consequently,
employing the Sentinel Array significantly reduces read amplification.

We have augmented FAST-FAIR and Circ-Tree with the idea of Sentinel Array.
Experimental results show that compared to the use of Sentinel Array boosts
the search performance by up to 42.6% and 48.4%, respectively, for FAST-FAIR
and Circ-Tree. Moreover, using Sentinel Arrays is up to 60.1% faster than
using the idea of fingerprints. Maintaining Sentinel Arrays, however, incurs
at most 4.1% and 6.5% more time for insertions with FAST-FAIR and Circ-Tree,
respectively.

The rest of this paper is organized as follows. The background and motivation
are discussed in Section 2. Section 3 shows the detailed design and
implementation of Sentinel Array. The evaluation and results are presented in
Section 4. In the last, we conclude this paper and look forward to future
work.

##  II Background and Motivation

###  II-A Background

Byte-addressable non-volatile memory (NVM), such as Phase Change Memory (PCM),
Spin-Transfer Torque RAM (STT-RAM), and Resistive RAM (RRAM), has an
increasingly important demand in high-performance computing and big data
analysis [15, 12]. It combines the durability of hard disk and the comparable
access performance of DRAM.

Data consistency between the volatile device and the non-volatile device is
always one of the most crucial features of storage systems. In the traditional
storage architecture, where DRAM acts as the main memory, the system needs to
ensure the data consistency from DRAM to the hard disk to prevent the data
from being lost if system failure happens [16]. As to NVM, the effort becomes
to keep data consistency between CPU and NVM. However, modern CPUs mainly
support an atomic write of no more than the memory bus width (8 bytes for
64-bit CPUs) when NVM is put on the memory bus [10, 17, 18]. Moreover,
compared to the known persistent status of each page in DRAM, the status of
data in cache-lines is unknown and the CPU or the memory controller may write
the cache-lines back to the main memory in a different order from the
programmed order. Make memory writes are in a certain order is critical to the
crash consistency of pointer-based data structures. For example, a new B+-tree
from a split operation must be written completely before the pointer to it is
added to its parental node. The pointer may turn to be dangling if the order
is reversed but the crash occurs [9] and it’s necessary to maintain the
writing order.

Some CPU instructions, such as cache line flush and memory fence, are provided
to maintain the written order from cache to memory. cache line flush
explicitly invalidates a dirty cache-line and flush it to memory. memory fence
makes a barrier that holds back the memory operations after the barrier until
those before the barrier complete. However, these instructions incur
performance overheads [9, 10, 19]. Therefore, it is necessary to reduce the
use of cache line flush and memory fence.

###  II-B Related Work

Researchers have proposed a number of mechanisms to provide data consistency
and optimize the performance of B+-trees developed for NVM. FAST-FAIR [10]
uses two features that CPU only ensures 8-bytes atomic write and duplicate
pointers are impossible in B+-tree node. It exploits the store dependencies in
shifting a sequence of KV pairs in a node on insertion/deletion. (e.g.
$KV_{i}\xrightarrow{}KV_{i+1}$), which reduces the cost of cache line flushes
and memory fence.

FAST-FAIR still follows the classic B+-tree node in a linear structure.
Circle-Tree [11] was proposed to reduce write amplification caused by shifting
KV pairs in the linear structure. The circular structure supports
bidirectional shifting. KV pair is inserted/deleted at the left or right side
by deciding which side requires fewer shifting operations.

Most B+-trees developed for NVM optimize insertion/deletion performance.
Although FP-Tree [14] use a technique named fingerprinting, which maps the key
of KV pair into a byte of fingerprint (hash value), to optimize the search
performance by checking the array of fingerprint instead of the original data
container. The efficacy of fingerprinting is low due to hash calculations and
the collisions of 1B hash values.

###  II-C Motivation

Figure 1: B+-tree with different search way comparison

Figure 2: Sentinel Array example

With further and further understanding of B+-tree’s properties, researchers
have switched from sorted nodes to unsorted nodes, and again to sorted nodes
now. More importantly, the researcher found that a large CPU cache of modern
processor makes linear search outperform binary search in a sorted B+-tree
node because of the former’s good cache line locality [20]. Fig. 1 verifies
that the latency of binary search is greater than that of linear search with
small node size (512B/1KB/2KB) when we searched 1 million KV pairs in a sorted
B+-tree node.

As the search is so important that used not only in independent operation, but
also in insertion and deletion, optimize the search performance is a key point
to boost the throughput. Accessing the key and check if that is the target one
is the cause of cache misses. This is the observation of read amplification in
searching a sorted B+-tree node. Our goal focus on reducing cache misses in
searching. Fig. 2 shows a sorted B+-tree node with a Sentinel Array which has
the smallest key from each cache line. When searching a target key 130,
instead of scanning 4 cache lines, we go through the Sentinel Array to
determine the 4th cache is the target key stays in. As a result, searching key
130 converts to searching the Sentinel Array and one actual cache line of KV
pairs.

##  III Design and Implementation of Sentinel Array

Sentinel Array is an auxiliary structure that contains sentinels, which
indicate the minimum key of each cache line of the origin array. The Sentinel
Array size is decided by the machine and the original array size. Given a
typical cache line size of 64 bytes and a sentinel (minimum key) in 8B, one
cache line contains eight sentinels, which are respectively from eight cache
lines of KV pairs. Therefore, a Sentinel Array in one cache line is able to
support a B+-tree node in 512B (64$\times$8). For a B+-tree node in 2KB with
32 cache lines, a Sentinel Array with four cache lines suffices. In the worst
case for searching over such a 2KB node, all four cache lines of Sentinel
Array and one cache line of KV pair are scanned. So in all, there are five
cache misses. However, in the worst case of searching the original 2KB node
incurs 32 cache misses. Concretely, employing a Sentinel Array significantly
reduces cache misses.

We note that Sentinel Array is not enforced consistency with cache line flush
and memory fence. Then reason is that Sentinel Array could be reconstructed
from keys in the original B+-tree node if the system crashes down. State-of-
the-art B+-trees developed for NVM embrace comprehensive mechanisms to
guarantee the consistency of KV pairs.

###  III-A Search of Sentinel Array

Figure 3: An illustration of searching 140 in a B+-tree node

(a) Average latency (1m data, search, $\mu$s)

(b) Average latency (1m data, insertion, $\mu$s)

(c) Average latency (10m data, search, $\mu$s)

(d) Average latency (10m data, insertion, $\mu$s)

Figure 4: A Comparison of Six Trees on Inserting and Searching 1/10 Million
Keys

(a) Average latency (100m data, search, $\mu$s)

(b) Average latency (100m data, insertion, $\mu$s)

Figure 5: A Comparison of Six Trees on Inserting and Searching 100 Million
Keys

Fig. 3 exemplifies the process of searching with Sentinel Array. The target
key is 140. We assume that four KV pairs occupy one cache line for 16B KV
pointer size in total and 64B cache line’s size. As shown by the left part of
Fig. 4, to reach 140, we need to traverse four cache lines, thereby causing
four cache misses. With a Sentinel Array shown on the right part of Fig. 4, we
only need to check the cache line of sentinels and the actual cache line
holding key 140 and its corresponding value, i.e., two cache misses. Algorithm
1 illustrates the algorithmic steps of searching with Sentinel Array. We note
that we add Sentinel Arrays to a B+-tree’s leaf nodes without loss of
generality.

Algorithm 1 Search With Sentinel Array

0: The target key.

0: The pointer of value of target key.

1: $begin\\_index=0$

2: for ($i=1$; $i<\frac{n}{count\\_in\\_line}$; $i++$) do

3: $//$ $count\\_in\\_line$ is number of KV pairs counted in one cache line

4: if $key<sentinel\\_array[i]$ then

5: break

6: end if

7: $begin\\_index=begin\\_index+count\\_in\\_line$

8: end for

9: Search the node’s array begin with $begin\\_index$

###  III-B Update of Sentinel Array

From time to time, keys are inserted and deleted in a B+-tree node. Therefore,
we need to update the node’s Sentinel Array, which is simply updated with the
original data array. Because the Sentinel Array is reconstructable from
consistent KV pairs, the update of it is efficient without using cache line
flush or memory fence.

###  III-C Multi-threading Access

Sentinel array is added to a B+-tree node as an auxiliary component. A node-
level lock shared with the original B+-tree node helps to support multi-
threading write/read with the Sentinel Array.

##  IV Evaluation

In this section, we implement and evaluate our Sentinel Array based on FAST-
FAIR and Circle-Tree in terms of insertion and search performance. We use a
synthetic benchmark YCSB [21] to evaluate the performance between original
FAST-FAIR/Circle-Tree and ones with Sentinel Array.

###  IV-A Evaluation Setup

All of our experiments are conducted on Linux Server (Kernel version 3.10)
with an Intel ® Xeon ® E5-2620v4 2.10GHz CPU with 512KB/2MB/20MB L1/L2/L3
cache, 8GB DRAM. We use DRAM to emulate the NVM space. We keep the read
latency of NVM as the same as DRAM, and emulate the write latency by adding an
extra delay after each clflushshopt instruction. We set the default write
latency of NVM as 300ns. [9, 15, 11]

(a) Average latency (1m data, search, $\mu$s)

(b) Average latency (1m data, insertion, $\mu$s)

Figure 6: A Comparison of Six Trees with Multi-Threading

(a) 99th percentile latency (insertion, $\mu$s)

(b) 99th percentile latency (search, $\mu$s)

(c) 99th percentile latency (update, $\mu$s)

Figure 7: A Comparison of six KV Stores (4KB Node) on SessionStore Workload of
YCSB

(a) 99th percentile latency (insertion, $\mu$s)

(b) 99th percentile latency (search, $\mu$s)

(c) 99th percentile latency (update, $\mu$s)

Figure 8: A Comparison of six KV Stores (4KB Node) with Multi-threading on
SessionStore Workload of YCSB

Fig. 4 and Fig. 5 show the insertion and search performance of multi-version
B+-tree variants. FAST-FAIR has open-source code and Circle-Tree is
implemented with respect to its original literature. Given the large gap
between the performance of the search function mentioned in the Circle-Tree
and the linear search in FAST-FAIR. The Circle-Tree’s search function was
implemented as a linear search starting from the logical base location of a
node to the last valid KV pair. Moreover, the name that has the ‘S’ suffix is
the Sentinel Array. The one that has the ‘F’ suffix contains a fingerprint
array proposed in the FP-Tree. All implementations have been compiled with -O
option.

Fig. 4 and Fig. 5 show average latencies by running 1/10/100 million insertion
and search operations. A shorter average latency means higher performance. We
have chosen the geometric mean to calculate the average latency.

From Fig. 4(a), a B+-tree with Sentinel Array outperforms the one without
Sentinel Array with much shorter search latencies. For example, using the idea
of Sentinel Array helps to reduce the average search latency of Circle-Tree
and FAST-FAIR by 48.4% and 42.6%, respectively, with 4KB node size. However,
as also shown in Fig. 4(a), the effect of using fingerprints is limited. For
example, the search latency of Circle-Tree with 4KB node size is shortened by
3.1%.

Although the Sentinel Array influences the insertion performance because of
its update needs, still with 4KB node, the average insertion latencies of
Circle-Tree and FAST-FAIR with the Sentinel Array are increased by 6.5% and
4.0% that of the same structure without the Sentinel Array in Fig. 4(b).

Another observation illustrated in Fig. 4 is that the search performance of
the B+-tree with the Sentinel Array gets better with larger node size. Search
latency of B+-tree with Sentinel Array is larger than the ones without
Sentinel Array with 512B node. The reason is that CPU combined with the
compiler has branch predict feature and could prefetch the data into cache-
line. The node with small size node could find the target with less cache
misses, but the Sentinel Array implementation needs at least one solid cache
miss to search Sentinel Array and more branch judgment between Array switch.
So it’s better to add the Sentinel Array to the B+-tree like data structure
with large node size. The same results are also observed in Fig. 5, which
refer to more data insertion and search.

###  IV-B Multi-threading Performance

We have performed multi-threading tests with 1/2/4/8 threads, to evaluate the
performance of the B+-tree with Sentinel Array. We evaluated the tests with
4KB node size and the data size is 1 million keys. Fig. 6(a) shows the
performance of each multi-threading search operation and we could observe that
FAST-FAIR and Circle-Tree with the Sentinel Array outperform the one without
the Sentinel Array. For example, using sentinels reduces the average search
latency by 46.1% and 47.6%, for FAST-FAIR and Circle-Tree with 2 threads.

The multi-thread insertion showed in Fig. 6(b), indicates that the Sentinel
Array has insignificant influence on insertion performance.

###  IV-C Evaluation on YCSB

For end-to-end comparison, we built an interface for each B+-tree like data
structure to receive and handle access requests issued by YCSB. Because the
key from a YCSB workload is a string with a prefix and a number, we removed
the prefix and treated the number as an 8B key in unsigned integer. The
default load value of YCSB has ten fields with 100 bytes per field, so the
pointer of a KV pair in a leaf node points to the address of stored value.

We used the ‘workloada’ workload with YCSB. It first inserts a predefined
number of KV pairs (load) and then follows a search/update ratio of 50%/50%
over keys selected in accordance with a Zipf distribution (run). YCSB reports
a series of latencies from which we choose the 99th percentile latency to rule
out the impact of very short or very long operations. It means that 99% of
overall write/read requests can be completed below such a latency.

In our experiment, we use 1 million KV pairs to test performance with YCSB.
Fig. 7(a) shows the average latencies of inserting 1 million data into each
tree. The structure with the Sentinel Array still gets lower latency than the
original implementation. Fig. 7(b) captures search performance and FAST-FAIR
and Circle-Tree with the Sentinel Array get 16.9% and 13.0% improvement on 4KB
node size. From Fig. 7(c), which represents the average latency of updating
the KV pair of the loaded data, the structure with the Sentinel Array still
gets a shorter latency. The reason is that update operation is like the search
operation and they only need to find the KV pair’s position without shifting
the KV pair.

We also tested multi-threads of the B+-tree with Sentinel Array with YCSB. As
showed in Fig. 8(a), Fig. 8(b) and Fig. 8(c) the average search latencies is
significantly reduced compared to ones without Sentinel Array.

##  V Conclusion

Search is of primary importance for B+-tree. In this paper, we revisited the
search operation of in-NVM B+-tree with sorted nodes and found that the search
performance is mainly affected by read amplification and cache misses. As a
result, we proposed to use a sentinel to stand for each cache line of a sorted
node. Before checking the actual array of KV pairs, a search operation first
goes through the probing Sentinel Array to determine which cache line the
target key would reside in. As the Sentinel Array is much smaller than the
actual array of KV pairs, the number of cache misses is significantly reduced.
Extensive experiments confirm that using sentinels significantly boosts the
search performance of state-of-the-art B+-trees developed for NVM, especially
with larger nodes.

## References

  * [1] L. M. Grupp, J. D. Davis, and S. Swanson, “The bleak future of NAND flash memory,” in _Proceedings of the 10th USENIX Conference on File and Storage Technologies_ , ser. FAST’12. USA: USENIX Association, 2012, p. 2.
  * [2] B. C. Lee, E. Ipek, O. Mutlu, and D. Burger, “Architecting phase change memory as a scalable dram alternative,” vol. 37, no. 3. New York, NY, USA: Association for Computing Machinery, Jun. 2009, p. 2–13. [Online]. Available: https://doi.org/10.1145/1555815.1555758
  * [3] M. K. Qureshi, V. Srinivasan, and J. A. Rivers, “Scalable high performance main memory system using phase-change memory technology,” in _Proceedings of the 36th Annual International Symposium on Computer Architecture_ , ser. ISCA ’09. New York, NY, USA: Association for Computing Machinery, 2009, p. 24–33. [Online]. Available: https://doi.org/10.1145/1555754.1555760
  * [4] P. Zhou, B. Zhao, J. Yang, and Y. Zhang, “A durable and energy efficient main memory using phase change memory technology,” p. 14–23, 2009. [Online]. Available: https://doi.org/10.1145/1555754.1555759
  * [5] S. Raoux, G. W. Burr, M. J. Breitwisch, C. T. Rettner, Y.-C. Chen, R. M. Shelby, M. Salinga, D. Krebs, S.-H. Chen, H.-L. Lung _et al._ , “Phase-change random access memory: A scalable technology,” _IBM Journal of Research and Development_ , vol. 52, no. 4.5, pp. 465–479, 2008.
  * [6] T. Kawahara, “Scalable spin-transfer torque RAM technology for normally-off computing,” _IEEE Design Test of Computers_ , vol. 28, no. 1, pp. 52–63, 2011.
  * [7] T.-C. Chang, K.-C. Chang, T.-M. Tsai, T.-J. Chu, and S. M. Sze, “Resistance random access memory,” _Materials Today_ , vol. 19, no. 5, pp. 254–264, 2016\.
  * [8] R. Bayer and E. M. Mccreight, “Organization and maintenance of large ordered indexes.” Berlin, Heidelberg: Springer-Verlag, Sep. 1972, vol. 1, no. 3, p. 173–189. [Online]. Available: https://doi.org/10.1007/BF00288683
  * [9] J. Yang, Q. Wei, C. Chen, C. Wang, K. L. Yong, and B. He, “NV-Tree: Reducing consistency cost for NVM-based single level systems,” in _13th USENIX Conference on File and Storage Technologies (FAST 15)_. Santa Clara, CA: USENIX Association, Feb. 2015, pp. 167–181. [Online]. Available: https://www.usenix.org/conference/fast15/technical-sessions/presentation/yang
  * [10] D. Hwang, W.-H. Kim, Y. Won, and B. Nam, “Endurable transient inconsistency in byte-addressable persistent B+-Tree,” in _16th USENIX Conference on File and Storage Technologies (FAST 18)_. Oakland, CA: USENIX Association, Feb. 2018, pp. 187–200. [Online]. Available: https://www.usenix.org/conference/fast18/presentation/hwang
  * [11] C. Wang, G. Brihadiswarn, X. Jiang, and S. Chattopadhyay, “Circ-tree: A B+-Tree variant with circular design for persistent memory,” _arXiv preprint arXiv:1912.09783_ , 2019.
  * [12] I. Moraru, D. G. Andersen, M. Kaminsky, N. Tolia, P. Ranganathan, and N. Binkert, “Consistent, durable, and safe memory management for byte-addressable non volatile main memory,” in _Proceedings of the First ACM SIGOPS Conference on Timely Results in Operating Systems_ , ser. TRIOS ’13. New York, NY, USA: Association for Computing Machinery, 2013. [Online]. Available: https://doi.org/10.1145/2524211.2524216
  * [13] S. Venkataraman, N. Tolia, P. Ranganathan, and R. H. Campbell, “Consistent and durable data structures for non-volatile byte-addressable memory,” in _Proceedings of the 9th USENIX Conference on File and Stroage Technologies_ , ser. FAST’11. USA: USENIX Association, 2011, p. 5.
  * [14] I. Oukid, J. Lasperas, A. Nica, T. Willhalm, and W. Lehner, “FPTree: A hybrid SCM-DRAM persistent and concurrent B-tree for storage class memory,” in _Proceedings of the 2016 International Conference on Management of Data_ , 2016, pp. 371–386.
  * [15] J. Zhao, S. Li, D. H. Yoon, Y. Xie, and N. P. Jouppi, “Kiln: Closing the performance gap between systems with and without persistence support,” in _Proceedings of the 46th Annual IEEE/ACM International Symposium on Microarchitecture_ , ser. MICRO-46. New York, NY, USA: Association for Computing Machinery, 2013, p. 421–432. [Online]. Available: https://doi.org/10.1145/2540708.2540744
  * [16] C. Mohan, D. Haderle, B. Lindsay, H. Pirahesh, and P. Schwarz, “ARIES: a transaction recovery method supporting fine-granularity locking and partial rollbacks using write-ahead logging,” _ACM Transactions on Database Systems (TODS)_ , vol. 17, no. 1, pp. 94–162, 1992.
  * [17] S. R. Dulloor, S. Kumar, A. Keshavamurthy, P. Lantz, D. Reddy, R. Sankaran, and J. Jackson, “System software for persistent memory,” in _Proceedings of the Ninth European Conference on Computer Systems_ , ser. EuroSys ’14. New York, NY, USA: Association for Computing Machinery, 2014. [Online]. Available: https://doi.org/10.1145/2592798.2592814
  * [18] Intel, “Intel® 64 and IA-32 architectures software developer’s manual,” _Volume 3A: System Programming Guide, Part_ , vol. 1, no. 64, p. 64, Sep. 2016.
  * [19] Y. Lu, J. Shu, L. Sun, and O. Mutlu, “Loose-ordering consistency for persistent memory,” in _2014 IEEE 32nd International Conference on Computer Design (ICCD)_ , 2014, pp. 216–223.
  * [20] A. Danowitz, K. Kelley, J. Mao, J. P. Stevenson, and M. Horowitz, “CPU DB: recording microprocessor history,” _Queue_ , vol. 10, no. 4, pp. 10–27, 2012\.
  * [21] B. F. Cooper, A. Silberstein, E. Tam, R. Ramakrishnan, and R. Sears, “Benchmarking cloud serving systems with YCSB,” in _Proceedings of the 1st ACM symposium on Cloud computing_ , 2010, pp. 143–154.


ISVHECRI 2018

11institutetext:  Institute of Nuclear Physics, Polish Acedemy of Sciences
Krakow, Poland
11email<EMAIL_ADDRESS>22institutetext: website:
http://credo.science

# Detection of Cosmic-Ray Ensembles with CREDO

Krzysztof W. Woźniak for the CREDO Collaboration The CREDO Collaboration: K.
Almeida-Cheminant, Ł. Bratek, D.A. Castillo, N. Dhital, A.R. Duffy, D. Góra,
B. Hnatyk, P. Homola, P. Jagoda, J. Jałocha-Bratek, M. Kasztelan, D. Lemański,
P. Kovács, V. Nazari, M. Niedźwiecki, K. Smelcerz, K. Smolek, J. Stasielak, S.
Stuglik, O. Sushchov, K.W. Woźniak, J. Zamora-Saá. 1122

###### Abstract

One of the main objectives of cosmic-ray studies are precise measurements of
energy and chemical composition of particles with extreme energies. Large and
sophisticated detectors are used to find events seen as showers starting in
the Earth’s atmosphere with recorded energies larger than 100 EeV. However, a
Cosmic-Ray Ensemble (CRE) developing before reaching the Earth as a bunch of
correlated particles may spread over larger areas and requires an extended set
of detectors to be discovered. The Cosmic-Ray Extremely Distributed
Observatory (CREDO) is a solution to find such phenomena. Even simple
detectors measuring the particle arrival time only are useful in this
approach, as they are sufficient both to provide candidate CRE events and to
determine the direction from which they are arriving.

##  1 Introduction

The primary ultra-high energy cosmic-ray (UHECR) particles arriving at the
Earth are very rare as their number falls steeply with increasing energy.
Those which interact strongly or electromagnetically do not reach the surface
of the Earth, but instead initiate a cascade of secondary particles in the
atmosphere. These secondary particles reach the ground as a shower of
relatively low energy muons, electrons and photons and are detectable in
cosmic-ray observatories. The most extended showers from a single cosmic-ray
particle may, at the ground level, have a diameter of 3-4 km uhecr_review .
The largest cosmic-ray observatories like Auger auger  and Telescope Array
telescope_array  have their detectors spread over very larger areas (3000 km2
and 700 km2, respectively) to ensure that such showers are fully contained.
This coverage is however not sufficient in case the shower starts not at the
top of the Earth atmosphere but at a significantly larger distance SPS_Erber ;
SPS_McBreen . Correlated bunches of particles may be initiated by UHECR
photons and decays of hypothetical heavy and energetic Dark Matter particles.

Figure 1: Detection of the Cosmic-Ray Ensemble (CRE) in distributed cosmic-ray
detectors.

Simulations of Super pre-showers (SPS) initiated by a photon and starting near
the Sun show that in such a scenario the secondary photons near the Earth
would form a disk with a diameter of hundreds of km but with thickness of a
few meters only dithal_icrc2017 . In this case many photons arrive at almost
the same time on top of the Earth atmosphere. Each of them has only a fraction
of the energy of the primary photon, but still they may form numerous
individual cascades. The energy of these cascades is at least an order of
magnitude lower than that of the primary particle and the distance between
them may be a few km or more. Detection of such small showers, even if
strongly correlated in time, by the traditional observatories like Auger may
be impossible. If in the pre-shower only a few particles are created it is
most probable that in the area covered by the observatory only one shower
appears, and it would not differ from showers generated by a single cosmic
ray. Several small showers may develop in the area covered by the observatory,
if in the pre-shower the number of particles is much larger. However, in this
case the energy of these particles has to be relatively small, thus the
signals in the detectors may have a pattern too different from a typical high
energy cascade to trigger recording of the event.

Figure 2: A smartphone with the CREDO application.

Figure 3: Examples of images taken with the completely covered camera of a
smartphone using the CREDO application. Figure 4: Estimation in two dimensions
of the direction of an incoming Cosmic-Ray Ensemble on the basis of position
and time registration in two detectors.

##  2 Cosmic-Ray Extremely Distributed Observatory

Particles from a Cosmic-Ray Ensemble may arrive in places hundreds or even
thousands km distant and produce many individual cascades. Their signature
would be the coincidence in time. Consequently, the detectors used to register
them need to cover the largest area possible, but may be quite sparse. It is
not necessary to register many signals from individual cascades, it’s better
to detect signals from several different cascades, even if it would be a
single measurement for each cascade. The best results can be achieved when the
information from detectors spread over a country, continent or even the whole
Earth are collected and analyzed together. The idea of CRE detection is
presented in Figure 1. The individual cascade leaves signal in one or a few
detector stations. There are stations in which no signal is registered in the
analyzed time interval and also the stations with a background signal - for
example cosmic rays not connected with CRE or from decays of radioactive
nuclei.

The detectors registering the arrival time of particles are sufficient for
testing the coincidences. However, rejection of backgrounds play an important
role in extraction of the signal and measurements of CRE properties described
in Section 3. Assuming that CRE initiates several small cascades it’s
reasonable to require a hardware based local coincidences of individual
detectors within a station or a coincidence of nearby stations. This way
contribution from uncorrelated background signals (mainly radioactive decays)
can be suppressed. It does not solve the problem of small cascades from cosmic
rays not belonging to the Cosmic-Ray Ensemble - but their contribution, at
least partially, can be removed in the analysis.

The minimal hardware requirements are so easy to fulfill that the simplest of
detectors capable of registering the time of a passage of a particle from a
cosmic shower can be useful. The most popular of such devices is a smartphone
equipped with a camera and a GPS module. It can work as a cosmic-ray detector
if the camera lens is covered (to prevent external light) and an appropriate
application is running.

The idea of using smartphones to detect cosmic rays is not new. There are two
groups CRAYFIS crayfis_web ; crayfis_arxiv  and DECO deco_web ; deco_icrc2015
which developed their custom applications for this purpose for several years.
However, their software is proprietary thus the Open Source CREDO application
(Figure 2) was also created credo_app .

The CREDO application works when the smartphone is attached to a power supply
(to prevent discharging of the battery) and the camera is protected from
light. The pictures similar to those presented in Figure 3 are constantly
taken and those with a pattern ressembling that from a cosmic ray are sent to
the server over the Internet. The information from the smartphones is stored
and a summary of the activities is provided in a WWW page credo_server .
Currently, the work on prototypes of simple and inexpensive desktop type
cosmic-ray detectors is also ongoing.

##  3 Time Based Triangulation

The Cosmic-Ray Ensemble created at a point very far from the Earth arrives as
a front of almost parallel moving particles. In this case simple measurement
of the position of the detector and the time when such particle arrives is
sufficient to determine the direction of the front, as shown in Figure 4.
Assuming the detected particles are traveling with the speed of light, $c$,
the registration of the times $t_{1}$ and $t_{2}$ in two detectors which are
placed at the distance $L$, enables calculation of the inclination of the
front, $\gamma$, as:

$cos(\gamma)=c(t_{1}-t_{2})/L.$ (1)

This angle calculated in two dimensions is in fact a half of the apex angle of
the cone that is formed by all possible directions of the front. Two detectors
(stations) are not sufficient to determine the direction of the front in 3D.
However it is enough to add a third station with a signal, to construct two
cones of possible directions of the front, with a common apex. Then there are
two lines - intersections of the cones - and that of them which is pointing
above the ground is the correct one. This way even only three station with
signals from CRE with a flat front after relatively simple calculations of the
Time Based Triangulation (TBT) method allow the estimation of the direction of
the CRE. However, in the presence of ia background signal in even one of these
stations the results of calculations are severely distorted. Frequently, the
estimated cones may have no intersection, either because a small cone is
inside a bigger one or the two small cones are separated. In most cases,
however, TBT calculations give a direction differing greatly from the correct
one. It is thus necessary to test all combinations of three stations and
analyze calculated directions to find the CRE direction.

Figure 5: Estimation of possible directions of an incoming Cosmic-Ray Ensemble
using the Time Based Triangulation method. Simulated signals from CRE were
detected in 5 stations while in another 10 stations randomly generated
background signals are present. In the upper plot all estimated directions in
geographic coordinates system are shown in Mollweide projection. The correct
CRE direction is denoted by the green circle, red points represent false
estimates, obtained from calculations using at least one background signal. In
the lower plot for selected range of latitude, $\phi$, and longitude,
$\lambda$, the calculated directions are shown as boxes, the correct ones are
in the green circle.

##  4 Simple CRE simulation

The accuracy of calculations in the Time Based Triangulation method was tested
using simple simulations of a Super-preshower starting very far from the Earth
and thus approximated by a flat front of particles. The positions of
hypothetical detectors are randomly generated in the area approximately
representing the territory of USA. The exact arrival time is randomly smeared
with a width of 1/3000000 s, an uncertainty coresponding to the 1 km position
error along the SPS direcion. Signals from the SPS is simulated always in 5
stations, but in addition in some number of other stations background signals
are randomly generated in the time interval around the arrival time of the
SPS. In Figure 5 results of TBT reconstruction for one event in which SPS
signals and background are present in 5 and 10 stations, respectively, are
shown. The calculations were performed for all 15!/(3!(15-3)!) = 455
combinations of three stations with signals, but only in 5!/(3!(5-3)!) = 10
combinations in all selected stations SPS signal is recorded. Even if one
station contains signal from a background the result of TBT calculations can
be completly wrong. In the grid of geographical coordinates (Figure 5 upper
plot) many calculated directions are thus false and are spread over almost a
half of the solid angle. In a closer examination of the angular area more
focused on the real direction of SPS we can see an evident cluster of results
near the correct direction and a less dense population around it. The
calculated values from real signals are quite precise, the errors do not
exceed 1∘ thus for small and moderate background fractions the cluster at the
SPS directions should be easy to identify.

In the presented in Figure 5 case of 5 real and 10 background signals the
correct SPS direction can be guessed after checking the plot "by eye", but in
more difficult cases much more elaborate data analysis is necessary. In
simulations with smaller number of stations with a background the density of
false directions is lower and the correct direction is more pronounced. For
larger number of stations with backgrounds the cluster showing the correct
direction becomes less visible and may be difficult to find. It is thus very
important to keep the ratio of stations with real signal to those with a
background signal as high as possible.

##  5 Summary

An obvious extension to the measurements of cosmic-ray showeres initiated by a
single high energy particle are studies of Cosmic-Ray Ensembles. They may give
an opportunity to discover new phenomena including new types of interactions
or decays of massive particles. The purpose of CREDO is to collect all
available cosmic-ray data and analyze correlations of signals in detectors
spread over large distances. In this case detectors as simple as smartphones
may be useful, thus an Open-source CREDO application was created. For
measurements of electromagnetic Super-preshowers information from even small
tabletop detectors may be sufficient. The Time Based Triangulation method
allows calculation of direction of a SPS if only the background rate is not
too large.

## Acknowledgments

This research has been supported in part by the Governments of Czechia,
Hungary, Poland and Slovakia through Visegrad Grant 21720040 from
International Visegrad Fund and by the PLGrid Infrastructure. We thank the
staff at ACC Cyfronet AGH-UST for their helpful support. CREDO application is
developed in Cracow University of Technology.

## References

  * (1) K. Kampert, A.A. Watson, arXiv:1207.4827 [physics.hist-ph] (2012).
  * (2) The Pierre Auger Collaboration, Nucl. Instr. Meth. A 798 (2015) 172.
  * (3) T. Abu-Zayyad, et al. Nucl. Instr. Meth. A 689 (2012) 87.
  * (4) T. Erber, Rev. Mod. Phys. 38 (1966) 626.
  * (5) B. McBreen and C. J. Lambert, Phys. Rev. D24 (1981) 2536.
  * (6) N. Dhital, et al. (CREDO Collaboration), PoS (ICRC2017) 1078 (arXiv:1709.05196 [astro-ph.IM]).
  * (7) https://crayfis.io
  * (8) D. Whiteson, M. Mulhearn, C. Shimmin, K. Cranmer, K. Brodie, D. Burns (CRAYFIS), arXiv:1410.2895 [astro-ph.IM].
  * (9) https://wipac.wisc.edu/deco/home
  * (10) J. Vandenbroucke et al. (DECO), PoS (ICRC2015) 691 (arXiv:1510.07665v1 [astro-ph.IM]).
  * (11) https://credo.science/credo-detector-mobile-app/
  * (12) https://api.credo.science/web/