# Coherence properties of NV-center ensembles in diamond coupled to an
electron-spin bath

Reyhaneh Ghassemizadeh  Fraunhofer Institute for Mechanics of Materials IWM,
Wöhlerstraße 11, 79108 Freiburg, Germany  Wolfgang Körner  Fraunhofer
Institute for Mechanics of Materials IWM, Wöhlerstraße 11, 79108 Freiburg,
Germany  Daniel F. Urban<EMAIL_ADDRESS>Fraunhofer Institute
for Mechanics of Materials IWM, Wöhlerstraße 11, 79108 Freiburg, Germany
Christian Elsässer  Fraunhofer Institute for Mechanics of Materials IWM,
Wöhlerstraße 11, 79108 Freiburg, Germany  University of Freiburg, Freiburg
Materials Research Center (FMF), Stefan-Meier-Straße 21, 79104 Freiburg,
Germany

###### Abstract

We investigate nitrogen-vacancy center (NV) ensembles in diamond under the
influence of strongly-correlated electron-spin baths. We thoroughly calculate
the decoherence properties of the NV central spin for bath concentrations of
$0.1-100$ ppm using the cluster-correlation expansion (CCE) method. We
systematically analyze possible origins of the significant deviations in the
values of the $T_{2}$ coherence time reported in literature. We demonstrate
that significant variations can originate from the choice of averaging and
fitting procedures used for the ensemble average and we point out the
respective aspects that need to be considered, when comparing the various
theoretical studies. Our study may ease readers to perform reliable and fast
simulations on the central spin problem. It provides an understanding and
interpretation of the outcome parameters describing the dynamics of the local
bath spins.

##  I Introduction

In the past two decades, solid-state based platforms for quantum technology
applications like quantum sensing and quantum computing became increasingly
popular and promising due to their outstanding properties like, e.g., their
long spin coherency even at room temperature [1, 2, 3, 4, 5, 6, 7, 8]. The
negatively charged nitrogen vacancy (NV) defect complex in diamond, which
consists of a substitutional nitrogen atom with a vacant nearest-neighbor
carbon site and an additional electron, is one of the candidates for the near
future qubits due to its extraordinary spin properties [9, 10]. The electronic
ground state of the NV center is a triplet state, $S=1$, with a spin life-time
of about 1 ms at room temperature [11, 12]. However, this long spin coherency
can be limited by the interaction of the NV electron spin with surrounding
spins and crystal defects. In highly purified diamond crystals, the main
source of the NV-spin decoherency is known to be nuclear spins of surrounding
13C isotopes [13, 14]. This decoherence mechanism is well understood within
the quantum mechanical framework [15, 12, 17, 16, 18, 19, 20]. However, during
the growth of diamond by chemical vapor deposition (CVD) or ion implementation
for creating the NV centers, paramagnetic impurities, e.g., substitutional 14N
(so-called P1 centers), neutral NV0, NVH-, or carbon-vacancy pairs VV0, are
frequently formed [21, 22, 24, 23, 25, 35]. These defects have unpaired
electron spins, which can strongly couple to the NV center and play a crucial
role in decreasing the spin coherency and, hence, the efficiency of the NV-
based devices. Consequently, many experimental [3, 27, 21, 34, 25, 22, 14, 26,
28, 35] and theoretical [31, 30, 14, 29, 32, 50] studies investigated the
decoherence of a central spin in a strongly coupled spin environment in the
past decade.

Generally, the normalized coherence decay signal $\mathcal{L}(t)$ of a central
spin in an interacting spin bath is well described by classical stochastic
noise, e.g., the Ornstein-Uhlenbeck process [33, 42, 34]. It can be fitted by
an exponential decay function

$\mathcal{L}_{\rm fit}(t)=\exp\left[-\left(\frac{t}{T_{2}}\right)^{p}\right],$
(1)

with the stretched exponential parameter $p$ and the coherence time $T_{2}$ as
fit parameters [33, 34, 14]. These contain important information on the NV
center spin decoherency. Here, $p$ determines the curvature and shape of the
decay and $T_{2}$ quantifies the relevant time scale for the decay. The spin-
echo decoherence signal of a single NV center in a weakly coupled bath of
nuclear spins is given by an approximately Gaussian-like decay with $p\sim
2$–3 [33, 42, 34, 17]. On the other hand, for an ensemble of NV centers in a
strongly correlated electron spin bath, the decay signal can be well
approximated by Eq. (1) with $p\approx$ 1–1.5 [36, 29, 14, 35, 22]. Recently
it has been reported, that $p$ can even take values below 1 in the case of a
quasi two-dimensional spin environment [37, 32]. Moreover, for the case of
strongly coupled P1 electron spins, an almost linear dependency of $T_{2}$ on
the inverse of the concentration of P1 centers ($\rho_{P1}$) has been reported
from both, theoretical calculations [29, 14, 30, 32, 50] and experimental
measurements [27, 14]. However, the reported $p$ and $T_{2}$ values vary over
a wide range, as demonstrated by the compilation in Tab. 1.

The broad range of $T_{2}$ and $p$ values obtained in different theoretical
studies motivates us to analyze the possible sources of this variety.
Therefore, we thoroughly analyze possible origins for the large variation of
50% and 100% in $T_{2}$ and $p$ values, by considering different analytical
and numerical criteria. We use the established CCE method to calculate the
Hahn-echo [19] decoherence signal of ensembles of NV centers in various
concentrations of the electron-spin bath ($\rho_{e}=0.1-100$ ppm). We chose
the pure electron bath environment because the existence of defects with an
extra electron have as an important impact on the NV center decoherence as the
P1 centers [29], and there is still a lack of studies in this regard in the
literature. Moreover, the conceptional sources of discrepancy in the $T_{2}$
and $p$ values are independent of the addition of hyperfine-field (HF)
coupling between nuclear and electronic bath spins.

By varying the spatial configuration of the spins surrounding the NV center we
obtained a wide range of coherence decay behaviors for individual bath
configurations, from an extremely fast decay to a slow decay. This leads to a
wide range of values of T2 and $p$ for each individual spatial spin
configuration at given $\rho_{e}$ with a large standard deviation. Hence,
including enough realizations in the ensemble average is crucial to obtain a
robust and reliable estimate of the coherence time $T_{2}$. Moreover, the
convergence parameters of the CCE method and the fit uncertainties in
numerical fits of the decay function of Eq 1 are further origins of the
inconsistency in the $T_{2}$ and $p$ values even at the same level of theory.

The manuscript is organized as follows. As a starting point, in Sec. II we
briefly review the available literature on theoretical studies of the
decoherence of a NV spin subject to an environmental spin bath. In Sec. III we
outline our theoretical model and the numerical settings used in the CCE
calculations. Our results and the related discussions in Sec. IV are
structured by first analysing the limitations in the ensemble average in Sec.
IV.1 and then discussing the relevance of the order of the CCE approximation
in Sec. IV.2. Our results for the decoherence signal of an ensemble of NV
centers in various concentrations of the electron-spin bath
($\rho_{e}=0.1-100$ ppm) are presented in Sec. IV.3 and augmented by a
detailed analysis of different fitting strategies in Sec. IV.4. We summarize
our findings in Sec. V.

##  II Review of literature results

| 1 [ppm] |  | 10 [ppm] |  | 100 [ppm]
---|---|---|---|---|---
|  $T_{2}$ [$\mu$s] | $p$ |  |  $T_{2}$ [$\mu$s] | $p$ |  |  $T_{2}$ [$\mu$s] | $p$
Ref.[32]\- JT(r) | 225 | 1.03 |  | - | - |  | - | -
Ref.[50]\- JT(r) | 230 | 1.2-1.3 |  | 25 | 1.2-1.3 |  | 3.5 | 1.2-1.3
Ref.[29]\- JT(r) | 407 | 0.8 |  | 39 | 0.9 |  | 3.3 | 0.8
Ref.[14]\- JT(r) | 70-140 | 1.5 |  | 8-10 | 1.5 |  | 0.6-0.8 | 1.5
Ref.[29]\- JT(a) | 250 | 0.9 |  | 25 | 1 |  | 2.5 | 0.8
Ref.[32]-no HF | 49.2 | 1.43 |  | - | - |  | - | -
Ref.[29]-no HF | 91.2 | 1.1 |  | 8.9 | 1.1 |  | 0.87 | 1.1
Ref.[30]∗ | 20 | 2-3 |  | 2.0 | 2-3 |  | 0.2 | 2-3
Table 1: Theoretical values of $T_{2}$ and $p$ for the coherence decay of the
NV ensemble in 1, 10 and 100 ppm concentration of P1 bath spins, as reported
in literature. Methods include random or aligned Jahn-Teller axes, denoted as
JT(r) and JT(a), respectively. ∗No changes in values with or without HF are
reported in Ref. [30].

Table 1 summarizes the theoretical results for $T_{2}$ and $p$ for NV center
ensembles in different concentrations of P1 centers, including nuclear and
electron spin baths, as reported in literature to the best of our knowledge.
In a P1-center spin-bath the hyperfine interaction (HF) between the electron
and the respective nuclei spin gives rise to an additional term in the
Hamiltonian of the bath which is not present in the case of a pure electron-
spin bath. Moreover, the splitting of the hyperfine levels of the P1 bath is
different based on the orientation of the Jahn-Teller (JT) axes of the P1
centers [29].

In an earlier simulation of Wang et al. [30] a so-called $P$-representation
method [38] was used to calculate the decoherence of the NV spin coupled to
the surrounding P1 spins. Although the derived linear dependency of the
$T_{2}$ time on the inverse of the P1 center concentration agrees to the more
recent available studies, the values of $T_{2}$ are highly underestimated
compared to the later works. Moreover, Ref. [30] reports no significant
difference in the $T_{2}$ times with and without considering the HF coupling
of the P1 centers.

The more recent work of Bauch et al. [14] employed a semi-classical approach
to estimate the values of $T_{2}$ as function of the concentration of P1
centers, $\rho_{P1}$. Here, P1 centers are found to be the dominating source
of the NV decoherence for $\rho_{P1}\geq$ 0.5 ppm where
$T_{2}\propto(\rho_{P1})^{-1.07}$. Including HF coupling with a random
distribution of JT axes in their theoretical model, they obtained $T_{2}$
values which are considerably below those obtained later by quantum mechanical
approaches [29, 32, 50]. They argue that their semi-classical approach only
poorly captures the dynamics of the bath electron spins. Furthermore, the
inaccuracy of the employed fitting caused a broad spread of their computed
values.

In a recent work, Park et al. [29] employ the quantum mechanical approach of
cluster correlation expansion (CCE) [19, 39, 40, 12, 42, 41] to accurately
determine the dynamics of strongly correlated spin baths. They considered the
cases where the HF coupling is included with either random or aligned JT axes,
or it is totally excluded, which mimics the behavior of a pure electron-spin
bath. They show that including HF interactions strongly suppresses the
electron-spin flip-flop transitions of the NV center and leads to
approximately 4.5 times larger $T_{2}$ values compared to the experimental
values. The addition of electron spins to the bath with at least the same
concentration as P1 centers decreases their $T_{2}$ values [29] to the
experimental values reported by Bauch et al. [14].

In the very recent studies of Schätzle et al. [32] and Marcks et al. [50] the
influence of the diamond layer thickness on the coherency of the NV ensembles
in the P1 spin baths is investigated using the CCE method. Both studies
indicate an increase in spin coherency for thin films of diamond. Their
reported T2 values for bulk diamond vary within less than 3% in comparison
with each other but they deviate by about 50% from the CCE results of Park et
al. [29].

##  III Theoretical Model

We use the generalized cluster correlation expansion method (gCCE) [43, 44] as
implemented in the python module PyCCE [43] to calculate the coherence
function of the central electron spin of the NV center in a bath of
paramagnetic spins. The CCE method quantum mechanically describes the
decoherence of a central spin (or qubit) in a bath of surrounding spins by
capturing the time evolution of the bath dynamics. In this approach, the
surrounding bath spins are grouped into clusters of $n$ spins (n = 1, 2, etc.)
in terms of the strength of their spin-spin dipolar interaction. The
contribution of each spin cluster is calculated quantum mechanically. The
total coherence function $\mathcal{L}(t)$ of the central spin is factorized
into sets of irreducible contributions from spin bath clusters [39, 40]. The
validity of the CCE method for solving the problem of the central spin in spin
baths of both nuclei [19, 12, 48, 47, 44, 43, 45, 46] and electrons [49, 31,
28, 29, 32] has been well established. However, in case of more strongly
correlated systems, like electron spin baths, more strict convergence criteria
need to be applied. The important convergence parameters in CCE are: (1) the
CCE order which defines the maximum number of spins interacting within a
cluster, (2) the maximum distance $r_{\rm dipole}$ between two bath spins that
can be assigned to the same cluster, and (3) the cut-off radius $r_{\rm bath}$
which determines the spherical region around the central spin in which bath
spins have an active influence on the central spin decoherency in the
considered timescale.

Our model systems are constructed as follows. We start from a cubic diamond-
crystal unit-cell with a lattice constant of 3.57 Å and expand it to a cubic
supercell with an edge length of approximately 1000 Å. The NV center is placed
at the center of this supercell with the $z$ axis set to be the trigonal
symmetry axis of the NV center along the [111] direction of the diamond. The
electron spins of the bath, with spin $S=1/2$ and gyro-magnetic ratio
$\gamma_{e}$ = -17608.597 rad ms${}^{-1}~{}G^{-1}$, are randomly distributed
at carbon sites in diamond with concentrations of $\rho_{e}=0.1-100$ ppm. An
external magnetic field of $B=100$ G is applied along the $z$ axis. We
examined the changes in coherence time for different values of the magnetic
field. Like Bauch et al. [14] we found that the NV center coherency is
insensitive to the value of the magnetic field if $B>10^{-3}$ G for low
concentrations of $\rho_{e}$ (e.g., 0.1 ppm) and $B>1$ G for higher
concentration (e.g., $\rho_{e}$ = 5 ppm).

The coherency function $\mathcal{L}(t)$ is calculated for the loss of the
relative phases of the qubit levels $m_{s}=0$ and $m_{s}=+1$ in presence of a
Hahn-echo $\pi$-pulse series after a free evolution time $t=2\tau$. Here,
$\tau$ is the time between successive $\pi$ pulses. The contribution of all
the other bath spins outside of a cluster is taken into account as a mean-
field effect via Monte-Carlo sampling over 80-300 random initial spin states
of the bath [44, 43]. Finally, the ensemble-averaged coherence function is
calculated by averaging over 500 spatially different configurations of bath
spins.

Figure 1: Converged values of $r_{\rm dipole}$ and $r_{\rm bath}$ with respect
to $\rho_{e}$ for the 2nd and 3rd cluster orders, abbreviated as gCCE2 and
gCCE3, respectively.

We performed convergence tests for the choice of parameters $r_{\rm dipole}$
and $r_{\rm bath}$ for both, the 2nd and the 3rd order cluster approximations
(gCCE2 and gCCE3). We consider the parameters to be converged when the
variations of the extracted fitting parameters T2 and $p$ are less than 1$\%$
and 5$\%$, respectively. The converged values of $r_{\rm dipole}$ and $r_{\rm
bath}$ as function of the electron-spin density $\rho_{e}$ are displayed in
Fig. 1 and selected values are summarized in Tab. 2. Regarding gCCE2, we
obtain the approximate relation $r_{\rm bath}\approx 1.2r_{\rm dipole}$ for
all the considered bath spin concentrations in agreement with previous works
(see e.g. the SI of Ref. [29]), while in the case of gCCE3, the relation
$r_{\rm bath}\approx 2r_{\rm dipole}$ holds for the different $\rho_{e}$.

| $\rho_{e}$ | 0.1 ppm | 1 ppm | 10 ppm | 100 ppm
---|---|---|---|---|---
gCCE2 |  $r_{\rm dipole}$ [nm] | 615 | 370 | 210 | 105
|  $r_{\rm bath}$ [nm] | 765 | 445 | 250 | 125
|  $\bar{N}_{b}$ ($\sigma_{\bar{N}_{b}})$ | 19 (0.4) | 67 (6) | 117 (10) | 145 (12)
gCCE3 |  $r_{\rm dipole}$ [nm] | 550 | 300 | 170 | 90
|  $r_{\rm bath}$ [nm] | 1000 | 600 | 340 | 180
|  $\bar{N}_{b}$ ($\sigma_{\bar{N}_{b}})$ | 19 | 141 (6) | 290 (16) | 423 (21)
Table 2: The converged values of $r_{\rm dipole}$ and $r_{\rm bath}$, given in
Å, for various bath spin concentrations ($\rho_{e}$) in ppm, the corresponding
averaged number of bath spins, ($\bar{N}_{b}$), and their standard deviation
($\sigma$).

##  IV Results & Discussion

###  IV.1 Ensemble average over spin positions

The decoherence behavior of an individual NV center is strongly dependent on
the specific spatial distribution of the surrounding bath spins [30, 17, 28].
In order to obtain the characteristic time $T_{2}$ describing an ensemble of
independent NV centers, it is therefore essential to average over a
sufficiently large number of bath configurations. In the literature, ensembles
of $80$ to $10^{4}$ individual NV+bath configurations were considered using
different theoretical approaches [30, 14, 29, 50]. In this section, we analyze
the effect of using a finite-size set of configurations and show how many bath
configurations are needed for an average that yields sufficiently small
numerical uncertainties in the extracted $T_{2}$ and $p$ values. To achieve
this, we first compute the decoherence for $N_{\rm max}=500$ random spatial
configurations for each spin bath concentration. Then, we evaluate the
averaged decoherency for ensembles of $N$ = 50 – 250 bath configurations,
picked randomly from the full set of 500 configurations, and repeat this 200
times for each $N$. For each of these ensemble averages we obtain the $T_{2}$
and the exponential power $p$ from a fitted exponential function using Eq.
(1). By this procedure, we obtain a histogram of $T_{2}^{(N)}$ and $p^{(N)}$
values that may be obtained when conducting the ensemble average with $N$
configurations, compared to using the full set of $N_{\rm max}$, at given spin
concentration $\rho_{e}$.

We plot the results of this analysis via the probability distribution function
(PDF) of the relative deviations, $\Delta
T_{2}=\left(T_{2}^{(N)}-T_{2}^{(N_{\rm max})}\right)/T_{2}^{(N_{\rm max})}$
and $\Delta p=\left(p^{(N)}-p^{(N_{\rm max})}\right)/p^{(N_{\rm max})}$, in
Fig. 2. This illustrates the spread of values and the uncertainty that
statistically can be obtained when taking an ensemble average of finite size
$N$. Increasing the number of configurations that enter the ensemble average
narrows down the spread of $T_{2}$ and $p$ values, and we find that an
averaging over 250 samples is sufficient to yield a maximum uncertainty of
approximately $\pm$ 5%.

Figure 2: Relative variation in $T_{2}$ time (top) and $p$ values (bottom)
when averaging over N = 50–250 spatial configurations of the spin bath with
concentration 1 ppm (left), 10 ppm (middle), and 100 ppm (right). Each panel
shows the probability distribution function (PDF) of errors relative to an
ensemble average obtained from 500 configurations (dashed lines).

###  IV.2 Order of the clusters expansion

(a)(b)

Figure 3: Decay of the coherence function $\mathcal{L}(t)$ of a NV center in
diamond subject to an electron spin bath for different orders of the CCE
approximation. The ensemble average is obtained from 500 individual spatial
configurations of electron spins with a concentration of (a) $\rho_{e}=$ 0.1
ppm and (b) $\rho_{e}=$ 5 ppm. Dashed lines indicate different ranges of the
decay as analyzed in Sec. IV.4.

Higher order cluster approximations of CCE are becoming computationally
demanding, especially for higher concentrations of bath spins. For the weakly
interacting spin baths of nuclei the use of the second order cluster
approximation in CCE calculations has been well established [44, 43, 45, 46].
Also in the case of strongly interacting spin baths of electrons it has been
shown that the second order cluster approximation can be sufficient for
isotope-enriched silicon [49] and for NV centers surrounded by very high
concentrations (60-200 ppm) [31] of P1 centers. A detailed study of Witzel et
al. [31, 42] on the 2nd, 3rd and 4rd order cluster approximations for baths of
electron spins indicates that although the second order (CCE2) result for the
averaged decay differs from the third order (CCE3) result, in the range of
short times the decay appears to be captured well by CCE2 approximation. Park
et al. [29] stated that the CCE2 approximation is sufficient for their CCE
simulations. However, the recent study of Marcks et al. [50] uses CCE4 for the
same system.

We performed gCCE calculations for the 2nd, 3rd and 4th order of spin clusters
for low concentrations of 0.1-0.2 ppm. For higher concentrations up to 5 ppm
we restricted our investigation to gCCE2 and gCCE3 due to the high
computational cost of gCCE4 at these high spin concentrations. For individual
bath configurations it may happen that the coherence function $\mathcal{L}(t)$
obtained by the gCCE2 and gCCE3 approximations are numerically almost
indistinguishable. However, in our results for an ensemble average over 250
different spatially distinct bath configurations we always see a difference in
their decay behavior for times of the order of $T_{2}$ or larger, as
illustrated in Fig. 3.

As already pointed out by Witzel et al. [42], the criterion of convergence
depends on which range of the decay time is taken into account (e.g., for
fitting the noise model). On shorter times the results of CCE2 and CCE3 may be
considered as converged. On the other hand, when considering longer times on
the order of $T_{2}$ or more, the CCE2 overestimates the coherence time
roughly by a factor of 1.2 compared to the CCE3 and CCE4. The sensitivity of
the decay curve on the order of the cluster expansion is more pronounced for a
dilute spin bath as can be seen from Fig. 3(a).

In the following, the importance of the choice of higher cluster order to
extract more stable T2 and $p$ values by fitting the analytic function of Eq.
1 to the numerical CCE data is demonstrated by our results.

###  IV.3 Coherence of NV ensembles

We evaluated the coherence decay $\mathcal{L}(t)$ for an ensemble average of
NV center spins obtained by gCCE2 and gCCE3 calculations for electron-spin
baths of $0.1-100$ ppm concentration, and gCCE4 calculations for the most
dilute spin baths of 0.1 – 0.2 ppm. In order to display and compare the
curves, we plot $-\ln(\mathcal{L}(t))$ vs. $t$ in a double-logarithmic graph
in Fig. 4. In the double-logarithmic graph the slope of the curve determines
the (local) value of $p$, since there is not a single value of $p$ which
describes the entire curves, i.e. the full time span [37, 32]. However, as
shown in Fig. 4, using the higher order clusters in the CCE calculations leads
to a nearly perfect linear behavior of the $-ln(\mathcal{L}(t))$ vs. $t$ in
logarithmic scale for the examined time ranges shown.

For the gCCE2 calculations, we observe a decoherence behavior on short times
characterized by a slope of $p\approx 1.5$ for lower concentrations and
$p\approx 1$ for very high concentrations of $\rho_{e}$ = 50 – 100 ppm (c.f.
short black lines in Fig. 4) On the other hand, for the gCCE3 calculations,
the slope of the linear fit is approximately 1.5, regardless of the bath spin
concentrations. The same is true within gCCE4 for $\rho_{e}=0.1$ and 0.2. At
this point we can conclude that the more linear the relationship of numerical
data for $-\ln(\mathcal{L}(t))$ gets, the more accurate a numerical fit using
the function Eq. (1) becomes.

Figure 4: Coherence decay of an NV center in electron-spin baths of
$\rho_{e}=$ 0.1 – 100 ppm plotted in double-logarithmic graphs for three
different cluster orders, gCCE2, gCCE3, and gCCE4. The ensemble averages were
obtained from 500 individual spatial spin configurations at each
concentration. Short solid black lines indicate the slopes for selected $p$
values, for comparison.

(a)(b)(c)(d)(e)(f)

Figure 5: Coherence decay signal $\mathcal{L}(t)$ of a NV center in an
electron bath of concentration (a)-(c) 0.1 ppm and (d)-(f) 1ppm. The different
fit results are shown as gray circles. Left column is an exponential fit to
$\mathcal{L}(t)$. Middle column is a power fit to the $-\ln(\mathcal{L}(t))$
and the right column is the linear fit on the double logarithmic scale. The
red circles indicate the data points included to obtain the fits.

###  IV.4 Extracting T2 and $p$

The values of $T_{2}$ and $p$ obtained for an ensemble average of spin
configurations can be quite sensitive to the choice of fitting procedure,
specifically if the dynamical the interactions are not sufficiently captured
due to the applied approximations in the CCE calculations. In this section we
first discuss the different fitting procedures used so far in literature and
show how sensitively the extracted fitting parameters $T_{2}$ and $p$ are
affected by the choice of the fit function. In the second part, we analyze the
influence of the time range which is taken into account for the fit.

The first step to obtain the $T_{2}$ from a decay signal is to fit a function
which best describes the noisy signal. There are two main approaches in the
literature for extracting $T_{2}$ and $p$ values from the coherence decay
signal $\mathcal{L}(t)$. Either one fits directly the exponential function,
Eq. (1), to the $\mathcal{L}(t)$ vs. $t$ data [49, 30, 28, 29, 14, 44, 35,
50], or one fits a linear function to the $\log(-\ln(\mathcal{L}_{t})))$ vs.
$\log(t)$ data [47, 37, 32]. Another possibility, which was not yet employed
in previous studies, is to fit a power function to the $-\ln(\mathcal{L}_{t})$
vs. $t$ data. In the following, we refer to these three fitting procedures as
_exponential fit_ , _linear fit_ , and _power fit_ , respectively.

In Fig. 5 we illustrate the differences in applying these methods to spin-bath
concentrations of 0.1 and 1 ppm calculated with the gCCE2 approximation.
Although the data points taken into account (red circles) are the same, the
values of $T_{2}$ and $p$ extracted by the three fit methods are different.
The difference between the first two fit methods (the exponential fit and the
power fit) are negligible with respect to the numerical uncertainty of the
CCE. However, the discrepancy between the linear fit and the other two fits is
much more pronounced. For instance, the mean values of the changes in T2 and
$p$ evaluated from an exponential fit compared to a power fit to the data
shown in Fig. 5 are 1.8 % and 2.0 %, respectively. The outcome of the linear
fit leads to a $T_{2}$ that is, on average, approximately 11 % smaller and a
$p$ that is approximately 18 % larger than the results of the exponential fit.
Considering the whole set of bath concentrations studied, we conclude that the
exponential fit allows larger $T_{2}$ and smaller $p$ values while in a linear
fit the T2 shrinks and the $p$ grows, in comparison. We point out, that the
inclusion of very few data points at the shortest times in the linear fit can
notably change its outcome, especially when the set of data is noisy. This is
not the case for the other two fitting procedures. In this regard, the choice
of an exponential fit function is beneficial to avoid the numerical
uncertainties of the noisy data set. However, this sensitivity on the choice
of the fit function becomes less pronounced for the higher order CCE
approximations.

Figure 6: Coherence time T2 obtained from either an exponential (circle
symbols) or a linear (square symbols) fit for two different fitting ranges
$1>\mathcal{L}(t)>\mathcal{L}_{f}$, as function of the electron spin bath
concentration $\rho_{e}$. The inset shows the corresponding results for the
parameter $p$. The black dashed lines mark the values taken from Ref. [29] for
a “pure electron” spin bath.

Next, we examine the range of the data points that need to be included for
fitting the decay function reliably. Given that in the strongly interacting
spin baths the complexity of the spin interactions leads to a more complex
noise signal, the fits outlined are not simply applicable to capture the decay
process over the full time span. In fact, the $T_{2}$ and $p$ values obtained
from fitting to different time intervals can be very different. As already
mentioned above, in order to avoid numerical uncertainties, the decay data
points at the very early times should be excluded in the linear fit.

In order to study the influence of the range of coherence decay (and thereby
the respective time span) which enters the fitting, we define a final value
$\mathcal{L}_{f}$ and truncate the data used for fitting accordingly. Here,
e.g., $\mathcal{L}_{f}=0.6$ means that the coherence decay in the range
$1<\mathcal{L}(t)<0.6$ is included in the fit. We determined the values of
$T_{2}$ and $p$ extracted from the gCCE4 data for densities $\rho_{e}=0.1$ ppm
and 0.2 ppm, and from gCCE3 data for $\rho_{e}$ = 0.5 – 100 ppm using the
exponential fit and the linear fit with $\mathcal{L}_{f}$ varying from 0.6
(i.e. only the first 40% of decay are considered) to 0.2 (i.e. 80% of the
decoherence curve is considered). The results are plotted in Fig. 6.

The $T_{2}$ and $p$ values obtained from either an exponential fit (circle
symbols) or a linear fit (square symbols) vary by up to $\approx$ 16% and
24%,respectively, for the dilute spin baths of 0.1 ppm and 0.2 ppm
concentrations, when $\mathcal{L}_{f}$ is varied between 0.6 and 0.2. For bath
concentrations with $\rho_{e}\geq 0.5$, the variation is smaller and in the
range of only 6% and 10% for $T_{2}$ and $p$, respectively. The large
variation of $T_{2}$ and $p$ values in the dilute spin baths is an indication
that even higher order CCE approximations than gCCE4 are required for a good
description of the coherence decay at larger times. On the other hand, gCCE3
appears to be a reasonably good approximation for the larger concentrations
$\rho_{e}\geq 0.5$ ppm, when considering the concomitant computational efforts
of the higher order cluster calculations.

We point out that the variation in the extracted $T_{2}$ and $p$ values
obtained from gCCE2 can be as large as 52% and 96%, respectively. Therefore,
one needs to be cautious using the second order cluster approximation in CCE
calculations of strongly interacting spin baths.

The almost linear dependence of $T_{2}$ on the inverse of the bath-spin
concentration (i.e. a linear dependence in the double logarithmic plot with a
slope of $\approx-1$) is well captured by our simulations. The sensitivity of
the slope of lines in Fig. 6 on the choice of $\mathcal{L}_{f}$ as well as the
choice of fitting method (exponential or linear) is very small. By decreasing
$\mathcal{L}_{f}$ from 0.6 to 0.2 the slope of the exponential (linear) fit
increase from -1.027 (-1.022) to -1.041 (-1.04). These values are in a good
agreement with the reported values of -1.01 [29] for the pure electron spin
bath and -1.06 [29], -1.02 [32] for the P1 center spin bath compared to the
experimental values -1.07 [14, 35].

Furthermore, a general trend in the relation between $T_{2}$ and $p$ is that
the covariance of these two fit parameters is always negative, meaning that
these two variables are moving in opposite directions, e.g., when changing
$\mathcal{L}_{f}$. If $T_{2}$ increases then $p$ decreases and vice versa.
This is in line with the decrements of the $p$ values in quasi two-dimensional
spin baths reported by Refs. [37, 32]. Including HF interactions or decreasing
the dimensionality of the spin baths prolongs the decay signal. As the
consequence of that, the best fit compromise results in a lowest value for
$p$.

##  V Summary

A profound knowledge of the coherence behavior of NV centers in a bath of
strongly interacting spins is of high interest for fundamental and applied
quantum technology. However, theoretical values for the coherence time $T_{2}$
and the related stretched exponential parameter $p$, as reported in
literature, show large variations. We have calculated $T_{2}$ and $p$ for NV
centers in ensembles of unpaired electron spins with concentrations of 0.1-100
ppm and analyzed the potential sources of discrepancy. We used the cluster
correlation expansion method to calculate the coherence decay. Our study
indicates that insufficiencies in the ensemble averaging as well as the choice
of parameters determining the numerical convergence of the CCE method may lead
to large variations in the resulting coherence parameters. In order to obtain
results with 5% accuracy, at least 250 configurations are required for the
ensemble average.

The results also strongly depend on the procedure used for fitting the time
evolution of the averaged coherence decay. We analyzed different fitting
approaches to obtain $T_{2}$ and $p$ values. We found that the results
obtained from fitting an exponential function to the numerical decay data may
vary up to $20\%$ compared to the results obtained from the linear fit to the
logarithmic data in a double logarithmic plot.

In order to accurately capture the dynamics of the central NV center spin in
strongly interacting electron spin baths, higher order CCE approximations than
CCE2 are required. As a compromise between accuracy and numerical workload we
chose CCE4 for dilute baths and CCE3 for higher bath concentrations. While the
use of CCE2 leads to variations of up to 50% and 100% in the extracted T2 and
$p$ values, respectively, using the higher order cluster approximations
reduces these fit uncertainties below 10%.

##  VI Acknowledgments

The authors would like to thank P. Schätzle, W. Hahn, D. Maile and J.
Stockburger for fruitful discussions. This work was funded by the Ministry of
Economic Affairs, Labour and Tourism Baden Württemberg in the context of the
Competence Center Quantum Computing Baden-Württemberg (project “SiQuRe 2”).

## References

  * [1] B. Naydenov, F. Dolde, L. T. Hall, C. Shin, H. Fedder, L. C. L. Hollenberg, F. Jelezko, J. Wrachtrup, Phys. Rev. B 83, 081201(2011).
  * [2] G. Balasubramanian, P. Neumann, D. Twitchen, M. Markham, R. Kolesov, N. Mizuochi, J. Isoya, J. Achard, J. Beck, J. Tissler, V. Jacques, P. R. Hemmer, F. Jelezko, J. Wrachtrup, Nature Materials 8, 383–387 (2009).
  * [3] N. Bar-Gill, L.M. Pham, A. Jarmola, D. Budker, R.L. Walsworth, Nature Communications 4, 1743 (2013).
  * [4] K. Saeedi, S. Simmons, J. Z. Salvail, P. Dluhy, H. Riemann, N. V. Abrosimov, P. Becker, H. J. Pohl, J. J. L. Morton, M. L. W. Thewal, Science 342, 830–833 (2013).
  * [5] G. Wolfowicz, F. J. Heremans, C. P. Anderson, S. Kanai, H. Seo, A. Gali, G. Galli, D. D. Awschalom, Nature Reviews Materials 6, 906–925 (2021).
  * [6] C. Becher, W. Gao, S. Kar, C. D. Marciniak, T. Monz, J. G. Bartholomew, P. Goldner, H. Loh, E. Marcellina, K. E. J. Goh, T. S. Koh, B. Weber, Z. Mu, J. Y. Tsai, Q. Yan, T. Huber-Loyola, S. Höfling, S. Gyger, S. Steinhauer, Val Zwiller, Mater. Quantum Technol. 3, 012501 (2023).
  * [7] F. Jelezko, T. Gaebel, I. Popa, A. Gruber, J.Wrachtrup, Phys. Rev. L 92, 076401 (2004).
  * [8] J. Cai, A. Retzker, F. Jelezko, M. B. Plenio, Nature Physics, 9 168–173 (2013).
  * [9] M. W. Doherty, F. Dolde, H. Fedder, F. Jelezko J. Wrachtrup, N. B. Manson, L. C. L. Hollenberg, Phys. Rev. B 85, 205203 (2012).
  * [10] E. Bernardi, R. Nelz, S. Sonusen, E. Neu, Phys. Rev. B 85, Crystals 7(5), 124 (2017).
  * [11] A. Gali, Nanophotonics 8, 11 (2019).
  * [12] N. Zhao, S. W. Ho, R. B. Liu, Phys. Rev. B 85, 115303 (2012).
  * [13] J. A. van Wyk, E. C. Reynhardtz, G. L. Highx, I. Kiflawik, J. Phys. D: Appl. Phys. 30 1790–1793 (1997).
  * [14] E. Bauch, S. Singh, J. Lee, C. A. Hart, J. M. Schloss, M. J. Turner, J. F. Barry, L. M. Pham, N. Bar-Gill, S. F. Yelin, R. L. Walsworth, Phys. Rev. B 102, 134210 (2020).
  * [15] L. T. Hall, J. H. Cole, L. C. L. Hollenberg, Phys. Rev. B 90 075201 (2014).
  * [16] P. L. Stanwix, L. M. Pham, J. R. Maze, D. Le Sage, T. K. Yeung, P. Cappellaro, P. R. Hemmer, A. Yacoby, M. D. Lukin, R. L. Walsworth, Phys. Rev. B 82, 201201 (R) (2010).
  * [17] J. R. Maze, J. M. Taylor, M. D. Lukin, Phys. Rev. B 78, 094303 (2008).
  * [18] F. Shi, X. Kong, P. Wang, F. Kong, N. Zhao, R. B. Liu, J. Du, Nature Physics 10, 21–25 (2014).
  * [19] W. M. Witzel, S. Das Sarma, Phys. Rev. B 74, 035322 (2006).
  * [20] W. Yang, W. L. Ma, R. B. Liu, Rep. Prog. Phys. 80, 016001 (2017).
  * [21] L. V. H. Rodgers, L. B. Hughes, M. Xie, P. C. Maurer, S. Kolkowitz, A. C. B. Jayich, N. P. de Leon, MRS Bull. 46, 623–633 (2021).
  * [22] C. Shinei, Y. Masuyama, M. Miyakawa, H. Abe, S. Ishii, S. Saiki, S. Onoda, T. Taniguchi, T. Ohshima, T. Teraji, J. Appl. Phys. 132, 214402 (2022).
  * [23] A. M. Ferrari, M. D’Amore, K. E. El-Kelany, F. S. Gentile, R. Dovesi, Journal of Physics and Chemistry of Solids 160, 110304 (2022).
  * [24] N. Nunn, S. Milikisiyants, E. O. Danilov, M. D. Torelli, L. Dei Cas, A. Zaitsev, O. Shenderova, A. I. Smirnov, J. Appl. Phys. 132, 075106 (2022).
  * [25] J. F. Barry, J. M. Schloss, E. Bauch, M. J. Turner, C. A. Hart, L. M. Pham, and R. L. Walsworth, Rev. Mod. Phys. 92, 015004 (2020).
  * [26] A. M. Edmonds, C. A. Hart, M. J. Turner, P. O. Colard, J. M. Schloss, K. S. Olsson, R. Trubko, M. L. Markha, A. Rathmill, B. Horne-Smith ,W. Lew, A. Manickam, S. Bruce, P. G. Kaup, J. C. Russo, M. J DiMario, J. T. South, J. T. Hansen, D. J Twitchen, R. L. Walsworth, Mater. Quantum Technol. 1, 025001 (2021).
  * [27] V. Stepanov, S. Takahashi, Phys. Rev. B 94, 024421 (2016).
  * [28] K. Hayashi, Y. Matsuzaki, T. Ashida, S. Onoda, H. Abe, T. Ohshima, M. Hatano, T.Taniguchi, H. Morishita, M. Fujiwara, N. Mizuochi, Journal of the Physical Society of Japan 89, 054708 (2020).
  * [29] H. Park, J. Lee, S. Han, S. Oh, H. Seo, npj Quantum Information 8, 95 (2022).
  * [30] Z. H. Wang, S. Takahashi, Phys. Rev. B 87, 115122 (2013).
  * [31] W. M. Witzel, M. S. Carroll, Lukasz Cywinski, S. D. Sarma, Phys. Rev. B 86, 035452 (2012).
  * [32] P. Schätzle, R. Ghassemizadeh, D. F. Urban, T. Wellens, P. Knittel, F. Reiter, J. Jeske, W. Hahn,  https://doi.org/10.48550/arXiv.2401.16169
  * [33] G. de Lange, Z. H. Wang, D. Riste, V. V. Dobrovitski, R. Hanson, Science 330, 60-63 (2010).
  * [34] C. Abeywardana, V. Stepanov, F. H. Cho, S. Takahashi, J. Appl. Phys. 120, 123907 (2016).
  * [35] S. Li, H. Zheng, Z. Peng, M. Kamiya, T. Niki, V. Stepanov, A. Jarmola, Y. Shimizu, S. Takahashi, A. Wickenbrock, D. Budker, Phys. Rev. B 104, 094307 (2021).
  * [36] R. Hanson, V. V. Dobrovitski, A. E. Feiguin, O. Gywat, D. D. Awschalom, Science 320, 352-355 (2008).
  * [37] E. J. Davis, B. Ye, F. Machado, S. A. Meynell, W. Wu, T. Mittiga, W. Schenken, M. Joos, B. Kobrin, Y. Lyu, Z. Wang, D. Bluvstein, S. Choi, C. Zu, A. C. Bleszynski Jayich, N. Y. Yao, Nat. Phys. 19, 836–844 (2023).
  * [38] W. Zhang, N. Konstantinidis, K. A. Al-Hassanieh,V. V. Dobrovitski, J. Phys.: Condens. Matter 19, 083202 (2007).
  * [39] W. Yang, R. B. Liu, Phys. Rev. B 78, 085315 (2008).
  * [40] W. Yang, R. B. Liu, Phys. Rev. B 79, 115320 (2009).
  * [41] G. L. Zhang, W. L. Ma, R. B. Liu,Phys. Rev. B 102, 245303 (2020).
  * [42] W. M. Witzel, K. Young, S. Das Sarma, Phys. Rev. B 90, 115431 (2014).
  * [43] M. Onizhuk, G. Galli, Adv. Theory Simul. 4, 2100254 (2021).
  * [44] M. Onizhuk, G. Galli, PRX Quantum 2, 010311 (2021).
  * [45] A. Sajid, Kristian S. Thygesen, Phys. Rev. B 106, 104108 (2022).
  * [46] M. Onizhuk, G. Galli, Phys. Rev. B 108, 075306 (2023).
  * [47] M. Ye , H. Seo, G. Galli, npj Computational Materials 5, 44 (2019).
  * [48] H. Seo, A. L. Falk, P. V. Klimov, K. C. Miao, G. Galli, D. D. Awschalom, Nature Communications 7, 12935 (2016).
  * [49] W. M. Witzel, M. S. Carroll, Lukasz Cywinski, S. D. Sarma, Phys. Rev. L 105, 187602 (2010).
  * [50] J. C. Marcks, M. Onizhuk, N. Delegan, Y. X. Wang, M. Fukami, M. Watts, A. A. Clerk, F. J. Heremans, G. Galli, D. D. Awschalom, Phys. Rev. Materials, (2024).


Computer Science • 20(4) 2019             https://doi.org/10.7494/csci.2019.20.4.3376

Track finding with Deep Neural Networks

Institute of Nuclear Physics PAN, Krakow, Poland,
e-mail<EMAIL_ADDRESS>Marcin Kucharczyk

Institute of Nuclear Physics PAN, Krakow, Poland,
e-mail<EMAIL_ADDRESS>Marcin Wolter

High-energy physics experiments require fast and efficient methods for reconstructing the tracks of charged particles. The commonly used algorithms are sequential, and the required CPU power increases rapidly with the number of tracks.
Neural networks can speed up the process due to their capability of modeling complex non-linear data dependencies and finding all tracks in parallel.

In this paper, we describe the application of a deep neural network for reconstructing straight tracks in a toy two-dimensional model.
It is planned to apply this method to the experimental data obtained by the MUonE experiment at CERN.

§ INTRODUCTION

Particle detectors used in high-energy physics (HEP) require the efficient and fast reconstruction of charged particle tracks. Effective tracking algorithms have been used for years in HEP experiments (for example, the Kalman Filter [10]), but they are sequential - the required CPU power increases rapidly with track density. Machine learning algorithms, especially Deep Neural Networks (DNN), can speed up the reconstruction process thanks to their capability to model complex non-linear data dependencies and parallelize easily.
Besides the reduction of CPU consumption, it is important how the efficiency and precision of DNN reconstruction compares with standard track reconstruction methods in the use of a neural network

Our goal is to check the potential and eventually apply the track-finding techniques based on machine learning. They are planned to be applied to the experimental testbeam data taken in 2018 by the MUonE experiment [3], which will be operating at the SPS accelerator at CERN. The experiment is supposed to measure a hadronic contribution to the anomalous muon magnetic moment in order to increase the sensitivity to the potential new physics phenomena. This might cause a discrepancy with respect to the standard model predictions, which would be a clear sign of new physics. For the MUonE experiment, the tracking is essential; deep neural networks may provide fast and efficient pattern recognition, which is the most crucial step in the track-reconstruction procedure.

§ TRACK RECONSTRUCTION

Our aim is to test whether we are able to reconstruct a particle track on a simple toy example without using sequential tracking algorithms while reconstructing the track parameters all at once using the deep neural network.
The previously trained neural network is fast and does not need much CPU power; however, the training of a DNN requires quite a bit of CPU power. For training bigger networks, graphics processor units (GPUs) are needed. Also, the amount of available memory might be a problem.
To ensure fast training on limited resources, a two-dimensional 28$\times$28 pixel toy model with straight particle tracks has been developed.
A random noise and pixel inefficiency has been added to make the model more realistic.

The input to our neural network is a two-dimensional pixel map representing the hits in the detector; no data preprocessing was done.

§.§ Reconstruction of single track

A straight track in two-dimensional space can be described by the equation representing a straight line (i.e., $y=ax+b$, with slope $a$ and intercept $b$ [see Figure <ref>]), so the neural network should return these two parameters. The parameters of the generated tracks were restricted in such a way that the track always crosses the upper and lower edges of the detector. For the described analysis, the training and test samples of 20,000 events were used for both single and multiple track events.

Single track event. Highlighted pixels are actual hits: red line represents true track; a represents track slope; and b represents intercept.

The neural network is trained to perform the regression using the entire 28$\times$28-pixel matrix as an input and returning the track parameters. Since this task is similar to the image processing, we opted to use the convolutional neural network (CNN) [15], which proved to be very powerful in the image-recognition tasks; this is routinely used in many applications.

Traditional multilayer perceptron (MLP) models have been also used for image recognition. However, due to the full connectivity between nodes, they suffer from the curse of dimensionality and do not scale well with higher-resolution images. Also, the feature extraction might be different in different regions of the image.

§.§.§ Convolutional layer

In the convolutional layer, the layer's parameters consist of a set of filters that have a relatively small receptive field. The filters are trained to activate when some specific type of feature is detected; this allows one single filter to find given features in the entire input image.
A stack of convolutional layers works as an auto-encoder, extracting more complex features in each layer.

§.§.§ Pooling layer

A pooling layer reduces the spatial size of the representation, as the exact location of a feature is less important than its approximate position relative to other features.
The use of pooling layers helps to control overfitting and reduce the amount of computation time.

§.§.§ Dropout layer

In the dropout approach at each training stage, a set of nodes randomly selected with a given probability are "dropped out" from the network. A reduced network is left and is trained during this stage. After this training stage, the removed nodes are returned to the net with their original weights.

This approach introduces some random noise to the network and decreases overtraining.

§.§.§ Network architecture

The CNN network is implemented using the KERAS [7] library with tensorflow [2] as a backend.

The network used to reconstruct single tracks consists of two convolutional layers (Conv2D) with 8 convolution windows with a size of 3$\times$3 followed by the MaxPooling layer with a pooling size of 2$\times$2. This is again followed by another two convolutional layers, each having 32 windows of a size of 3$\times$3. The single Dropout layer with a 0.25 dropout fraction suppresses the overtraining, and the final regression is performed in the Dense layer with 400 nodes and the rectified linear unit (ReLU) activation function.
This dense layer is followed by another Dropout layer with a 0.5 dropout fraction and the final Dense output layer with two nodes, which returns the reconstructed track parameters. This layer has a hyperbolic tangent activation function.

The neural network (presented in Figure <ref>) has more than 800,000 trainable parameters in total.

Convolutional deep neural network used to reconstruct track parameters for single track events.

The network is trained using the Adam optimizer [14], and the mean squared error is used as a loss function. The training is performed for 12 epochs (with a batch size of 128), which is sufficient to reach a plateau of the loss function. The evolution of the loss function during training is shown in Figure <ref>.

Evolution of loss function during training for single track events.

The model with the convolutional neural network and long/short-term memory (LSTM) layer is based on the research presented by the HEP.TrkX project [6] in Refs. [8, 9].

§.§ Reconstruction of multiple tracks

A network described above is capable of reconstructing events with a single track. For events with multiple tracks, an additional layer is needed; this transforms the output of the neural network into a sequence of track parameters. The number of reconstructed tracks might vary from event to event.

§.§.§ Long/Short-Term Memory layer

The multiple track finding problem might be treated in a similar fashion to the image-captioning, where the descriptions of the tracks (i.e., track parameters) are analogous to the text captions assigned to the various patterns seen in the image [20].
For this purpose, the long/short-term memory (LSTM) [12] layer is used.

LSTM is frequently used to process time signals; i.e., a flow of sequential inputs. In the case of track reconstruction, the CNN network sequentially reconstructs consecutive tracks, which are the sequential input for the LSTM layer.

§.§.§ Network architecture

While reconstructing multiple track events, the convolutional neural network (which is identical to the one used for single track events) extracts all relevant features from the input image. The LSTM layer stacked on top of the convolutional and dense layers emits the track parameters for each track in a sequence, thus finding the parameters of consecutive tracks. The diagram of a complete network is shown in Figure <ref>.

Convolutional deep neural network with convolutional layers followed by dense and LSTM layers. Network is trained to reconstruct track parameters for multiple track events.

The first part of the network (i.e., the convolutional layers and the first dense layer) remain unchanged. The LSTM layer with 400 units, the hyperbolic tangent activation function, and the sigmoid recurrent activation function is used. This is followed by a Dense with two units and the linear activation function. It is wrapped in the TimeDistributed wrapper, which applies this dense layer to each temporal slice of an input.

The test of multiple track reconstruction is performed for events with exactly three simulated tracks. Therefore, the last layer of the network is a Reshape layer, which forms the output into a 3$\times$2 array containing the parameters of the three tracks. During the training process, the network output is compared to the true parameters of the three simulated tracks.

The neural network described above has more than two million trainable parameters.
The network is trained in the same manner as with the single tracks. The evolution of the loss function during training is shown in Figure <ref>. The training should be stopped at this point to avoid overtraining.

Evolution of loss function during training for multiple track events.

§.§.§ Track fitting

After performing the neural network reconstruction, the hits close to the track returned by the network are used to perform the straight line fit. In order to include only relevant hits to the fit, a cut-off of five pixels around the DNN track (Euclidian distance) is chosen. Such a cut-off ensures a reasonable balance between including many noisy hits and not losing too many hits originating from a real track, since the reconstruction performed by a neural network is never perfect. The accuracy of a hit is limited to the pixel size, so all of the hits are given the same error equal to one for the fitting.

All of these hits are used to fit a track using the robust fit [16] implemented in the ROOT [5] framework. The robust fit allows us to drop the outliers and redo the fit with fewer hits. The robust fit is performed with a limit of at least 70% of the original input hits to be used in the final fit.
The fit (performed after the neural network pattern recognition) is implemented to improve the resolution of the track reconstruction done by a neural network. The fitting algorithm in the ROOT package is highly optimized, so the robust fit uses only about 5% CPU more than the standard fitting procedure does.

§ RESULTS

The track reconstruction was first performed for events with a single track. In the next step, the analysis was repeated for events with three tracks.

§.§ Single-track events

In the first step of our analyses, we tested the performance of the convolutional neural network without an LSTM layer on single-track events. The pixel efficiency was set to 70%, and the noise varied within a range of 0 to 30%.
The noise level is defined as the probability of a single pixel not belonging to the track to generate a false signal.
It should be emphasized that 30% noise is quite extreme and unrealistic for real detectors. It is only included here to test the robustness of the pattern recognition based on neural networks.

The hit positions are not smeared, but the hit accuracy is limited to the size of a pixel.
The example events with various noise levels are shown in Figure <ref>.

Single-track events with pixel efficiency of 70% and noise levels of 0 (upper left), 10% (upper right), 20% (bottom left), and 30% (bottom right). Units are pixel numbers.

Figure <ref> shows the distribution of the difference between the reconstructed and true track parameters (slope and intercept) for events with no noise and with a high 30% noise. The Gaussian is fitted to these distributions; the width of the Gaussian is taken as a measure of the resolution.
The fit range is chosen to get the width of the central part of the distribution without including the non-Gaussian tails. This is arbitrary and can influence the resolution measurement; however, it is enough to show the trend of how the resolution changes with increasing noise for the toy model described here.
The results for various noise levels are shown in Figures <ref> and <ref> and summarized in Table <ref>.

Resolution of slope (left column) and intercept (right column) for single-track events at pixel efficiency 70% and noise levels of 10% (upper row) and 30% (bottom row).

Track slope resolution obtained by neural network and robust fit (upper plot), track intercept resolution (middle plot), and reconstruction efficiency for single-track events (bottom plot) for noise levels of 0, 10, 20, and 30%.

2c| Slope     2c| Intercept [pixels]
Efficiency [%]

CNN     Fit     CNN     Fit

0% noise	    0.020      0.007     0.25     0.11     99

10% noise	    0.057      0.053     0.74     0.77     99

20% noise	    0.106      0.091     1.37     1.22     94

30% noise	    0.161      0.152     2.32     2.03     84

Slope and intercept resolution obtained by neural network (CNN) and fit (Fit) together with reconstruction efficiency for single-track events at pixel efficiency of 70% and different noise levels.

The track-reconstruction accuracy is very good for events with no noise; it degrades at higher noise levels. The robust fit applied after the neural network pattern recognition does not significantly improve the resolution.

The track is assumed to be properly reconstructed when the deep neural network finds a track, which is not further than five pixels from the true track at any point. This cut is arbitrary, but this approach enables us to show how the resolution changes with the increasing noise (which is more important than the exact measurement of the reconstruction efficiency in the toy model). The reconstruction efficiency is defined as the ratio of the number of reconstructed tracks (within the five pixel cut-off) to the total number of generated tracks.
The reconstruction efficiency is measured for the fitted tracks.
For the events without noise, the reconstruction efficiency (shown in Figure <ref>) reaches to above 99% and drops down to 84% for 30% noise.

§.§ Events with multiple tracks

In the next step, the LSTM layer was added to the network as described in Section <ref>. This layer allows us to reconstruct events with multiple tracks. Example events with three tracks are shown in Figure <ref>.

Three-track events with pixel efficiency of 70% and noise levels of 0% (upper left), 10% (upper right), 20% (bottom left), and 30% (bottom right). Units are pixel numbers.

Pattern recognition with three tracks at high noise levels is not an easy task, as the tracks are barely visible by the eye. So, it is not surprising that the reconstruction efficiency goes down to about 50% for 30% noise. At a 10% noise level, the efficiency is still reasonable (76%); for events with no noise, this reaches 92%. The resolution is significantly worse than for single-track events (see Figures <ref> and <ref> and Table <ref>).

Resolution of slope (left column) and intercept (right column) for three-track events with pixel efficiency of 70% and noise levels of 10% (upper row) and 30% (bottom row).

Track slope resolution obtained by neural network and robust fit (upper plot), track intercept resolution (middle plot), and reconstruction efficiency for events with three tracks (bottom plot) for noise levels of 0, 10, 20, and 30%.

2c| Slope     2c| Intercept [pixels]
Efficiency [%]

CNN     Fit     CNN     Fit

0% noise	    0.055      0.039     0.53     0.41     92

10% noise	    0.143      0.111     1.45     1.25     76

20% noise	    0.209      0.170     2.32     1.81     57

30% noise	    0.272      0.258     3.05     2.47     48

Slope and intercept resolution obtained by neural network (CNN) and fit (Fit) together with reconstruction efficiency for multiple-track events at pixel efficiency of 70% and different noise levels.

The example events with three tracks and various noise levels are shown in Figure <ref>.
The tracks found by the neural network (and afterward by a fit) can be further than five pixels away from the true track. In this case, they are not counted as “reconstructed” for the efficiency measurement.

For events with higher noise levels, some tracks shown in Figure <ref> remain unreconstructed or are reconstructed poorly. However, the high 30% noise level and relatively low efficiency of 70% are extreme conditions for the pattern-recognition algorithm.

Example three-track events with tracks found by neural network (blue line) and robust fit (red line) together with true tracks (green line). Hits are denoted as crosses; those marked in red were used for fitting. Events were generated with no noise (upper left plot), 10% noise (upper right plot), 20% noise (bottom left plot), and 30% noise (bottom right plot). Units are pixel numbers.

§ CPU CONSUMPTION

In the toy model presented here, the standard pattern-recognition methods (like the Kalman filter) are not applied; therefore, they cannot be compared to the track finding provided by the neural network. The comparison can be done only between the CPU's time needed to perform the DNN pattern recognition and the CPU used to perform the fit.

Table <ref> shows the wall-clock time on our PC needed for the network training, the actual pattern recognition, and the following fit for the three-track events and 10% noise.

Wall-time [sec]

Neural network training	    417

Neural network pattern recognition	     10

ROOT robust fit	    467

CPU consumption (wall-clock time on 4-core PC) by neural network pattern recognition and robust fit for three-track events with 10% noise (data samples of 20,000 events used for training and analysis)

The Tensorflow package in which the DNN is implemented is a multi-threaded application and can use all four cores of our PC. The fit implemented in ROOT can use only a single core.

The neural network requires a lot of CPU time for training; however, this time is needed only once (before starting the actual analysis). This could be further shortened if graphics processor units (GPUs) were used. Once the network is trained, the pattern recognition needs about a 50-times-lower wall-clock time than the robust fit. This corresponds to about a 12-times-lower CPU time, since the fit uses a single core only.

§ FUTURE IMPROVEMENTS OF TRACKING ALGORITHM

The most obvious extension of the presented track-finding algorithm is the reconstruction of tracks in three dimensions. This could be done by replacing the 2D convolutional layers in a neural network with 3D  layers.
The size of the network should be also increased to accomodate more complex 3D input events.

The more important improvement would be the use of the mixture density network (MDN) already proposed in 1994 [4]. In this approach, the output of the network in the case of regression is not a vector of the most probable outputs for a given input vector but a probability density of the possible outputs.

This approach makes MDNs similar to Bayesian Neural Networks (BNN) [18, 17], which also return the probability density distribution rather than the actual most probable output vectors. The advantage of MDNs is that, in contrast to BNNs, they do not require the training of multiple neural networks; this makes them much faster and easier to train and apply.

In the case of track reconstruction,
the obvious advantage of MDN is that it learns the probability distribution, therefore learning both the track position and its uncertainty at the same time in a single network. The network then returns the uncertainties of the track parameters, allowing only for the selection of well-reconstructed tracks and improving the performance of the track reconstruction.

The less obvious advantage for track finding is that the usage of the mixture density network (MDN) makes the long/short-term memory (LSTM) layer no longer necessary, therefore simplifying the network. The tracks do not need to be ordered like in the case of LSTM (which was always slightly artificial). In the MDN approach, a single Gaussian might be used to describe the probability density of each parameter of a track, and the number of Gaussians sets the limit on the number of tracks. Therefore, the network still stays flexible, being able to reconstruct any number of tracks below this limit.

§ APPLICATION IN FUTURE EXPERIMENTS

As the applied algorithms have been proven to be successful on a toy Monte Carlo model with three tracks, they are planned to be used in more-realistic conditions using testbeam data for the MUonE experiment planned to be operating at the SPS accelerator at CERN [3]. This is characterized by a relatively low occupancy; however, due to the limited acceptance and very low angles of the outgoing particles, a separation between an electron and a muon after elastic scattering is not straightforward and needs efficient and very precise track-reconstruction algorithms. Therefore, in order to achieve the physics goal (which is a measurement of the hadronic contribution to the anomalous muon magnetic moment with an uncertainty at a level of 0.05%), the deep neural network methods are going to be applied and tested with respect to the commonly used pattern-recognition and tracking procedures based on histogramming or simple line fitting.

During the CERN 2018 run, the MUonE experiment performed a feasibility test at the COMPASS experiment site, collecting a data sample size of about 1.4 billion events. Figure <ref> shows the concept of the MUonE setup; it is a modular system consisting of 16 silicon micro-strip layers, 2 carbon targets, and an electromagnetic calorimeter. The modules are separated by a relative distance of 1 m from each other and spaced by air.

Concept of MUonE detector's final layout, where input beam enters detector from left. Thick black lines denote tracking modules. Figure adopted from [3].

There are three tracks to be reconstructed in the detector; i.e., the direction and momentum of the incident muon as well as the directions of the outgoing electron and muon. Therefore, according to no CPU time limit and the expected relatively low detector occupancy, the goal of the pattern recognition is to maximize the possible track reconstruction efficiency. Standard pattern-recognition and track-finding procedures were applied, and the obtained results are going to be compared with the outcome of the parallel deep neural network method (which is a natural candidate for further improvements in the reconstruction performance). In principle, this allows for performing the pattern recognition “at once” while at the same time having great potential for increasing the precision and efficiency.

§ CONCLUSIONS

The neural network algorithm was proven to be successful when applied to the toy model.
However, it is difficult to make any judgment on whether the network described in this article can be successfully used in a real experiment for track reconstruction. This requires more studies and further development of the tracking system.

In the next step, we plan to test the track reconstruction in more-realistic conditions using testbeam data taken by the MUonE experiment. Therefore, in order to reduce the CPU timing, increasing the efficiency of the overall tracking reconstruction procedure and improving the precision the deep neural network methods are going to be applied and tested with respect to the commonly used procedures. Such tests could give an answer for whether pattern recognition based on deep neural networks can compete with standard track-finding procedures.

The toy model described in this article does not require much computing power; however, if extended to real-size events and to three dimensions, more-powerful machines or a GPU cluster would be useful for speeding up the network training.

The work was supported in part by
Polish Government Grant No. 2017/25/B/ST2/01234
and ACK CYFRONET Computing Grant No. lhcbflav08.

[1]
Keras LSTM tutorial.

[2]
Abadi M., Barham P., Chen J., Chen Z., Davis A., Dean J., Devin M., Ghemawat S., Irving G., Isard M., et al.: Tensorflow: A system for large-scale machine learning.
In: 12th $\{$USENIX$\}$ Symposium on Operating Systems Design and Implementation ($\{$OSDI$\}$ 16), pp. 265–283. 2016.

[3]
Abbiendi G., al.: Measuring the leading hadronic contribution to the muon $g-2$ via $\mu e$ scattering.
In: European Physical Journal C, vol. 77, p. 139, 2017.

[4]
Bishop C.M.: Mixture density networks.
In: , 1994.

[5]
Brun R., Rademakers F.: ROOT: An object oriented data analysis framework.
In: Nucl. Instrum. Meth., vol. A389, pp. 81–86, 1997.

[6]
Calafiura P.: HEP advanced tracking algorithms with cross-cutting applications (Project HEP.TrkX).

[7]
Chollet F., et al.: Keras: Deep learning library for theano and tensorflow, 2015.

[8]
Farrell S., Anderson D., Calafiura P., Cerati G., Gray L., Kowalkowski J., Mudigonda M., Spentzouris P., Spiropoulou M., Tsaris A., et al.: The HEP. TrkX Project: deep neural networks for HL-LHC online and offline tracking.
In: EPJ Web of Conferences, vol. 150, p. 00003. EDP Sciences, 2017.

[9]
Farrell S., Calafiura P., Mudigonda M., Anderson D., Vlimant J.R., Zheng S., Bendavid J., Spiropulu M., Cerati G., Gray L., et al.: Novel deep learning methods for track reconstruction.
In: arXiv preprint arXiv:1810.06111, 2018.

[10]
Frühwirth R.: Application of Kalman filtering to track and vertex fitting.
In: Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 262(2-3), pp. 444–450, 1987.

[11]
Graves A., Wayne G., Danihelka I.: Neural Turing Machines.
In: arXiv preprint arXiv:1410.5401, 2014.

[12]
Hochreiter S., Schmidhuber J.: Long short-term memory.
In: Neural computation, vol. 9(8), pp. 1735–1780, 1997.

[13]
James F., Roos M.: MINUIT: a system for function minimization and analysis of the parameter errors and corrections.
In: Comput. Phys. Commun., vol. 10(CERN-DD-75-20), pp. 343–367, 1975.

[14]
Kingma D.P., Ba J.: Adam: A method for stochastic optimization.
In: arXiv preprint arXiv:1412.6980, 2014.

[15]
LeCun Y., Bengio Y., Hinton G.: Deep learning.
In: nature, vol. 521(7553), p. 436, 2015.

[16]
Motulsky H.J., Brown R.E.: Detecting outliers when fitting data with nonlinear regression–a new method based on robust nonlinear regression and the false discovery rate.
In: BMC bioinformatics, vol. 7(1), p. 123, 2006.

[17]
Neal R.M.: Bayesian learning for neural networks, vol. 118.
Springer Science & Business Media, 2012.

[18]
Pearl J.: Markov and Bayesian Networks, chap. 3 Probabilistic Reasoning in Intelligent Systems, 1988.

[19]
Srivastava N., Hinton G., Krizhevsky A., Sutskever I., Salakhutdinov R.: Dropout: a simple way to prevent neural networks from overfitting.
In: The journal of machine learning research, vol. 15(1), pp. 1929–1958, 2014.

[20]
Vinyals O., Toshev A., Bengio S., Erhan D.: Show and tell: A neural image caption generator.
In: Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 3156–3164. 2015.