# Rocks, Water and Noble Liquids: Unfolding the Flavor Contents of Supernova
Neutrinos

Sebastian Baum  Francesco Capozzi  Shunsaku Horiuchi

###### Abstract

Measuring core-collapse supernova neutrinos, both from individual supernovae
within the Milky Way and from past core collapses throughout the Universe (the
diffuse supernova neutrino background, or DSNB), is one of the main goals of
current and next generation neutrino experiments. Detecting the heavy-lepton
flavor (muon and tau types, collectively $\nu_{x}$) component of the flux is
particularly challenging due to small statistics and large backgrounds. While
the next galactic neutrino burst will be observed in a plethora of neutrino
channels, allowing to measure a small number of $\nu_{x}$ events, only upper
limits are anticipated for the diffuse $\nu_{x}$ flux even after decades of
data taking with conventional detectors. However, paleo-detectors could
measure the time-integrated flux of neutrinos from galactic core-collapse
supernovae via flavor-blind neutral current interactions. In this work, we
show how combining a measurement of the average galactic core-collapse
supernova flux with paleo detectors and measurements of the DSNB electron-type
neutrino fluxes with the next-generation water Cherenkov detector Hyper-
Kamiokande and the liquid noble gas detector DUNE will allow to determine the
mean supernova $\nu_{x}$ flux parameters with precision of order ten percent.

##  1 Introduction

Supernovae play an important role for the workings of our Universe as well as
for our understanding thereof. The rapid injection of energy from supernovae
is crucial for the properties of galaxies, and the extreme environment of a
supernova makes them the birthplaces of much of the heavy elements in our
Universe as well as prolific sources of energetic cosmic rays. Neutrinos play
a crucial role in the explosion mechanisms of supernovae, not least because
the matter density during the initial stages of a core-collapse supernova is
so high that neutrinos are the only (known) particles which can escape this
environment. While core-collapse supernovae produce enormous fluxes of all
flavors of neutrinos, to date, the only supernova neutrinos measured by
terrestrial experiments are the historical $\bar{\nu}_{e}$ events detected
from SN1987A [1, 2, 3, 4]. Detecting all flavors of supernova neutrinos is an
important goal, for example to infer the total energetics of the collapse, to
test complex neutrino oscillation phenomena, and to search for physics beyond
the Standard Model [5].

In the upcoming decades, the next generation of neutrino detectors, in
particular Hyper-Kamiokande [6] and DUNE [7], will provide unprecedented
sensitivity to supernova neutrinos. Utilizing charged-current interactions,
Hyper-Kamiokande is predominantly sensitive to $\bar{\nu}_{e}$ neutrinos while
DUNE is mostly sensitive to the $\nu_{e}$ fluxes. In order to detect heavy
lepton flavor neutrinos
$\nu_{x}=\\{\nu_{\mu},\bar{\nu}_{\mu},\nu_{\tau},\bar{\nu}_{\tau}\\}$, one
must rely on processes with significantly smaller event rates. One strategy to
detect $\nu_{x}$ is to use scattering on protons in scintillator detectors [8,
9, 10, 11, 12]. Another strategy is to exploit flavor-blind channels, for
example, elastic scattering events in large direct dark matter detection
experiments, and subtract the (anti-)electron-type neutrino contribution by
merging data sets.

If next-generation neutrino detectors are operational at the time of the next
galactic core-collapse event, supernova neutrinos of all flavors will be
measured with relatively large statistics - Hyper-Kamiokande, DUNE and JUNO,
would detect $\mathcal{O}(10^{5})$ $\bar{\nu}_{e}$, $\mathcal{O}(10^{3})$
$\nu_{e}$ and $\mathcal{O}(10^{2})$ $\nu_{x}$ events, respectively. In
addition, the current generation of tonne-scale dark matter detectors such as
XENONnT, LZ, or PandaX-4T, would make flavor-blind detections of neutrinos
from such a galactic core-collapse event [13, 14]. While neutrinos from the
next galactic core-collapse supernova will teach us invaluable lessons about
neutrinos and supernova physics, this approach also comes with drawbacks: The
galactic core-collapse supernova rate is estimated to be $2-3$ per century
[15, 16, 17, 18, 19, 20] and the only remedy is patience. Furthermore, one
would measure the neutrinos from an individual core-collapse supernova which
may well be atypical. Both of these issues are mitigated by attempting to
detect neutrinos from the Diffuse Supernova Neutrino Background (DSNB), the
neutrino flux from distant core-collapse supernovae throughout the Universe
[21, 22]. However, detecting the DSNB is challenging due to its relatively
small neutrino flux and because the DSNB spectrum is significantly redshifted
- the DSNB flux is dominated by supernovae at redshifts $z\sim 1-2$. The next
generation of neutrino detectors, in particular Hyper-Kamiokande and DUNE, are
expected to detect $\bar{\nu}_{e}$ and $\nu_{e}$ neutrinos from the DSNB and
measure the corresponding neutrino flux parameters. However, no (conventional)
neutral-current detector with the required combination of exposure and energy
threshold is available to measure the $\nu_{x}$ fluxes. Even future direct
dark matter detectors such as DARWIN [23] are expected to set only upper
limits on the DSNB $\nu_{x}$ flux roughly an order of magnitude above the
expected flux [24].

In this work, we show the prospects of paleo detectors to measure the
$\nu_{x}$ supernova neutrino flux. Rather than operating a real-time
laboratory experiment to search for neutrinos, paleo detectors [25, 26, 27]
would use nuclear damage tracks recorded in natural minerals of geological
ages (100 Myr$-$1 Gyr) to search for neutrino-induced nuclear recoils [28].
Many minerals are excellent solid state track detectors [29, 30, 31, 32] and
retain the damage tracks caused by nuclear recoils for timescales which can
exceed the age of the Earth by many orders of magnitude. Modern microscopy
techniques such as Small Angle X-Ray scattering [33, 34, 35] potentially allow
for the readout of these damage tracks in macroscopic samples. Paleo detectors
have been considered to study atmospheric [36] and solar neutrinos [37, 38] in
addition to supernova neutrinos [28].

There are two crucial differences between detecting supernova neutrinos in
paleo detectors and in Hyper-Kaomiokande/DUNE: First, the geological time
scales paleo detectors record neutrino-interactions over are many orders of
magnitude longer than the (inverse) supernova rate of the Milky Way ($2-3$ per
century [15, 16, 17, 18, 19, 20]). Physically then, paleo detectors are
sensitive to the mean supernova neutrino properties from many past core
collapses, similar to the DSNB modulo different distances. When averaged over
timescales much longer than a century, the supernova neutrino flux on Earth is
no longer dominated by the DSNB, but by neutrinos from supernovae in our
galaxy, and this time-averaged flux is a factor of order 100 larger than the
DSNB flux [28]. Second, paleo detectors would detect supernova neutrinos via
Coherent Elastic neutrino-Nucleus Scattering (CE$\nu$NS), a process mediated
by neutral currents and hence sensitive to the total neutrino flux rather than
the $\bar{\nu}_{e}$/$\nu_{e}$ fluxes Hyper-Kamiokande/DUNE would measure. In
this regard, paleo detectors utilize the same interaction physics as direct
dark matter experiments. Combined with the significantly larger time-
integrated galactic supernova neutrino flux, paleo detectors offer a unique
opportunity to measure all flavors of neutrinos emitted from a population of
past core-collapse events.

In order to extract the $\nu_{x}$ flux parameters, a measurement of the total
galactic supernova neutrino flux in paleo detectors must be combined with
flavor-sensitive measurements of supernova neutrinos. To this end, we consider
a combination of the paleo-detector measurement with measurements of the DSNB
$\bar{\nu}_{e}$/$\nu_{e}$ parameters with Hyper-Kamiokande/DUNE. The main
background for supernova neutrinos in paleo detectors stems from radiogenic
neutrons originating from the 238U decay chain, and we project the precision
to which the $\nu_{x}$ flux parameters can be extracted for different
assumptions on the uranium concentration in paleo-detector samples. We find
that for a 238U concentration of $10^{-11}\,$g/g ($10^{-12}\,$g/g), the
average energy as well as the total energy of the $\nu_{x}$ flux can be
measured with few tens of percent ($\lesssim 10\%$) precision.

The remainder of this work is structured as follows: In section 2, we discuss
the calculation of the time-averaged flux of neutrinos from galactic
supernovae at Earth, the calculation of the DSNB and the calculation of the
neutrino-induced signal in paleo detectors. A discussion of the most relevant
background sources in paleo detectors can be found in section 2.1. In section
3 we explain the statistical analysis we use to extract the neutrino
parameters from mock data in Hyper-Kamiokande, DUNE, and paleo detectors. We
show our results for the projected precision of the reconstruction of these
parameters in section 4. In section 5, we summarize and discuss our findings.

##  2 Paleo detectors as (supernova) neutrino detectors

In this section, we discuss the signal from supernova neutrinos in paleo
detectors as well as the most relevant backgrounds. For a more detailed
discussion of paleo detectors as (CE$\nu$NS) neutrino detectors, see Ref.
[28], and for discussions of the relevant background sources in paleo
detectors see also Refs. [26, 27].

We model the spectrum of the neutrino species $\nu_{i}$ from a supernova as a
pinched Fermi-Dirac distribution [39]

$\left(\frac{dn}{dE}\right)_{\nu_{i}}=E_{\nu}^{\rm
tot}\frac{\left(1+\alpha\right)^{1+\alpha}}{\Gamma(1+\alpha)}\frac{E^{\alpha}}{\langle
E_{\nu}\rangle^{2+\alpha}}e^{\left[-\left(1+\alpha\right)\frac{E}{\langle
E_{\nu}\rangle}\right]}\;,$ (2.1)

where $E_{\nu}^{\rm tot}$ is the total energy radiated in $\nu_{i}$, $\langle
E_{\nu}\rangle$ is the average energy of the neutrinos, and $\alpha$ is a
(dimensionless) spectral shape parameter. For flavor-blind interactions such
as CE$\nu$NS, the relevant neutrino flux is the sum over all flavors,

$\frac{dn}{dE_{\nu}}=\left(\frac{dn}{dE}\right)_{\nu_{e}}+\left(\frac{dn}{dE}\right)_{\bar{\nu}_{e}}+4\left(\frac{dn}{dE}\right)_{\nu_{x}}\;.$
(2.2)

The time-averaged neutrino spectrum from core-collapse supernovae in our
Galaxy is then obtained by integrating over the probability density $f(R_{E})$
for a supernova to occur at a distance $R_{E}$ from Earth and by multiplying
with the galactic core-collapse supernova rate, $\dot{N}_{\rm CC}^{\rm gal}$,

$\left(\frac{d\phi}{dE_{\nu}}\right)^{\rm gal}=\dot{N}_{\rm CC}^{\rm
gal}\frac{dn}{dE_{\nu}}\int_{0}^{\infty}dR_{E}\,\frac{f(R_{E})}{4\pi
R_{E}^{2}}\;.$ (2.3)

For $f(R_{E})$ we assume that supernovae occur predominantly in the stellar
disk and model the distribution of core-collapse supernovae as a double-
exponential in galactocentric coordinates,

$\rho\propto e^{-R/R_{d}}e^{-\left|z\right|/H}\;,$ (2.4)

where $R$ is the galactocentric radius, $z$ the height above the galactic
plane, and we use $R_{d}=2.9\,$kpc and $H=95\,$pc [20].

Similarly, the DSNB neutrino flux is obtained by integrating over the
volumetric cosmic supernova rate, $\dot{n}_{\rm CC}^{\rm cosmo}$, [21, 22]

$\left(\frac{d\phi}{dE_{\nu}}\right)^{\rm
DSNB}=\int_{0}^{\infty}\frac{dz}{H_{0}\sqrt{\Omega_{\Lambda}+\Omega_{m}\left(1+z\right)^{3}}}\dot{n}_{\rm
CC}^{\rm cosmo}(z)\frac{dn}{dE^{\prime}_{\nu}}\;,$ (2.5)

where $H_{0}\simeq 67\,$km/s/Mpc is the Hubble constant,
$\Omega_{\Lambda}\simeq 0.68$ and $\Omega_{m}\simeq 0.32$ are the cosmological
density parameters for dark energy and matter, respectively, and $z$ is the
cosmological redshift. Note that the supernova neutrino spectrum,
$dn/dE_{\nu}$, must be evaluated at the redshifted neutrino energy
$E_{\nu}^{\prime}=(1+z)E_{\nu}$. The cosmic supernova rate is related to the
(volumetric) star formation rate $\dot{\rho}_{*}(z)$ via the conversion factor

$\dot{n}_{\rm CC}^{\rm cosmo}(z)=\dot{\rho}_{*}(z)\frac{\int_{M_{\rm
min}=8}^{100}\psi(M)\,dM}{\int_{0.1}^{100}M\psi(M)\,dM},$ (2.6)

which is equivalent to assuming that all massive stars with mass above $M_{\rm
min}=8M_{\odot}$ undergo core collapse. This is a reasonable assumption for
our estimates (see, e.g., discussions in [40]) – although a fraction of
massive stars are expected to collapse to black holes, with systematically
higher temperature neutrino emissions [41], we neglect their contribution in
this work. Firstly, their occurrence rate remains highly uncertain, and
secondly, their contributions would enhance the detectability of the DSNB so
we remain conservative by omitting their occurrence. The star-formation rate
is observationally measured [42, 43], especially well in the low redshift
range of importance for estimating the DSNB, and we adopt the functional fit
described in, e.g., Refs. [41, 44]. Note that alternatively one could also
take $\dot{n}_{\rm CC}^{\rm cosmo}(z)$ directly from supernova surveys, see
for example Ref. [45]. However, the limited depth of current supernova surveys
means that in the most relevant redshift range, $z\lesssim 2$, the errors of
such direct measurements of $\dot{n}_{\rm CC}^{\rm cosmo}(z)$ are larger than
those of inferring $\dot{n}_{\rm CC}^{\rm cosmo}(z)$ from measurements of the
cosmic star formation rate via eq. (2.6).

Supernova neutrinos (predominantly111Additional signals arise via quasi-
elastic charged-current interactions and, for more energetic neutrinos, via
(deep) inelastic neutrino-nucleus interactions [36]. However, the contribution
from quasi-elastic interactions is suppressed compared to the CE$\nu$NS one by
the lack of coherent enhancement, and the contribution of inelastic
interactions of more enegetic neutrinos is suppressed by the quickly falling
neutrino flux.) induce signals in paleo detectors via CE$\nu$NS interactions.
The differential recoil spectrum (per unit target mass) for a target nucleus
$T$ to recoil with energy $E_{R}$ is given by [46, 47]

$\left(\frac{dR}{dE_{R}}\right)_{T}=\frac{1}{m_{T}}\int_{E_{\nu}^{\rm
min}}dE_{\nu}\,\frac{d\sigma}{dE_{R}}\frac{d\phi}{dE_{\nu}}\;,$ (2.7)

where $m_{T}$ is the mass of $T$ and $E_{\nu}^{\rm min}=\sqrt{m_{T}E_{R}/2}$
is the minimum neutrino required to induce a nuclear recoil with $E_{R}$. The
differential CE$\nu$NS cross section is

$\frac{d\sigma}{dE_{R}}(E_{R},E_{\nu})=\frac{G_{F}^{2}}{4\pi}Q_{W}^{2}m_{T}\left(1-\frac{m_{T}E_{R}}{2E_{\nu}^{2}}\right)F^{2}(E_{R})\;,$
(2.8)

with the Fermi coupling constant $G_{F}$, the nuclear form factor
$F(E_{R})$222We use the Helm form factor [48] as in Ref. [28], and

$Q_{W}\equiv\left(A_{T}-Z_{T}\right)-\left(1-4\sin^{2}\theta_{W}\right)Z_{T}\;,$
(2.9)

where $\theta_{W}$ is the weak mixing angle, and $A_{T}$ ($Z_{T}$) the number
of nucleons (protons) in $T$.

The observable in a paleo detector is the track length spectrum, which is
obtained from the recoil energy spectrum by summing over all target nuclei and
weighting with the stopping power of $T$ in the material, $dE_{R}/dx_{T}$,

$\frac{dR}{dx}=\sum_{T}\xi_{T}\frac{dE_{R}}{dx_{T}}\left(\frac{dR}{dE_{R}}\right)_{T}\;,$
(2.10)

where $\xi_{T}$ is the mass fraction of $T$ in the target material. Similarly,
the length $x$ of a track from a nucleus $T$ recoiling with energy $E_{R}$ is
obtained by integrating the stopping power,

$x(E_{R})=\int_{0}^{E_{R}}dE\,\left|\frac{dE}{dx_{T}}\right|^{-1}\;.$ (2.11)

We use SRIM [49, 50] to calculate stopping powers in the target materials used
in this work, epsomite [Mg(SO${}_{4})\\!\cdot\\!7($H2O)], gypsum
[Ca(SO${}_{4})\\!\cdot\\!2$(H2O)], and sinjarite
[CaCl${}_{2}\\!\cdot\\!2($H2O)].

###  2.1 Background Sources

Let us briefly discuss the most relevant background sources for supernova
neutrino searches in paleo detectors. These backgrounds have been discussed in
detail in Refs. [26, 28] see also Ref. [27]. Note that all relevant
backgrounds stem from nuclear recoils. Natural defects in minerals are either
single-site or span across the entire mono-crystalline volume, and thus do not
resemble the damage tracks induced by CE$\nu$NS of supernova neutrinos.

  * •

Cosmogenics: In order to suppress cosmic-ray induced backgrounds, minerals
used as paleo detectors must have been shielded by a sufficiently large
overburden for the entire time they have been recording nuclear damage tracks.
The modest amounts of target materials required for paleo detectors, at most a
few kg, can be sourced from much greater depths than those of existing
underground laboratories. For example, boreholes drilled for geological R&D or
oil exploration are promising sources of samples. For an overburden of
$\gtrsim 5\,$km rock, cosmogenic backgrounds, in particular the cosmogenic
muon induced neutron flux, are suppressed to a negligible level.

  * •

Radiogenics: Any mineral used as a paleo detector will be contaminated by
trace amounts of radioactive material. In order to mitigate the associated
radiogenic backgrounds, it is crucial to use as radiopure minerals as possible
as paleo detectors. The most important radioactive isotope in natural minerals
is 238U. Minerals formed in the Earth’s crust have typical 238U concentrations
of $C^{238}\sim 10^{-6}\,$g/g, which would lead to prohibitively large
radiogenic backgrounds. The 238U concentration in so-called Ultra Basic Rocks
(UBRs), formed from the material of the Earth’s mantle, and in Marine
Evaporites (MEs), salts formed from sea water, are much lower, making them
attractive as paleo detectors, see Refs. [51, 52, 53, 54, 55, 56, 57] as well
as the discussion in Refs. [26, 28]. As in previous work on paleo detectors
[25, 26, 58, 28, 36, 27], we will assume benchmark 238U concentrations of
$C^{238}=10^{-11}\,$g/g for the ME examples epsomite
[Mg(SO${}_{4})\\!\cdot\\!7($H2O)], gypsum [Ca(SO${}_{4})\\!\cdot\\!2$(H2O)],
and sinjarite [CaCl${}_{2}\\!\cdot\\!2($H2O)] we use as target materials in
this work. We will also show prospects for measuring supernova neutrinos for a
more optimistic assumption of $C^{238}=10^{-12}\,$g/g. The most relevant
radiogenic background in paleo detectors are nuclear recoils induced by the
elastic scattering of radiogenic neutrons off the nuclei constituting the
mineral. These neutrons are produced by spontaneous fission of heavy
radiogenic isotopes as well as by $(\alpha,n)$-reactions of the
$\alpha$-particles from the 238U decay chain. We use SOURCES-4A [59] to model
these neutron spectra, and calculate the induced nuclear recoil spectrum using
TENDL-2017 [60, 61, 62, 63] neutron-nucleus cross sections as tabulated in the
JANIS4.0 database [64].333We take only elastic neutron-nucleus scattering into
account; this yields a conservative estimate of the background because
neutrons typically lose a larger fraction of their energies through inelastic
processes than in elastic scattering. Note also that our Monte Carlo
simulation of the nuclear recoils induced by radiogenic neutrons has been
validated with a calculation of the nuclear recoils induced by the same
neutron spectra with FLUKA [65, 66, 67] in Ref. [36] for the particular case
of a halite target.

  * •

Astrophysical neutrinos: Just like supernova neutrinos, neutrinos from other
sources give rise to nuclear recoils and, in turn, nuclear damage tracks in
paleo detectors. The most relevant backgrounds to supernova neutrino induced
recoils are solar and atmospheric neutrinos. We take the corresponding
neutrino fluxes from Ref. [68] and model the corresponding nuclear recoil
spectrum as for the supernova neutrinos described above.

Figure 1: Number of events per bins of track length in a Epsomite
[Mg(SO${}_{4})\\!\cdot\\!7($H2O)] paleo detector with 100 gram Gyr exposure.
Throughout this work, we use 70 log-spaced bins in the range $7.5\,{\rm
nm}\leq x\leq 10^{3}\,$nm. The red lines show the events induced by CE$\nu$NS
of neutrinos from core-collapse supernovae in our galaxy. The blue line
indicates the spectrum induced by scattering of solar, atmospheric, and DSNB
neutrinos, and the black line shows the spectrum induced by radiogenic
neutrons. As we show in this work, by utilizing the entire spectrum the
galactic supernova signal and its properties can be measured although it
produces fewer tracks than the backgrounds for all track lengths.

In figure 1 we show the signal induced by the coherent scattering of neutrinos
from supernovae in our Galaxy in epsomite [Mg(SO${}_{4})\\!\cdot\\!7($H2O)]
together with the backgrounds induced by other astrophysical neutrinos and
radiogenic neutrons, assuming $C^{238}=10^{-11}\,$g/g. Note that we have
included contributions from solar and atmospheric neutrinos as well as from
the DSNB in the spectrum labeled “$\nu$ bkg”. As we will see, the properties
of the galactic supernova neutrino flux can be inferred in a paleo detector
despite the fact that the overall normalization of the associated track length
spectrum is smaller than that of the backgrounds. This is possible thanks to
the enormous exposure achievable in paleo detectors: Reading out 100 g of a
mineral which has been recording damage tracks for 1 Gyr corresponds to an
exposure of $\varepsilon=100\,{\rm g\,Gyr}=100\,$kt yr. For such an exposure,
there would be about $2.8\times 10^{3}$ $\nu_{x}$ tracks with length between
100 and 500 nm. In the same range, there would be $1.8\times 10^{5}$ tracks
from background events for our fiducial assumption of $C^{238}=10^{-11}\,$g/g.
Therefore, considering only statistical uncertainty, the $\nu_{x}$ signal can
be observed with a statistical significance of $2.8\times
10^{3}/\sqrt{1.8\times 10^{5}}=6.6\sigma$. For the same exposure in gypsum
[Ca(SO${}_{4})\\!\cdot\\!2$(H2O)] and sinjarite
[CaCl${}_{2}\\!\cdot\\!2($H2O)], there would be a significance of $3\sigma$
and $2\sigma$, respectively. Note that even for 238U concentrations as large
as $C^{238}=5\times 10^{-11}\,$g/g, the (statistical) significance for seeing
$\nu_{x}$ induced tracks in an epsomite paleo detector would still be
$3\sigma$.

##  3 Analysis

In order to calculate the track length distributions in paleo detectors for
both supernova neutrino and background events we use the public code presented
in Ref. [27]. In our analysis we employ three rock samples: Epsomite
[Mg(SO${}_{4})\\!\cdot\\!7($H2O)], gypsum [Ca(SO${}_{4})\\!\cdot\\!2$(H2O)]
and sinjarite [CaCl${}_{2}\\!\cdot\\!2($H2O)]; these minerals provide
particularly good detection statistics for supernova neutrinos. We assume that
100 g of each material can be read out with track length resolution of 15
nm.444In order to include the effects of finite track length resolution, we
compute the number of tracks in the $i$-th bin with reconstructed length
$x\in\left[x_{i}^{\rm min},x_{i}^{\rm max}\right]$ in a sample of mass $M$
which has been recording tracks for a time $t_{\rm age}$ as $N_{i}=M\times
t_{\rm age}\int dx^{\prime}\,W(x^{\prime};x_{i}^{\rm min},x_{i}^{\rm
max})\left(dR/dx^{\prime}\right)$, where $dR/dx^{\prime}$ is the true track
length spectrum and $W$ is a window function, see, for example, Ref. [27].
Such a sample size and readout resolution may, for example, be achievable with
small angle X-ray scattering microscopy [33, 34, 35], see the discussion in
Ref. [26]. We consider both a default 238U concentration of
$C^{238}=10^{-11}\,$g/g and a more optimistic one of $C^{238}=10^{-12}\,$g/g.
These choices are summarized in table 1. Concerning the neutrino energy
spectrum, the parameters used for generating our set of mock data are given in
table 2. These have been extracted from figure 1 in [69], where the first 3D
supernova simulation up to 7 seconds post bounce has been presented, providing
a more reliable estimate of the total energies of each neutrino species. In
particular, the values for the average energies reported in table 2 refer to
the time-averaged quantities reported in figure 1 of [69], whereas the total
energies are obtained by integrating over time.

Table 1: Summary of paleo detector properties used in our simulation. Track Resolution |  $15\,$nm
---|---
Track Length Min |  $7.5\,$nm
Track Length Max |  $1000\,$nm
Number of Log Bins | 70
Rock Sample Mass |  $100\,$g
Rock Sample Age |  $1\,$Gyr
Default Uranium Concentration ($C^{238}$) |  $10^{-11}\,$g/g
Optimistic Uranium Concentration ($C^{238}$) |  $10^{-12}\,$g/g
Table 2: Summary of neutrino properties used to generate our set of mock data. Note that $E^{\rm tot}_{\nu_{x}}$ refers to the total energy of one of the four $\nu_{x}$ neutrino species ($\nu_{\mu}$, $\bar{\nu}_{\mu}$, $\nu_{\tau}$, $\bar{\nu}_{\tau}$), see eq. (2.2), and $\alpha_{\nu}$ is the pinching-parameter of the energy distribution of neutrinos [39], see eq. (2.1). $\langle E_{\nu_{e}}\rangle$ |  $11.2$ MeV
---|---
$\langle E_{\bar{\nu}_{e}}\rangle$ |  $13.5$ MeV
$\langle E_{\nu_{x}}\rangle$ |  $13.4$ MeV
$E^{\rm tot}_{\nu_{e}}$ |  $6.8\times 10^{52}$ erg
$E^{\rm tot}_{\bar{\nu}_{e}}$ |  $6.6\times 10^{52}$ erg
$E^{\rm tot}_{\nu_{x}}$ |  $6.2\times 10^{52}$ erg
$\alpha_{\nu_{e}}=\alpha_{\nu_{x}}=\alpha_{\bar{\nu}_{e}}$ | 3
Galactic Supernova Rate |  $0.023$ per year

Since the contribution of galactic supernova neutrinos to the track length
distribution is sub-dominant with respect to the radiogenic and solar neutrino
background, paleo detectors alone cannot (completely) break the degeneracy
between the three supernova neutrino components
($\nu_{e},\bar{\nu}_{e},\nu_{x}$). To alleviate this problem we include an
independent measurement of the DSNB with the future neutrino detectors DUNE
[70, 7, 71, 72] and Hyper-Kamiokande [6], which will be sensitive to $\nu_{e}$
and $\bar{\nu}_{e}$, respectively. Let us note that observations of the
neutrinos from a single future galactic supernova in neutrino and dark matter
detectors can, in general, provide much higher precision measurements of the
supernova neutrino spectrum than the combined measurement of the DSNB in DUNE
and Hyper-Kamiokande and galactic supernova neutrinos in paleo detectors we
consider here. However, our approach has two important advantages. First, we
do not need to wait for the next galactic core-collapse supernova. Second, our
approach would measure neutrinos from a large population of supernovae rather
than from one specific (possibly atypical) supernova.

In order to project the precision with which DUNE and Hyper-Kamiokande can
constrain the parameters controlling the supernova neutrino spectrum, we
compute the DSNB flux as in eq. (2.5), except that we use the spectrum from
the appropriate single neutrino flavor $\nu_{i}$ rather than the total
neutrino flux, see eqs. (2.1) and (2.2). For detection at DUNE, we consider a
40 kton liquid argon detector and detection of electron neutrinos through the
charged-current interaction $\nu_{e}+{\rm Ar}\to e^{-}+{\rm K}^{+}$. We adopt
the cross section based on the random phase approximation scheme of Ref. [73].
Since DUNE is under construction, we assume a detection efficiency of 86%
based on the DUNE design report [7, 74]. For detection at Hyper-Kamiokande, we
consider a gadolinium-enhanced water detector with fiducial volume of 187 kton
and detection of anti-electron neutrinos through inverse-beta interaction on
water $\bar{\nu}_{2}+p\to n+e^{+}$. The cross section is well-known [75, 76],
and we assume a detection efficiency of 67% [74]. We consider a 20-year run
time for both Hyper-Kamiokande and DUNE. We note that although other
detectors, e.g., JUNO [77], have sensitivity to DSNB neutrinos [78], DUNE and
Hyper-Kamiokande will dominate the sensitivity to $\nu_{e}$ and
$\bar{\nu}_{e}$, respectively.

In our analysis we perform a scan over a 12-dimensional parameter space. Six
parameters correspond to the average and total energies of
$\nu_{e},\bar{\nu}_{e}$ and $\nu_{x}$ (the pinching parameters
$\alpha_{\nu_{i}}$ are fixed to the values in table 2). The floating
parameters related to paleo detectors are the normalization of the radiogenic
and solar neutrino background. Those related to DSNB measurements are the
normalization of the backgrounds in DUNE and Hyper-Kamiokande, for which we
adopt the spectra modeled in Refs. [78, 74]. The last two parameters are the
galactic supernova rate and the overall normalization of the DSNB signal,
which parametrize the uncertainty on the star formation rate. We introduce
Gaussian priors with (relative) width of 10% on the galactic supernova rate
and the normalization of the DSNB, and priors with 20% width on the background
normalization in DUNE and Hyper-Kamiokande. To put our choices of priors on
the galactic supernova rate and the DSNB in contex, we note that the
uncertainty on the cosmic star formation rate is already at the level of $\sim
20\%$ at the low redshifts most relevant to the DSNB and local core collapses
[42, 21, 43]. We expect considerable improvements on this uncertainty in the
next decade, not least from the upcoming surveys with the James Webb Space
Telescope [79] and the Vera C. Rubin Observatory [80]. While the time
evolution of the galactic star formation rate on gigayear timescales remains
uncertain with tens-of-percent fluctuations (see, e.g., Refs. [81, 82, 83]),
we emphasize that the relevant quantity for our analysis is the local
supernova rate integrated over the past $\sim 1\,$Gy rather than the
instantaneous rate. Note that we are not including any priors on the
normalization of the solar neutrino and radiogenic backgrounds in paleo
detectors nor on the total and mean energies of the supernova neutrino fluxes.
We scan over the parameter space using a Markov Chain Monte Carlo through the
public software PyMultinest [84, 85]. We define the following likelihood:

$\begin{split}\log(\mathcal{L})=&-\sum_{i=1}^{N_{\rm bins}^{\rm
paleo}}\left[N_{i}^{\rm th}-N_{i}^{\rm data}+N_{i}^{\rm
data}\,\left(\log(N_{i}^{\rm data})-\log(N_{i}^{\rm
th}\right)\right]-\frac{1}{2}\sum_{p}\left(\frac{f_{p}}{\sigma_{p}}\right)^{2}\\\
&-\sum_{\rm exp}\sum_{i=1}^{N_{\rm bins}^{\rm exp}}\left[M_{{\rm exp},i}^{\rm
th}-M_{{\rm exp},i}^{\rm data}+M_{{\rm exp},i}^{\rm data}\,\left(\log(M_{{\rm
exp},i}^{\rm data})-\log(M_{{\rm exp},i}^{\rm th}\right)\right]\,,\end{split}$
(3.1)

where $N_{i}^{\rm data}$ is the number of simulated data events in the $i$-th
track length bin for paleo detectors and $N_{i}^{\rm th}$ the number of
expected events, defined as $N_{i}^{\rm th}=N_{i}^{\rm
0,th}(1+\sum_{p}f_{p})$, where the $f_{p}$ are the free normalization
parameters discussed above and $N_{i}^{\rm 0,th}$ is the expected number of
events calculated with the default values of parameters. $M_{{\rm exp},i}^{\rm
data,th}$ is defined analogously to $N_{i}^{\rm th}$, but refers to DSNB
events for exp$=$(Hyper-Kamiokande,DUNE). Concerning the binning, we use 70
log-spaced bins between $7.5\,{\rm nm}$ and $10^{3}\,$nm.

For the DSNB we take into account also possible neutrino oscillations inside
the supernova.555Since paleo detectors measure supernova neutrinos
(predominantly) via CE$\nu$NS, flavor oscillations have no effect on the
signal in paleo detectors. In particular, focusing only on the MSW effect and
neglecting collective and turbulence effects, we consider three cases: no
oscillations, normal and inverted mass ordering. For normal mass ordering the
survival probability of $\nu_{e}$ is $P(\nu_{e}\to\nu_{e})=0$, whereas the one
for $\bar{\nu}_{e}$ is $P(\bar{\nu}_{e}\to\bar{\nu}_{e})=0.7$. For inverted
mass ordering the survival probability of $\nu_{e}$ is
$P(\nu_{e}\to\nu_{e})=0.3$, whereas the one for $\bar{\nu}_{e}$ is
$P(\bar{\nu}_{e}\to\bar{\nu}_{e})=0$. Note that when we report our projected
constraints on the neutrino parameters below, neutrino flavors always refer to
the those emitted in the supernova explosion, i.e., prior to any flavor
oscillations.

##  4 Results

Only DSNB (HyperK + DUNE), no oscillations

Figure 2: Posterior probability distribution functions of the supernova
neutrino parameters obtained by observing the DSNB in DUNE and Hyper-
Kamiokande, assuming no oscillations. Left and right panels refer to average
and total energy of $\nu_{e}$ and $\bar{\nu}_{e}$, respectively.

No oscillations, $C^{238}=10^{-11}$

Figure 3: Posterior probability distribution functions of the supernova
neutrino parameters obtained assuming no oscillations and $C^{238}=10^{-11}$,
using both paleo detectors and DSNB observations in Hyper-Kamiokande and DUNE.
Left, central and right panels refer to average and total energy of $\nu_{e}$,
$\bar{\nu}_{e}$ and $\nu_{x}$, respectively.

Figure 2 displays the allowed regions at 68% and 99.7% confidence levels for
the average energy and total energy of $\nu_{e}$ and $\bar{\nu}_{e}$ assuming
no oscillations, obtained using only the measurement of DSNB neutrinos in DUNE
and Hyper-Kamiokande for 20-years of operation. A degeneracy between average
energy and total energy is observed, which is due to the fact that the
statistics of the signal events is proportional to the product $\langle
E_{\nu_{\alpha}}\rangle^{2}E_{\nu_{\alpha}}^{\rm tot}$. A better precision is
reached for Hyper-Kamiokande (right panel) due to its larger statistics.
Figure 3 displays the allowed regions at 68 and 99.7% confidence levels for
the average energy and total energy of $\nu_{e}$, $\bar{\nu}_{e}$ and
$\nu_{x}$, assuming $C^{238}=10^{-11}$ and no oscillations, obtained using
both paleo detectors and DSNB neutrinos in Hyper-Kamiokande and DUNE. The main
result here is that the best precision is obtained for the $\nu_{x}$
parameters since the
$\nu_{x}=\\{\nu_{\mu},\bar{\nu}_{\mu},\nu_{\tau},\bar{\nu}_{\tau}\\}$
contribute roughly four times more to the total neutrino flux than $\nu_{e}$
or $\bar{\nu}_{e}$, see eq. (2.2). In particular, the values $E^{\rm
tot}_{\nu_{x}}=0$ and $\langle E_{\nu_{x}}\rangle=0$ are disfavored at more
than 99.7% confidence level. Moreover, the constraints on the parameters for
the electron flavors are slightly improved with respect to the case with only
DSNB neutrinos. Figure 4 is equivalent to figure 3, but it includes
oscillation effects for both normal and inverted mass ordering and, for
brevity, we show only the results for the constraints on the $\nu_{x}$
parameters. For both mass orderings, the bounds on $\nu_{x}$ parameters are
slightly improved with respect to the no oscillation case. Although
oscillations do not affect the neutral current coherent scattering signal in
paleo detectors, oscillations do lead to signals from the $\nu_{x}$ component
in either Hyper-Kamiokande or DUNE, depending on the mass ordering under
consideration. The best precision is reached for inverted ordering, since in
this case, $\nu_{x}$ can be observed in Hyper-Kamiokande, which has a larger
statistics compared to DUNE.

Normal ordering, $C^{238}=10^{-11}$ Inverted ordering, $C^{238}=10^{-11}$

Figure 4: Same as figure 3, but for normal ordering (left panel) and inverted
ordering (right panel) of $\nu_{x}$.

No oscillations, $C^{238}=10^{-12}$ Normal ordering, $C^{238}=10^{-12}$
Inverted ordering, $C^{238}=10^{-12}$

Figure 5: Posterior probability distribution functions of the supernova
neutrino parameters obtained assuming no $C^{238}=10^{-12}$. Left, central and
right panels refer to the case of no oscillations, normal ordering and
inverted case of $\nu_{x}$, respectively.

In figure 5 we show the projected constraints on the $\nu_{x}$ parameters for
the three different assumptions of the flavor oscillations considered above
but for a more optimistic 238U concentration of $C^{238}=10^{-12}\,$g/g in the
paleo detectors rather than the $C^{238}=10^{-11}\,$g/g assumed in figures 3
and 4. Comparing these figures we can see that using paleo detectors with
uranium concentration of $C^{238}=10^{-12}\,$g/g would allow to constrain
$\langle E_{\nu_{i}}\rangle$ and $E_{\nu_{i}}^{\rm tot}$ with (relative)
precision $\lesssim 10\%$ rather than the few tens of percent precision we
project for $C^{238}=10^{-11}\,$g/g. This improvement is due to the reduction
of the dominant background in paleo detectors, the (elastic) scattering of
neutrons originating from the 238U decay chain.

If we modify the assumption of a 10% prior on both the galactic supernova rate
and the normalization of the DSNB, our results are significantly affected (not
shown). In particular, when increasing the width of these priors to 20%, the
99.7% confidence level region extends towards $\langle E_{\nu_{x}}\rangle=0$
and large values of $E_{\rm tot}^{\nu_{x}}$. This happens because widening
this prior loosens the link between the DSNB measurements in Hyper-Kamiokande
and DUNE with the paleo-detector measurement of neutrinos from galactic core-
collapse events, making it possible for the contribution of $\nu_{x}$ to be
compensated by $\nu_{e}$ and $\bar{\nu}_{e}$. On the other hand, increasing
the uncertainty on the normalization of backgrounds for DSNB experiments does
not appreciably affect the final constraints on $\langle E_{\nu_{x}}\rangle$
and $E_{\nu_{x}}^{\rm tot}$.

A final remark is in order. Recently, the possibility of measuring the DSNB in
conventional dark matter detectors and extracting the $\nu_{x}$ parameters has
been carefully studied in Ref. [24]. Analogously to the signal in paleo
detectors, in conventional dark matter detectors the supernova neutrino signal
is dominated by background, but the signal-to-noise ratio is smaller for two
reasons: Conventional dark matter detectors have much smaller exposures than
paleo detectors (e.g, the DARWIN proposal has an envisaged exposure on the
order of $\varepsilon\sim 100\,$t yr [23]). Second, conventional dark matter
detectors would be sensitive to DSNB neutrinos, while paleo detectors could
measure the (time-averaged) signal from galactic supernova neutrinos; the
time-averaged galactic supernova flux is roughly two orders of magnitude
larger than the DNSB flux. This difference in the signal-to-noise ratio is the
reason why in Ref. [24], only prospective upper limits on $\nu_{x}$ parameters
are provided, whereas in our work, we show goodness-of–fit regions.
Nevertheless, the two approaches to measure the $\nu_{x}$ component are
complementary. Paleo detectors allow only to probe galactic supernovae,
whereas dark matter detectors can measure the DSNB, which includes
extragalactic supernova neutrinos. Furthermore, neutrinos propagating over
cosmological distances, like those of the DSNB, allow studying different
fundamental physics (see, for example, Refs. [86, 87]).

##  5 Conclusions

Supernova neutrinos carry invaluable information about the astrophysical
processes occurring deep inside a stellar core collapse. The few tens of
events observed from SN1987a have already demonstrated the potential of
detecting a supernova neutrino burst, despite the small statistics collected.
Next generation experiments, such as Hyper-Kamiokande, DUNE and JUNO, would
detect $\mathcal{O}(10^{5})$ $\bar{\nu}_{e}$, $\mathcal{O}(10^{3})$ $\nu_{e}$
and $\mathcal{O}(10^{2})$ $\nu_{x}$ events, respectively, from a galactic
supernova, thus covering all neutrino flavors with differing precision.
However, only a few supernovae are expected to occur per century in our galaxy
on average, and the only remedy is patience. A possible way around this
problem is probing the DSNB, i.e., the integrated flux of neutrinos produced
by all the supernova explosions occurring throughout the Universe. Another
benefit of the DSNB is that it allows one to be sensitive to the entire
population of supernovae and not only to a specific progenitor, which might
have peculiar characteristics. Hyper-Kamiokande, DUNE and JUNO have the
capability to detect the DSNB, but with much smaller statistics compared to a
single galactic supernova and with a significant background to deal with. For
$\bar{\nu}_{e}$ events in Hyper-Kamiokande, such background can be mitigated
by adding gadolinium in the water Cherenkov detector, which can provide a 90%
efficiency in tagging inverse beta decay events induced by neutrinos. For
$\nu_{x}$, the detection prospects are less promising, even considering dark
matter detectors where the background is intrinsically lower. For instance,
the expected number of DSNB events in XENONnT is $\mathcal{O}(10^{-3})$ per
year [24], whereas the background is expected to add up to a few events per
year [88]. In this case only upper limits can be set on the parameters
describing the $\nu_{x}$ flux, i.e., the average energy and total energy,
$\langle E_{\nu_{x}}\rangle$ and $E_{\nu_{x}}^{\rm tot}$ [24].

In this work we have shown that one could measure the $\nu_{x}$ parameters
using paleo-detectors. Differently from conventional detectors which measure
neutrino events in real time, paleo detectors integrate events over time
scales as large as a billion years. Thus, paleo detectors could measure the
integrated neutrino flux from core-collapse supernovae in our Galaxy via the
nuclear recoils induced by (flavor-blind) CE$\nu$NS reactions of neutrinos
with the atomic nuclei in paleo detectors. Leveraging the gigayear exposure
times, with only 100 g of a epsomite paleo detector, we expect about 3$\times
10^{3}$ $\nu_{x}$ events and a background of $2\times 10^{5}$ events, assuming
a uranium-238 concentration of $C^{238}=10^{-11}\,$g/g. Although the signal-
to-background ratio is similar to what is expected for DSNB events in dark
matter detectors, the statistics is orders of magnitude higher, leading to a
much better signal-to-noise ratio – the ratio between signal events and the
square root of background is 6.6$\sigma$. Combining such a measurement of the
galactic supernova neutrino flux with measurements of the DSNB via charged
current interactions in Hyper-Kamiokande and DUNE would allow to provide
closed contours of the 68% and 99.7% confidence level regions for the
parameters controlling the supernova neutrino spectra of all flavors,
$\nu_{e}$, $\overline{\nu}_{e}$, and $\nu_{x}$. With a more optimistic
assumption on the concentration of uranium, $C^{238}=10^{-12}\,$g/g, we
project that the parameters $\langle E_{\nu_{x}}\rangle$ and $E_{\nu_{x}}^{\rm
tot}$ could be measured with a precision of $\lesssim 10\%$ on each
parameters. However, the precision on some of these supernova neutrino
parameters can be lost for uncertainties on the DSNB and galactic supnernova
rate larger than the 10% assumed here. Current uncertainties on the cosmic
star formation rate at low redshifts are already at the level of $\sim 20\%$
[42, 21, 43], and while reaching $\sim 10\%$ uncertainties is ambitious,
upcoming surveys with, for example, the James Webb Space Telescope or the Vera
C. Rubin Observatory are expected to lead to considerable improvements in our
knowledge of the cosmic and local star formation rate in the near future. We
emphasize that the use of paleo detectors is complementary to the one of dark
matter detectors, since they are sensitive to a different population of
supernova progenitors and potentially to different effects of physics beyond
the Standard Model.

## Acknowledgements

SB is supported in part by NSF Grant PHY-2014215, DOE HEP QuantISED award
#100495, and the Gordon and Betty Moore Foundation Grant GBMF7946. The work of
FC at Virginia Tech is supported by the U.S. Department of Energy under the
award number DE-SC0020250 and DE-SC0020262. The work of FC at IFIC is
supported by GVA Grant No.CDEIGENT/2020/003. The work of SH is supported by
the US Department of Energy under the award number DE-SC0020262 and NSF Grant
numbers AST-1908960 and PHY-1914409. This work was supported by World Premier
International Research Center Initiative (WPI Initiative), MEXT, Japan.

## References

  * [1] Kamiokande-II Collaboration, K. Hirata et al., _Observation of a Neutrino Burst from the Supernova SN 1987a_ , _Phys. Rev. Lett._ 58 (1987) 1490.
  * [2] R. M. Bionta et al., _Observation of a Neutrino Burst in Coincidence with Supernova SN 1987a in the Large Magellanic Cloud_ , _Phys. Rev. Lett._ 58 (1987) 1494.
  * [3] M. Aglietta et al., _Supernova 1987A in the Large Magellanic Cloud_ , _IAU Circular_ No. 4323 (1987) .
  * [4] E. N. Alekseev, L. N. Alekseeva, I. V. Krivosheina and V. I. Volchenko, _Detection of the Neutrino Signal From SN1987A in the LMC Using the Inr Baksan Underground Scintillation Telescope_ , _Phys. Lett._ B205 (1988) 209.
  * [5] S. Horiuchi and J. P. Kneller, _What can be learned from a future supernova neutrino detection?_ , _J. Phys. G_ 45 (2018) 043002 [1709.01515].
  * [6] Hyper-Kamiokande Collaboration, K. Abe et al., _Hyper-Kamiokande Design Report_ , 1805.04163.
  * [7] DUNE Collaboration, B. Abi et al., _Deep Underground Neutrino Experiment (DUNE), Far Detector Technical Design Report, Volume II: DUNE Physics_ , 2002.03005.
  * [8] J. F. Beacom, W. M. Farr and P. Vogel, _Detection of Supernova Neutrinos by Neutrino Proton Elastic Scattering_ , _Phys. Rev. D_ 66 (2002) 033001 [hep-ph/0205220].
  * [9] B. Dasgupta and J. F. Beacom, _Reconstruction of supernova $\nu_{\mu}$, $\nu_{\tau}$, anti-$\nu_{\mu}$, and anti-$\nu_{\tau}$ neutrino spectra at scintillator detectors_, _Phys. Rev. D_ 83 (2011) 113006 [1103.2768].
  * [10] H.-L. Li, Y.-F. Li, M. Wang, L.-J. Wen and S. Zhou, _Towards a complete reconstruction of supernova neutrino spectra in future large liquid-scintillator detectors_ , _Phys. Rev. D_ 97 (2018) 063014 [1712.06985].
  * [11] H.-L. Li, X. Huang, Y.-F. Li, L.-J. Wen and S. Zhou, _Model-independent approach to the reconstruction of multiflavor supernova neutrino energy spectra_ , _Phys. Rev. D_ 99 (2019) 123009 [1903.04781].
  * [12] B. Chauhan, B. Dasgupta and V. Datar, _A deuterated liquid scintillator for supernova neutrino detection_ , _JCAP_ 11 (2021) 005 [2106.10927].
  * [13] C. J. Horowitz, K. J. Coakley and D. N. McKinsey, _Supernova observation via neutrino - nucleus elastic scattering in the CLEAN detector_ , _Phys. Rev. D_ 68 (2003) 023005 [astro-ph/0302071].
  * [14] R. F. Lang, C. McCabe, S. Reichard, M. Selvi and I. Tamborra, _Supernova neutrino physics with xenon dark matter detectors: A timely perspective_ , _Phys. Rev. D_ 94 (2016) 103009 [1606.09243].
  * [15] E. Cappellaro, R. Barbon and M. Turatto, _Supernova statistics_ , _Springer Proc. Phys._ 99 (2005) 347 [astro-ph/0310859].
  * [16] R. Diehl et al., _Radioactive Al-26 and massive stars in the galaxy_ , _Nature_ 439 (2006) 45 [astro-ph/0601015].
  * [17] A. Strumia and F. Vissani, _Neutrino masses and mixings and…_ , hep-ph/0606054.
  * [18] W. Li, R. Chornock, J. Leaman, A. V. Filippenko, D. Poznanski, X. Wang, M. Ganeshalingam and F. Mannucci, _Nearby supernova rates from the Lick Observatory Supernova Search - III. The rate-size relation, and the rates as a function of galaxy Hubble type and colour_ , _Mon. Not. Roy. Astron. Soc._ 412 (2011) 1473 [1006.4613].
  * [19] M. T. Botticella, S. J. Smartt, R. C. Kennicutt, Jr., E. Cappellaro, M. Sereno and J. C. Lee, _A comparison between star formation rate diagnostics and rate of core collapse supernovae within 11 Mpc_ , _Astron. Astrophys._ 537 (2012) A132 [1111.1692].
  * [20] S. M. Adams, C. S. Kochanek, J. F. Beacom, M. R. Vagins and K. Z. Stanek, _Observing the Next Galactic Supernova_ , _Astrophys. J._ 778 (2013) 164 [1306.0559].
  * [21] J. F. Beacom, _The Diffuse Supernova Neutrino Background_ , _Ann. Rev. Nucl. Part. Sci._ 60 (2010) 439 [1004.3311].
  * [22] C. Lunardini, _Diffuse supernova neutrinos at underground laboratories_ , _Astropart. Phys._ 79 (2016) 49 [1007.3252].
  * [23] DARWIN Collaboration, J. Aalbers et al., _DARWIN: towards the ultimate dark matter detector_ , _JCAP_ 11 (2016) 017 [1606.07001].
  * [24] A. M. Suliga, J. F. Beacom and I. Tamborra, _Towards probing the diffuse supernova neutrino background in all flavors_ , _Phys. Rev. D_ 105 (2022) 043008 [2112.09168].
  * [25] S. Baum, A. K. Drukier, K. Freese, M. Górski and P. Stengel, _Searching for Dark Matter with Paleo-Detectors_ , _Phys. Lett. B_ 803 (2020) 135325 [1806.05991].
  * [26] A. K. Drukier, S. Baum, K. Freese, M. Górski and P. Stengel, _Paleo-detectors: Searching for Dark Matter with Ancient Minerals_ , _Phys. Rev. D_ 99 (2019) 043014 [1811.06844].
  * [27] S. Baum, T. D. P. Edwards, K. Freese and P. Stengel, _New Projections for Dark Matter Searches with Paleo-Detectors_ , _Instruments_ 5 (2021) 21 [2106.06559].
  * [28] S. Baum, T. D. P. Edwards, B. J. Kavanagh, P. Stengel, A. K. Drukier, K. Freese, M. Górski and C. Weniger, _Paleodetectors for Galactic supernova neutrinos_ , _Phys. Rev. D_ 101 (2020) 103017 [1906.05800].
  * [29] R. L. Fleischer, P. B. Price, R. M. Walker and E. L. Hubbard, _Track registration in various solid-state nuclear track detectors_ , _Phys. Rev._ 133 (1964) A1443.
  * [30] R. L. Fleischer, P. B. Price and R. M. Walker, _Tracks of charged particles in solids_ , _Science_ 149 (1965) 383.
  * [31] R. L. Fleischer, P. B. Price and R. M. Walker, _Solid-state track detectors: applications to nuclear science and geophysics_ , _Ann. Rev. Nucl. Part. Sci._ 15 (1965) 1.
  * [32] S.-L. Guo, B.-L. Chen and S. Durrani, _Chapter 4 - solid-state nuclear track detectors_ , in _Handbook of Radioactivity Analysis (Third Edition)_ (M. F. L’Annunziata, ed.), pp. 233 – 298.  Academic Press, Amsterdam, third edition ed., 2012.  DOI.
  * [33] M. Rodriguez, W. Li, F. Chen, C. Trautmann, T. Bierschenk, B. Afra, D. Schauries, R. Ewing, S. Mudie and P. Kluth, _Saxs and tem investigation of ion tracks in neodymium-doped yttrium aluminium garnet_ , _Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms_ 326 (2014) 150 . 17th International Conference on Radiation Effects in Insulators (REI).
  * [34] F. Schaff, M. Bech, P. Zaslansky, C. Jud, M. Liebi, M. Guizar-Sicairos and F. Pfeiffer, _Six-dimensional real and reciprocal space small-angle x-ray scattering tomography_ , _Nature_ 527 (2015) 353.
  * [35] M. Holler, A. Diaz, M. Guizar-Sicairos, P. Karvinen, E. Färm, E. Härkönen, M. Ritala, A. Menzel, J. Raabe and O. Bunk, _X-ray ptychographic computed tomography at 16 nm isotropic 3d resolution_ , _Scientific Reports_ 4 (2014) 3857.
  * [36] J. R. Jordan, S. Baum, P. Stengel, A. Ferrari, M. C. Morone, P. Sala and J. Spitz, _Measuring Changes in the Atmospheric Neutrino Rate Over Gigayear Timescales_ , _Phys. Rev. Lett._ 125 (2020) 231802 [2004.08394].
  * [37] J. I. Collar, _Comments on ’limits on dark matter using ancient mica’_ , _Phys. Rev. Lett._ 76 (1996) 331 [astro-ph/9511055].
  * [38] N. Tapia-Arellano and S. Horiuchi, _Measuring solar neutrinos over gigayear timescales with paleo detectors_ , _Phys. Rev. D_ 103 (2021) 123016 [2102.01755].
  * [39] M. T. Keil, G. G. Raffelt and H.-T. Janka, _Monte Carlo study of supernova neutrino spectra formation_ , _Astrophys. J._ 590 (2003) 971 [astro-ph/0208035].
  * [40] S. Horiuchi, J. F. Beacom and E. Dwek, _The Diffuse Supernova Neutrino Background is detectable in Super-Kamiokande_ , _Phys. Rev. D_ 79 (2009) 083013 [0812.3157].
  * [41] S. Horiuchi, K. Sumiyoshi, K. Nakamura, T. Fischer, A. Summa, T. Takiwaki, H.-T. Janka and K. Kotake, _Diffuse supernova neutrino background from extensive core-collapse simulations of $8$-$100{\rm M}_{\odot}$ progenitors_, _Mon. Not. Roy. Astron. Soc._ 475 (2018) 1363 [1709.06567].
  * [42] A. M. Hopkins and J. F. Beacom, _On the normalisation of the cosmic star formation history_ , _Astrophys. J._ 651 (2006) 142 [astro-ph/0601463].
  * [43] P. Madau and M. Dickinson, _Cosmic Star Formation History_ , _Ann. Rev. Astron. Astrophys._ 52 (2014) 415 [1403.0007].
  * [44] S. Horiuchi, T. Kinugawa, T. Takiwaki, K. Takahashi and K. Kotake, _Impact of binary interactions on the diffuse supernova neutrino background_ , _Phys. Rev. D_ 103 (2021) 043003 [2012.08524].
  * [45] L.-G. Strolger, T. Dahlen, S. A. Rodney, O. Graur, A. G. Riess, C. McCully, S. Ravindranath, B. Mobasher and A. K. Shahady, _The Rate of Core Collapse Supernovae to Redshift 2.5 From The CANDELS and CLASH Supernova Surveys_ , _Astrophys. J._ 813 (2015) 93 [1509.06574].
  * [46] J. Billard, L. Strigari and E. Figueroa-Feliciano, _Implication of neutrino backgrounds on the reach of next generation dark matter direct detection experiments_ , _Phys. Rev. D_ 89 (2014) 023524 [1307.5458].
  * [47] C. A. J. O’Hare, _Dark matter astrophysical uncertainties and the neutrino floor_ , _Phys. Rev. D_ 94 (2016) 063527 [1604.03858].
  * [48] R. H. Helm, _Inelastic and Elastic Scattering of 187-Mev Electrons from Selected Even-Even Nuclei_ , _Phys. Rev._ 104 (1956) 1466.
  * [49] J. F. Ziegler, J. P. Biersack and U. Littmark, _The Stopping and Range of Ions in Solids_. Pergamon Press, New York, 01, 1985.
  * [50] J. F. Ziegler, M. D. Ziegler and J. P. Biersack, _SRIM – The stopping and range of ions in matter (2010)_ , _Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms_ 268 (2010) 1818 . 19th International Conference on Ion Beam Analysis.
  * [51] S. Thomson and G. Wardle, _Coloured natural rocksalts: a study of their helium contents, colours and impurities_ , _Geochimica et Cosmochimica Acta_ 5 (1954) 169 .
  * [52] K. C. Condie, C. S. Kuo, R. M. Walker and V. R. Murthy, _Uranium Distribution in Separated Clinopyroxenes from Four Eclogites_ , _Science_ 165 (1969) 57.
  * [53] J. A. S. Adams, J. K. Osmond and J. J. W. Rogers, _The geochemistry of thorium and uranium_ , _Phys. Chem. Earth_ 3 (1959) 298.
  * [54] M. Seitz and S. Hart, _Uranium and boron distributions in some oceanic ultramafic rocks_ , _Earth and Planetary Science Letters_ 21 (1973) 97 .
  * [55] W. E. Dean, _Section 5 trace and minor elements in evaporites_ , in _Marine Evaporites_.  SEPM Society for Sedimentary Geology, 1987.  DOI.
  * [56] M. Yui, Y. Kikawada, T. Oi, T. Honda and T. Nozaki, _Abundance of Uranium and Thorium in Rock Salts_ , _Radioisotopes_ 47 (1998) 488.
  * [57] W. Sanford, M. Doughten, T. Coplen, A. Hunt and T. Bullen, _Evidence for high salinity of Early Cretaceous sea water from the Chesapeake Bay crater_ , _Nature_ 503 (2013) 252.
  * [58] T. D. P. Edwards, B. J. Kavanagh, C. Weniger, S. Baum, A. K. Drukier, K. Freese, M. Górski and P. Stengel, _Digging for dark matter: Spectral analysis and discovery potential of paleo-detectors_ , _Phys. Rev. D_ 99 (2019) 043541 [1811.10549].
  * [59] _SOURCES 4A: A Code for Calculating (alpha,n), Spontaneous Fission, and Delayed Neutron Sources and Spectra_ , Tech. Rep. LA-13639-MS, Los Alamos National Lab, 1999.  10.2172/15215.
  * [60] A. J. Koning and D. Rochman, _Modern Nuclear Data Evaluation with the TALYS Code System_ , _Nucl. Data Sheets_ 113 (2012) 2841.
  * [61] D. Rochman, A. J. Koning, J. C. Sublet, M. Fleming et al., _The TENDL library: hope, reality and future_ , in _Proceedings of the International Conference on Nuclear Data for Science and Technology, September 11-16, 2016, Bruges, Belgium_ , 2016, https://tendl.web.psi.ch/bib_rochman/tendl.nd2016.pdf.
  * [62] J. C. Sublet, A. J. Koning, D. Rochman, M. Fleming and M. Gilbert, _TENDL-2015: Delivering Both Completeness and Robustness_ , in _Advances in Nuclear Nonproliferation Technology and Policy Conference, Sept. 25-30, Santa Fe, NM, USA_ , https://tendl.web.psi.ch/bib_rochman/ANTPC_TENDL.pdf.
  * [63] M. Fleming, J. C. Sublet, J. Kopecky, D. Rochman and A. J. Koning, _Probing experimental and systematic trends of the neutron-induced TENDL-2014 nuclear data library_ , in _CCFE report UKAEA-R(15)29, October 2015_ , 2015, https://tendl.web.psi.ch/bib_rochman/UKAEA-R(15)29_30.9.15.pdf.
  * [64] N. Soppera, M. Bossant and E. Dupont, _JANIS 4: An Improved Version of the NEA Java-based Nuclear Data Information System_ , _Nucl. Data Sheets_ 120 (2014) 294.
  * [65] A. Ferrari, P. R. Sala, A. Fasso and J. Ranft, _FLUKA: A multi-particle transport code (Program version 2005)_ , . CERN-2005-010, SLAC-R-773, INFN-TC-05-11.
  * [66] T. T. Böhlen, F. Cerutti, M. P. W. Chin, A. Fassò, A. Ferrari, P. G. Ortega, A. Mairani, P. R. Sala, G. Smirnov and V. Vlachoudis, _The FLUKA Code: Developments and Challenges for High Energy and Medical Applications_ , _Nucl. Data Sheets_ 120 (2014) 211.
  * [67] G. Battistoni, A. Ferrari, M. Lantz, P. R. Sala and G. I. Smirnov, _A neutrino-nucleon interaction generator for the fluka monte carlo code_ , in _CERN-Proceedings-2010-001_ , pp. 387–394, 2010.
  * [68] C. A. J. O’Hare, _Can we overcome the neutrino floor at high masses?_ , _Phys. Rev. D_ 102 (2020) 063024 [2002.07499].
  * [69] R. Bollig, N. Yadav, D. Kresse, H. T. Janka, B. Müller and A. Heger, _Self-consistent 3D Supernova Models From $-$7 Minutes to +7 s: A 1-bethe Explosion of a $\sim$19 $M_{\odot}$ Progenitor_, _Astrophys. J._ 915 (2021) 28 [2010.10506].
  * [70] DUNE Collaboration, B. Abi et al., _Deep Underground Neutrino Experiment (DUNE), Far Detector Technical Design Report, Volume I Introduction to DUNE_ , _JINST_ 15 (2020) T08008 [2002.02967].
  * [71] DUNE Collaboration, B. Abi et al., _Deep Underground Neutrino Experiment (DUNE), Far Detector Technical Design Report, Volume III: DUNE Far Detector Technical Coordination_ , _JINST_ 15 (2020) T08009 [2002.03008].
  * [72] DUNE Collaboration, B. Abi et al., _Deep Underground Neutrino Experiment (DUNE), Far Detector Technical Design Report, Volume IV: Far Detector Single-phase Technology_ , _JINST_ 15 (2020) T08010 [2002.03010].
  * [73] E. Kolbe, K. Langanke, G. Martinez-Pinedo and P. Vogel, _Neutrino nucleus reactions and nuclear structure_ , _J. Phys. G_ 29 (2003) 2569 [nucl-th/0311022].
  * [74] Z. Tabrizi and S. Horiuchi, _Flavor Triangle of the Diffuse Supernova Neutrino Background_ , _JCAP_ 05 (2021) 011 [2011.10933].
  * [75] P. Vogel and J. F. Beacom, _Angular distribution of neutron inverse beta decay, anti-neutrino(e) + p — $>$ e+ + n_, _Phys. Rev. D_ 60 (1999) 053003 [hep-ph/9903554].
  * [76] A. Strumia and F. Vissani, _Precise quasielastic neutrino/nucleon cross-section_ , _Phys. Lett. B_ 564 (2003) 42 [astro-ph/0302055].
  * [77] JUNO Collaboration, Z. Djurcic et al., _JUNO Conceptual Design Report_ , 1508.07166.
  * [78] K. Møller, A. M. Suliga, I. Tamborra and P. B. Denton, _Measuring the supernova unknowns at the next-generation neutrino telescopes through the diffuse neutrino background_ , _JCAP_ 05 (2018) 066 [1804.03157].
  * [79] J. P. Gardner et al., _The James Webb Space Telescope_ , _Space Sci. Rev._ 123 (2006) 485 [astro-ph/0606175].
  * [80] LSST Science Collaboration Collaboration, P. A. Abell et al., _LSST Science Book, Version 2.0_ , 0912.0201.
  * [81] O. Snaith, M. Haywood, P. Di Matteo, M. D. Lehnert, F. Combes, D. Katz and A. Gómez, _Reconstructing the star formation history of the Milky Way disc(s) from chemical abundances_ , _Astron.Astrophys_ 578 (2015) A87 [1410.3829].
  * [82] M. Haywood, M. D. Lehnert, P. Di Matteo, O. Snaith, M. Schultheis, D. Katz and A. Gómez, _When the Milky Way turned off the lights: APOGEE provides evidence of star formation quenching in our Galaxy_ , _Astron.Astrophys_ 589 (2016) A66 [1601.03042].
  * [83] R. Mor, A. C. Robin, F. Figueras, S. Roca-Fabregas and X. Luri, _Gaia DR2 reveals a star formation burst in the disc 2-3 Gyr ago_ , _Astron.Astrophys_ 624 (2019) L1 [1901.07564].
  * [84] F. Feroz, M. P. Hobson and M. Bridges, _MultiNest: an efficient and robust Bayesian inference tool for cosmology and particle physics_ , _Mon. Not. Roy. Astron. Soc._ 398 (2009) 1601 [0809.3437].
  * [85] J. Buchner, A. Georgakakis, K. Nandra, L. Hsu, C. Rangel, M. Brightman, A. Merloni, M. Salvato, J. Donley and D. Kocevski, _X-ray spectral modelling of the AGN obscuring region in the CDFS: Bayesian model selection and catalogue_ , _Astron. Astrophys._ 564 (2014) A125 [1402.0004].
  * [86] Y. S. Jeong, S. Palomares-Ruiz, M. H. Reno and I. Sarcevic, _Probing secret interactions of eV-scale sterile neutrinos with the diffuse supernova neutrino background_ , _JCAP_ 06 (2018) 019 [1803.04541].
  * [87] A. De Gouvêa, I. Martinez-Soler, Y. F. Perez-Gonzalez and M. Sen, _Fundamental physics with the diffuse supernova background neutrinos_ , _Phys. Rev. D_ 102 (2020) 123012 [2007.13748].
  * [88] XENON Collaboration, E. Aprile et al., _Projected WIMP sensitivity of the XENONnT dark matter experiment_ , _JCAP_ 11 (2020) 031 [2007.08796].


#  Ab initio determination of pseudospin for paramagnetic defects in SiC

András Csóré  Department of Atomic Physics, Budapest University of Technology
and Economics, Budafoki út 8., H-1111, Budapest, Hungary  Adam Gali  Wigner
Research Centre for Physics, PO. Box 49, Budapest H-1525, Hungary  Department
of Atomic Physics, Budapest University of Technology and Economics, Budafoki
út 8., H-1111, Budapest, Hungary

(August 27, 2024)

###### Abstract

Paramagnetic point defects in solids may exhibit a rich set of interesting and
not yet fully resolved physics. In particular, character of wavefunctions and
electron-phonon coupling in these defects may highly influence their
interaction with external magnetic fields. Complex interplay between the
electronic orbitals, phonons and electron spin determines the effective
pseudospin of the system that we demonstrate on vanadium and molybdenum
defects in hexagonal silicon carbide (SiC) by means of ab initio calculations.
In this Letter, we find a giant anisotropy in the $g$-tensor of these defects
with Kramers doublet spin ground state, resulting in reduced and vanishing
interaction with the magnetic field in parallel and transverse directions,
respectively. The consequences of our finding in the application of these
defects for quantum information processing are briefly discussed.

Point defects may introduce levels in the fundamental band gap of
semiconductors or insulators that radically change the optical and magnetic
properties of the host material. In particular, these point defects could be
paramagnetic, i.e. the electron spin is greater than zero. The simplest case
introducing paramagnetic electronic structure is the Kramers doublet (KD)
electron spin state. In particular, neutral vanadium (V) Kaufmann  _et al._
(1997); Kunzer  _et al._ (1993); Spindlberger  _et al._ (2019) and singly
positive charge state of molybdenum (Mo) point defects Csóré  _et al._ (2016);
Gällström  _et al._ (2009); Bosma  _et al._ (2018) in hexagonal SiC are
identifed as KD systems. Recently, highly anisotropic interaction with
magnetic field governed by the corresponding $g$-tensors has been observed or
tentatively proposed for these KD systems Bosma  _et al._ (2018); Baur  _et
al._ (1997), however the underlying physics has not been understood so far.
Here we note that, although parallel component ($g_{\parallel}$) of $g$-tensor
was thoroughly studied for V in 6H SiC by invoking crystal field theory, the
emerging strong electron-phonon coupling is not included Kaufmann  _et al._
(1995); Dörnen  _et al._ (1992); Maier  _et al._ (1992); Reinke  _et al._
(1993); Kunzer  _et al._ (1993) or considered as minor effect Kaufmann  _et
al._ (1997). The strong anisotropy in $g$-tensor can be described by an
effective Hamiltonian and pseudospins in which parameters cannot be predicted
by applying simple models, and rather fitting procedure to known experimental
data is applied. In order to identify and understand these spin-related
phenomena and correct interpretation of experimental data, there is an urgent
need to determine these spin-Hamiltonian parameters from first principles
calculations.

In this Letter, we demonstrate that ab initio calculations can reproduce the
anisotropy of the observed $g$-tensors for V and Mo defects in hexagonal SiC.
We show that both electron-phonon coupling manifested as dynamic Jahn-Teller
(DJT) effect Bersuker (2006) and the character of the wavefunction will
determine the pseudospin of the system, i.e. its interaction with the external
magnetic field. We show that the complex interplay of electronic orbitals,
phonons and spins results in modified parallel ($g_{\parallel}$) component
with respect to the free electron $g$ factor ($g_{0}=2.003$) and vanishing
transverse ($g_{\perp}$) component. We discuss the relevance of our results in
the light of realization of telecom wavelength solid state qubits.

Both transition metal (TM) atoms substitute a Si atom in the SiC lattice as
found in earlier studies Gällström  _et al._ (2009); Ivády  _et al._ (2011).
However, lattice structures of 4H and 6H polytypes offer inequivalent lattice
sites implying TM${}_{\text{Si}}$ (TM = {Mo+,V}) defects to form two
configurations – a hexagonal ($h$) and a quasicubic ($k$) one – in 4H, and
three configurations – a hexagonal ($h$) and two quasicubic ($k_{1},k_{2}$)
ones – in 6H SiC [cf. Fig. 1(a)] all exhibiting C${}_{\text{3v}}$ symmetry.
For V${}_{\text{Si}}$ all configurations have been observed in both 4H
Spindlberger  _et al._ (2019) and 6H SiC Kaufmann  _et al._ (1997); Kunzer
_et al._ (1993); Kaufmann  _et al._ (1995); Dörnen  _et al._ (1992); Maier
_et al._ (1992); Reinke  _et al._ (1993). Recently, we have conclusively
identified V${}_{\text{Si}}$ configurations in 4H SiC Spindlberger  _et al._
(2019). Identification of V${}_{\text{Si}}$ configurations in 6H SiC is
possible via the corresponding spin-orbit (SO) splittings and $g$-factors that
is provided in this Letter. On the other hand, only a single signal has been
detected for Mo${}^{+}_{\text{Si}}$ in both hexagonal SiC polytypes as
reported in recent PL studies Gällström  _et al._ (2009); Bosma  _et al._
(2018). Since there is uncertainty with respect to the Mo${}^{+}_{\text{Si}}$
configuration being responsible for the single PL signal, we investigate both
defect models in 4H SiC, i.e. Mo${}^{+}_{\text{Si}}(h)$ and also
Mo${}^{+}_{\text{Si}}(k)$. Our detailed results on Mo${}^{+}_{\text{Si}}$
defect models are provided in the Supplementary Information Sup .

TM${}_{\text{Si}}$ defects were embedded in a 576-atom 4H supercell and a
432-atom 6H supercell. For sampling the Brillouin-zone we used $\Gamma$-point
which ensures the correct degeneracy of orbitals in C${}_{\text{3v}}$
symmetry. Plane wave expansion of Kohn-Sham wavefunctions with a cutoff of 420
eV was applied as a natural choice for supercell-method. Relaxed geometries
were achieved by minimizing the total energy with respect to the coordinates
of the ions where the corresponding quantum mechanical forces are prescribed
to fall below 0.01 eV/$\AA$. We treated the core electrons within the
framework of Projector Augmented Wave (PAW) method Blöchl (1994); Sup  as
implemented in the VASP code Kresse and Furthmüller (1996). In order to
compute the spin-orbit (SO) splitting in the ground state we employed
noncollinear approach Steiner  _et al._ (2016) with fixed spin quantization
axis along the crystal axis ($c$-axis), where the geometry was fixed in
C${}_{\text{3v}}$ configurations as obtained from spinpolarized calculations.
From SO calculations orbitally reduced angular momentum (discussed later on)
for each KS orbital can be directly read out Steiner  _et al._ (2016). We
employed density functional theory (DFT) to calculate the electronic structure
within the hybrid-DFT + Vw scheme introduced by Ivády et al. Heyd  _et al._
(2003); Ivády  _et al._ (2014). For the corresponding $w$ values in the ground
state we found $w_{\text{Mo}}\approx 0$ eV for Mo${}_{\text{Si}}^{+}$ Csóré
_et al._ (2016) and $w_{\text{V}}=2.2$ eV for V${}_{\text{Si}}$ Spindlberger
_et al._ (2019). We briefly describe this technique in Ref. Sup .

Both TM${}_{\text{Si}}$ defects introduce spin doublet ($S=\frac{1}{2}$), i.e.
KD ground state is formed by a single electron residing on a degenerate in-gap
$e$ level. In addition, higher-energy empty $a_{1}$ and $e$ levels also occur
in the band gap as plotted in Fig. 1(b), however their energy order is site
dependent, i.e. $a_{1}(0)e(0)$ for TM${}_{\text{Si}}(h)$ and $e(0)a_{1}(0)$
for TM${}_{\text{Si}}(k)$ as already reported in Ref. Spindlberger  _et al._ ,
2019. As a result the ground state all-electron wavefunction transforms as
${}^{2}E$. However, the lower degenerate $e$ level is split by the SO coupling
resulting in two SO sublevels in the ground state denoted by GS1 and GS2 in
energy order [cf. Fig. 1(c)]. All in-gap one-electron levels exhibit closely
atomic-like $d$-orbital character, thus the symmetry of the GS1-2 KD
wavefunctions may be determined by the atomic KD states for each
TM${}_{\text{Si}}$ as constructed from $d$-orbitals (see Table 1) Bosma  _et
al._ (2018). In particular, states can be established as proper linear
combinations of $\Psi_{1-4}$ (see Table 1) providing that GS1-2 transform
either as $E_{\frac{1}{2}}$ or $E_{\frac{3}{2}}$ representations in the
C${}_{\text{3v}}$ double group notation. Here we note that, a close inspection
on the local environment of the different sites implies that the $h$ site
exhibits the weaker C${}_{\text{3v}}$ (and thus stronger T${}_{\text{d}}$)
character than the $k$ or $k_{1}/k_{2}$ sites, while $k$ and $k_{2}$ show the
strongest C${}_{\text{3v}}$ nature. Consequently, SO splitting is expected to
be the lowest for TM${}_{\text{Si}}(h)$ defects, since in T${}_{\text{d}}$
symmetry it is entirely quenched in the first order based on group theory
considerations.

Figure 1: (a) Employed 4H (left) and 6H (right) SiC supercells in orthographic view embedding a TM atom at a Si(h) site. Crystal axis ($c$-axis), bilayer structures and labeling of atoms are indicated. (b) One-electron structure and (c) splittings in $d$-orbitals of TM${}_{\text{Si}}(h)$ defects in both 4H and 6H SiC. In the level structure of TM${}_{\text{Si}}(k)$ defects order of $a_{1}^{*}$ and $e^{*}$ levels are reversed and hence ${}^{2}A_{1}$ and the upper ${}^{2}E$ levels are also swapped. Energy of SO splittings in the ground and excited states are denoted by $E_{\text{SO}}^{\text{GS}}$ and $E_{\text{SO}}^{\text{ES}}$, respectively. Table 1: Kramers doublets formed by $d$-orbitals and the corresponding single and double group irreducible representations under C${}_{\text{3v}}$ symmetry. We give widespread notations for double group irreducible representations (irreps) and also the corresponding $m_{j}$ values. labels | orbitals | irreps. | $m_{j}$
---|---|---|---
single | double
$\Psi_{1}$ | $\ket{{d_{+2},+\frac{1}{2}}};\ket{{d_{-2},-\frac{1}{2}}}$ | ${}^{2}E$ |  $E_{\frac{1}{2}}$ ($\Gamma_{4}$) | $\pm\frac{5}{2}$
$\Psi_{2}$ | $\ket{{d_{+2},-\frac{1}{2}}};\ket{{d_{-2},+\frac{1}{2}}}$ | ${}^{2}E$ |  $E_{\frac{3}{2}}$ ($\Gamma_{5,6}$) | $\pm\frac{3}{2}$
$\Psi_{3}$ | $\ket{{d_{+1},+\frac{1}{2}}};\ket{{d_{-1},-\frac{1}{2}}}$ | ${}^{2}E$ |  $E_{\frac{3}{2}}$ ($\Gamma_{5,6}$) | $\pm\frac{3}{2}$
$\Psi_{4}$ | $\ket{{d_{+1},-\frac{1}{2}}};\ket{{d_{-1},+\frac{1}{2}}}$ | ${}^{2}E$ |  $E_{\frac{1}{2}}$ ($\Gamma_{4}$) | $\pm\frac{1}{2}$
$\Psi_{5}$ | $\ket{{d_{0},+\frac{1}{2}}};\ket{{d_{0},-\frac{1}{2}}}$ | ${}^{2}A_{1}$ |  $E_{\frac{1}{2}}$ ($\Gamma_{4}$) | $\pm\frac{1}{2}$

In this work, we aim to reveal the pivotal role of shape of the wavefunction
and the interaction of the described electronic structure with the SiC phonon
bath resulting in the highly anisotropic $g$-tensors. To this end, we derive
the parallel and transverse elements of the diagonalized $g$-tensor
($\mathbf{g}$) from the following spin Hamiltonian including the SO coupling
($\hat{H}_{\text{SO}}$), the Zeeman effect ($\hat{H}_{\text{Zee}}$) upon
external magnetic field and the hyperfine (HF) interaction
($\hat{H}_{\text{HF}}$) as
$\hat{H}_{\text{eff}}=\hat{H}_{\text{SO}}+\hat{H}_{\text{Zee}}+\hat{H}_{\text{HF}}$,
where

$\displaystyle\hat{H}_{\text{SO}}=-p\lambda_{0}\hat{L_{z}}\hat{S_{z}},$
$\displaystyle\hat{H}_{\text{Zee}}=-\mathbf{B}\hat{\boldsymbol{\mu}},$
$\displaystyle\hat{H}_{\text{HF}}=\hat{\mathbf{I}}\mathbf{A}\hat{\mathbf{S}}\text{.}$
(1)

In the expression of $\hat{H}_{\text{SO}}$ (Eq. 1) $\hat{S}_{z}$ and
$\hat{L}_{z}$ are the $z$-components of the spin and angular momentum
operators, $\hat{\mathbf{S}}$ and $\hat{\mathbf{L}}$, respectively,
$\lambda_{0}$ is the SO constant and $p$ stands for the so-called Ham
reduction factor Ham (1968); Bersuker (2006) arising from the electron-phonon
coupling. HF interaction is definitely active because of the presence of the
impurity atom with non-zero nuclear spin. In the HF Hamiltonian
$\hat{\mathbf{I}}$ is the nuclear spin operator and $\mathbf{A}$ is the HF
tensor which can be separated into a Fermi-contact and a dipolar spin-spin
interaction where the latter can rotate the electron spin because of the spin
ladder operators ($S^{+}I^{-}+S^{-}I^{+}$). In the formula of
$\hat{H}_{\text{Zee}}$ (Eq. 1) $\mathbf{B}$ represents the external magnetic
field and the dipole momentum operator ($\hat{\boldsymbol{\mu}}$) can be
expressed as

$\hat{\boldsymbol{\mu}}=-(\mu_{\text{B}}pr\hat{\mathbf{L}}+\mu_{\text{B}}g_{0}\hat{\mathbf{S}})=-\mu_{\text{B}}\mathbf{g}\mathbf{\tilde{S}}\text{,}$
(2)

where the contributions of $\hat{\mathbf{L}}$ and $\hat{\mathbf{S}}$ are
separated Thiering and Gali (2018) and $\hat{\boldsymbol{\mu}}$ is also
expressed in the phenomenological pseudospin ($\mathbf{\tilde{S}}$) formalism
Abragam and Bleaney (1970), where $\mathbf{\tilde{S}}=\frac{1}{2}$ for KD
systems Chibotaru and Ungur (2012). In Eq. 2 $\mu_{\text{B}}$ is the Bohr
magneton, $g_{0}$ is the free electron spin $g$-factor and $r$ represents the
Stevens orbital reduction factor Stevens and Pryce (1953). From Eq. 2 parallel
($g_{\parallel}$) and transverse ($g_{\perp}$) components of the $g$-tensor
($\mathbf{g}$) can be expressed as

$\displaystyle g_{\parallel}$
$\displaystyle=2(g_{0}S_{z}+L_{z}^{\text{eff}})=\frac{2\mu_{z}}{\mu_{B}}\text{,}$
(3) $\displaystyle g_{\perp}$
$\displaystyle=\frac{\mu_{+}+\mu_{-}+i(\mu_{-}-\mu_{+})}{\mu_{B}}\text{,}$ (4)

where we use expectation values of the ladder dipole moment operators
($\mu_{\pm}$) to express $g_{\perp}$ (Eq. 4) Bosma  _et al._ (2018). In Eq. 3
$S_{z}$ and $L_{z}^{\text{eff}}$ are expectation values of $\hat{S}_{z}$ and
the effective angular momentum operator,
$\hat{L}_{z}^{\text{eff}}=pr\hat{L}_{z}$, respectively.

Foremost, we calculate $g_{\parallel}$ based on Eq. 3 via obtaining
$\hat{L}_{z}^{\text{eff}}$ arising from the double reduction of the atomic
angular moment, $L_{z}$. Eigenvalue ($L_{z}^{\text{o}}$) of the orbitally
reduced angular momentum operator $\hat{L}_{z}^{\text{o}}=r\hat{L}_{z}$ can be
directly read out from first principles calculations (see Table 2). To
determine the eigenvalue ($L_{z}^{\text{eff}}$) of the effective angular
moment operator $\hat{L}_{z}^{\text{eff}}=p\hat{L}_{z}^{\text{o}}$, we account
for the emerging strong electron-phonon coupling by calculating $p$. Finally,
we determine $g_{\perp}$ according to Eq. 4.

Table 2: Parameters of the corresponding quadratic DJT APES ($E_{\text{JT}}$, $\delta_{\text{JT}}$) allowing the calculation of the effective phonon energy ($\hbar\omega$) and Ham reduction factor ($p$). Intrinsic (reduced) SO splitting energies ($\Delta E_{\text{SOC}}^{\text{(red)}}$) are also presented. Orbitally reduced (effective) angular momenta ($L_{z}^{\text{o(eff)}}$) for GS1 and GS2 and difference between them $L_{z}^{\text{eff}}$ values ($\Delta L_{z}^{\text{eff}}$) are also provided for V${}_{\text{Si}}$ in 4H and 6H SiC. Polytype | Site | $E_{\text{JT}}$ | $\delta_{\text{JT}}$ | $\hbar\omega$ | $p$ | $\Delta E_{\text{SOC}}$ | $\Delta E_{\text{SOC}}^{\text{red}}$ | $\Delta E_{\text{SOC}}^{\text{exp}}$ | $L_{z}^{\text{o}}$ (GS1, GS2) | $L_{z}^{\text{eff}}$ (GS1, GS2) | $\Delta L_{z}^{\text{eff}}$
---|---|---|---|---|---|---|---|---|---|---|---
(meV) | (meV) | (meV) | (GHz) | (GHz) | (GHz)
4H | $h$ | 9.4 | 5.6 | 60.19 | 0.63 | 9.91 | 6.29 | 43 ($\beta$) | -0.022, -0.013 | -0.014,-0.009 | 0.005
$k$ | 13.1 | 7.1 | 49.81 | 0.60 | 819.21 | 490.37 | 529 ($\alpha$) | -0.125,0.094 | -0.069,0.059 | 0.128
6H | $h$ | 9.3 | 5.8 | 47.19 | 0.57 | 24.18 | 13.78 | 16 ($\gamma$) | -0.012,-0.018 | -0.007,-0.010 | 0.003
$k_{1}$ | 11.4 | 5.8 | 49.43 | 0.55 | 82.94 | 45.62 | 25 ($\beta$) | -0.015,-0.016 | -0.008,-0.009 | 0.001
$k_{2}$ | 11.9 | 6.1 | 65.91 | 0.61 | 808.58 | 493.23 | 524 ($\alpha$) | -0.117,0.086 | -0.071,0.052 | 0.123

Ground state electronic structure introduced by the TM${}_{\text{Si}}$ defects
[cf. Fig. 1(b)], i.e. the half-filled orbitally degenerate $e$ level may split
by coupling with $e$ phonon modes as manifestation of $E\otimes e$ JT effect
Bersuker (2006); Ham (1968). Experimental results imply no symmetry reduction
thus dynamic JT (DJT) is expected for these systems Spindlberger  _et al._
(2019); Bosma  _et al._ (2018). To estimate the magnitude of DJT effect
originating from the electron-phonon coupling we calculated the trivial points
of the quadratic DJT adiabatic potential energy surface (APES) Bersuker
(2006), i.e. the three minima (C${}_{\text{1h}}$), the three barrier
(C${}_{\text{1h}}$) and the high symmetry (C${}_{\text{3v}}$) points. Energy
separation between the C${}_{\text{3v}}$ and the three minima is the JT energy
($E_{\text{JT}}$), while barrier points are separated by the barrier energy
($\delta_{\text{JT}}$). A general quadratic DJT APES is shown in Fig. 2.

Figure 2: Schematic illustration of a general quadratic DJT APES, where
$Q_{x/y}$ denote the configurational coordinates. Trivial points and energy
separations between them are indicated. C${}_{\text{1h}}$(min) represents the
three global minima separated by three barriers denoted by C${}_{\text{1h}}$
(barrier).

As a result angular momentum might be severely reduced by the persisting DJT
effect known as Ham effect Ham (1968); Bersuker (2006) resulting in the
reduction of the spin-orbit coupling (SOC) and the $g$-tensor elements.
Reduction can be expressed as $p\hat{\mathbf{L}}$, where $p$ is the already
introduced Ham reduction factor. For calculation of $p$ corresponding APES
have to be determined as described by the Hamiltonian Bersuker (2006) of

$\displaystyle\hat{H}_{\text{DJT}}$
$\displaystyle=\hbar\omega(a_{x}^{{\dagger}}a_{x}+a_{y}^{{\dagger}}a_{y}+1)+F(x\sigma_{z}+y\sigma_{x})$
(5) $\displaystyle+G[(x^{2}-y^{2})\sigma_{z}+2xy\sigma_{x}],$

where $a_{x/y}^{({\dagger})}$ represent annihilation (creation) operators of
two-dimensional $e$ modes vibrating in the $xy$ plane and electrons are
represented by the Pauli matrices $\sigma_{x}$ and $\sigma_{z}$. In Eq. 5,
$\hbar\omega$ stands for the effective energy of the $e$ modes, while
$F=\sqrt{E_{\text{JT}}\cdot 2\hbar\omega}$ and
$G=\delta_{\text{JT}}\hbar\omega/2E_{\text{JT}}$ govern the linear and
quadratic nature of the APES, respectively. In this way all parameters in DJT
Hamiltonian can be directly readout from the corresponding APES enabling the
numerical solution of Eq. 5 and thus the determination of the polaronic
wavefunctions that can be expanded in complex basis as

$\ket{\Psi_{\pm}}=\sum_{nm}(c_{nm}\ket{E_{\pm}}\otimes\ket{n,m}+d_{nm}\ket{E_{\mp}}\otimes\ket{n,m})\text{.}$
(6)

In Eq. 6 $E_{\pm}$ represents the complex components of the ${}^{2}E$ ground
state, where the subscript denotes the corresponding $L_{z}=\pm 1$ formed by
two $d$-orbitals as $E_{-}$:{-2,+1} and $E_{+}$:{+2,-1}. States of $E_{\pm}$
are mixed by the vibronic wavefunctions of $\ket{n,m}$, where $n+m\leq 4$
basis set provides convergent $\ket{\Psi_{\pm}}$. In this way mixing
coefficients of $c_{nm}$ and $d_{nm}$ can be calculated that enables us to
determine $p$ via the formula of

$p=\sum_{nm}(c^{2}_{nm}-d^{2}_{nm})$ (7)

as derived and implemented by Thiering et al. Thiering and Gali (2018). The
corresponding calculated values of $E_{\text{JT}}$, $\delta_{\text{JT}}$ and
$p$ are listed in Table 2. We report the ground state polaronic wavefunctions
expanded in symmetry adapted basis for all vanadium defect configurations in
Ref. Sup .

We also report the corresponding intrinsic and reduced ground state SO
splittings ($E_{\text{SOC}}$) in Table 2 as obtained by our calculations.
Reduced SOC can be calculated as

$\displaystyle\Delta E_{\text{SOC}}^{\text{red}}$ $\displaystyle=p\Delta
E_{\text{SOC}}=p(E_{\text{SOC}}^{\text{GS2}}-E_{\text{SOC}}^{\text{GS1}})$ (8)
$\displaystyle=p(\lambda_{0}^{\text{GS2}}\langle\hat{L_{z}^{\text{o}}}\hat{S_{z}}\rangle^{\text{GS2}}-\lambda_{0}^{\text{GS1}}\langle\hat{L_{z}^{\text{o}}}\hat{S_{z}}\rangle^{\text{GS1}}),$

where $\Delta E_{\text{SOC}}^{\text{(red)}}$ stands for the intrinsic
(reduced) SOC splitting of the ground state, $\lambda_{z}^{\text{GS1-2}}$ is
the intrinsic SOC constant of the SO sublevels GS1-2 and
$\langle\hat{L_{z}^{\text{o}}}\hat{S_{z}}\rangle^{\text{GS1-2}}$ represents
the expectation value of the $\hat{L_{z}^{\text{o}}}\hat{S_{z}}$ product. In
the calculations we employ $S_{z}=+\frac{1}{2}$ for both GS1 and GS2. Since
the $d$-orbitals are well-localized on the V impurity and almost the entire SO
splitting originates from the V atom one may deduce from the atomic SO
splitting formula that $\Delta E_{\text{SOC}}^{\text{red}}\leavevmode\nobreak\
\sim\leavevmode\nobreak\ \Delta L_{z}^{\text{eff}}$, where $\Delta
L_{z}^{\text{eff}}\leavevmode\nobreak\ =\leavevmode\nobreak\
p(L_{z}^{\text{o,GS2}}\leavevmode\nobreak\ -\leavevmode\nobreak\
L_{z}^{\text{o,GS1}})$ Sup .

We list the calculated values of $E_{\text{SOC}}$ and
$E_{\text{SOC}}^{\text{red}}$ along with the experimental values Wolfowicz
_et al._ (2020) in Table 2. Accordingly, both calculated and experimental
values for V${}_{\text{Si}}(k)$ in 4H SiC and for V${}_{\text{Si}}(k_{2})$ in
6H SiC are at least one order of magnitude higher than that for
V${}_{\text{Si}}(h)$ in 4H and V${}_{\text{Si}}(h/k_{1})$ in 6H SiC. This
significant difference may be explained via the calculated $\Delta
L_{z}^{\text{eff}}$ values also listed in Table 2 where similar trend occurs
obeying the previosuly deduced linear relationship between $\Delta
E_{\text{SOC}}^{\text{red}}$ and $\Delta L_{z}^{\text{eff}}$. Furthermore, SO
splitting for V${}_{\text{Si}}(k_{1})$ is slightly larger than that of
V${}_{\text{Si}}(h)$ arising from the local symmetry exhibiting stronger
C${}_{\text{3v}}$ character for V${}_{\text{Si}}(k_{1})$ than
V${}_{\text{Si}}(h)$.

At this point we are ready to calculate $g_{\parallel}$ according to Eq. 3:
values are included in Table 3. We find good agreement between the trends of
the calculated and that of the experimental $g_{\parallel}$ values. In
particular, $g_{\parallel}$ for V${}_{\text{Si}}(k)$ in 4H and for
V${}_{\text{Si}}(k_{2})$ in 6H SiC are well-separated from those of
V${}_{\text{Si}}(h)$ in 4H and V${}_{\text{Si}}(h/k_{1})$ in 6H SiC,
respectively, supporting our defect identification based on the SO splittings.
For the calculation of $g_{\perp}$ (Eq. 4), we consider the ladder magnetic
dipole operator, $\hat{\mu}_{\pm}$ that can couple state $\ket{m_{j}}$ to
state $\ket{m_{j}\pm 1}$, where $m_{j}\leavevmode\nobreak\
=\leavevmode\nobreak\ m_{l}\leavevmode\nobreak\ +\leavevmode\nobreak\ m_{s}$.
However, GS1 and GS2 transform as either $E_{1/2}$ (linear combination of
$\Psi_{1}$ and $\Psi_{4}$) or $E_{3/2}$ (linear combination of $\Psi_{2}$ and
$\Psi_{3}$) with the $m_{j}$ values given in Table 1. Consequently,
$\hat{\mu}_{\pm}$ cannot couple neither $\Psi_{1}$ and $\Psi_{4}$, nor
$\Psi_{2}$ and $\Psi_{3}$ therefore $g_{\perp}=0$ (cf. Table 3) in each case.
Here we note that second order contributions to $g_{\perp}$ might occur by
mixing $\Psi_{2}$ or $\Psi_{3}$ with $\Psi_{5}$, since
$\bra{\Psi_{5}}\hat{\mu}_{\pm}\ket{\Psi_{5}}\neq 0$. This may occur
vibronically or by the hyperfine interaction. Experimental ground state HF
parameters are available in the literature which are around 15-70 MHz for
Mo${}_{\text{Si}}^{+}$ in 6H SiC Baur  _et al._ (1997); 160-230 MHz for
V${}_{\text{Si}}(h)$ and 100-190 MHz for V${}_{\text{Si}}(k)$ Wolfowicz  _et
al._ (2020) in 4H SiC. HF will mix the corresponding wavefunctions only in the
second order, thus it is expected that the final $g_{\perp}$ factor will be at
least two orders of magnitude smaller than that of $g_{\parallel}$.

Calculated SO splitting and $g_{\parallel}$ values agree well with the
experimental ones giving the possibility of defect identification.
Accordingly, we identify the $\alpha$ and $\beta$ centers Wolfowicz  _et al._
(2020); Spindlberger  _et al._ (2019) as V${}_{\text{Si}}(k)$ and
V${}_{\text{Si}}(h)$ in 4H SiC, respectively, supporting the earlier
considerations reported in Ref. Spindlberger  _et al._ (2019). Regarding 6H
SiC, we identify the $\alpha$ center as V${}_{\text{Si}}(k_{2})$ while SO
splittings and values for $g_{\parallel}$ of V${}_{\text{Si}}(h)$ and
V${}_{\text{Si}}(k_{1})$ are not well-separated Wolfowicz  _et al._ (2020) for
immediate identification from the calculated $g$-constants and further
considerations are needed. To this end, we calculated the corresponding zero-
phonon lines (ZPL) that resulted in larger energy for $k_{1}$ than that for
$h$ configuration Sup . As a result we associate the $\alpha$ signal with
V${}_{\text{Si}}(k/k_{2})$ and the $\beta$ line with
V${}_{\text{Si}}(h/k_{1})$ in 4H/6H SiC and the $\gamma$ center is identified
as V${}_{\text{Si}}(h)$ in 6H SiC.

Table 3: Experimental and calculated values of $\mathbf{g}$-tensor elements ($g_{\parallel},g_{\perp}$) for GS1 and GS2 of all possible V${}_{\text{Si}}$ configurations in 4H and 6H SiC. Polyt. | Site | Experiment | Theory
---|---|---|---
GS1 | GS2 | GS1 | GS2
(PL center) |  $g_{\parallel}$,$g_{\perp}$ |  $g_{\parallel}$,$g_{\perp}$ |  $g_{\parallel}$,$g_{\perp}$ |  $g_{\parallel}$,$g_{\perp}$
4H | $h(\beta)$ | 1.870,<1111Ref. Wolfowicz  _et al._ (2020) | 2.035,<1111Ref. Wolfowicz  _et al._ (2020) | 1.975,0 | 1.987,0
$k(\alpha)$ | 1.748,0222Ref. Kaufmann  _et al._ (1997); Baur  _et al._ (1997) | 2.160,0222Ref. Kaufmann  _et al._ (1997); Baur  _et al._ (1997) | 1.866,0 | 2.106,0
6H | $h(\gamma)$ | 1.933,<1111Ref. Wolfowicz  _et al._ (2020) | 1.972,<1111Ref. Wolfowicz  _et al._ (2020) | 1.989,0 | 1.983,0
$k_{1}(\beta)$ | 1.95,-111Ref. Wolfowicz  _et al._ (2020) | 2.00,-111Ref. Wolfowicz  _et al._ (2020) | 1.987,0 | 1.985,0
$k_{2}(\alpha)$ | 1.749,0222Ref. Kaufmann  _et al._ (1997); Baur  _et al._ (1997) | -,- | 1.860,0 | 2.108,0

Our results have implications on quantum information processing based on solid
state defect qubits. It has been proposed Bosma  _et al._ (2018); Spindlberger
_et al._ (2019) that Mo and V dopants in SiC with optical transitions near or
inside telecom wavelength bands make possible to integrate these solid state
qubits to telecommunication technology. Indeed, all-optical identification and
coherent control of ensemble Mo center have been realized Bosma  _et al._
(2018). Parallel to our study, vanadium defects have been isolated and
coherent control of single spins have been demonstrated with showing all the
ingredients required for a highly efficient spin-photon interface Wolfowicz
_et al._ (2020). Our study shows the nature of the ground state spin of these
systems, namely, the order of spin levels and the origin of zero-field-
splitting. This knowledge is crucial in optimizing the quantum optics
protocols. Our results reveal the microscopic mechanism behind the phenomena
of giant anisotropy in the interaction of the electron spin with the external
stray magnetic fields which is only observable in the parallel component but
minor in the transverse components.

In summary, we carried out hybrid-DFT calculations in order to reveal
microscopic origin behind the highly anisotropic magnetic properties of KD
systems as observed in experiments Bosma  _et al._ (2018). To this end, we
shed light on the ground state electronic structure and we calculated the
corresponding SO splittings. We found that V${}_{\text{Si}}(k)$ and
V${}_{\text{Si}}(k_{2})$ exhibit one order of magnitude larger SO splittings
than that of V${}_{\text{Si}}(h)$ in 4H and V${}_{\text{Si}}(h)$,
V${}_{\text{Si}}(k_{1})$ in 6H SiC, respectively. This allowed us to identify
the $\alpha$ centers in 4H and 6H SiC as V${}_{\text{Si}}(k)$ and
V${}_{\text{Si}}(k_{2})$, respectively, while the $\beta$ center in 4H SiC
corresponds to V${}_{\text{Si}}(h)$. From the spin Hamiltonian we derived
$g_{\perp}$ and $g_{\parallel}$ (Eqs. 3 and 4) and found good agreement
between the known experimental and calculated values for V${}_{\text{Si}}$
supporting the significant role of electron-phonon coupling and character of
the wavefunction in evolving the interaction of the electron spin with the
magnetic field.

The support from ÚNKP-19-3 New National Excellence Program of the Ministry of
Human Capacities of Hungary is acknowledged by A. Cs. A. G. acknowledges the
Hungarian NKFIH grants No. KKP129866 of the National Excellence Program of
Quantum-coherent materials project, No. 127902 of the EU QuantERA Nanospin
project, EU H2020 Quantum Technology Flagship project ASTERIQS (Grant No.
820394), and the National Quantum Technology Program (Grant No.
2017-1.2.1-NKP-2017-00001).

## References

  * Kaufmann  _et al._ (1997) B. Kaufmann, A. Dörnen, and F. S. Ham, Phys. Rev. B 55, 13009 (1997).
  * Kunzer  _et al._ (1993) M. Kunzer, H. D. Müller, and U. Kaufmann, Phys. Rev. B 48, 10846 (1993).
  * Spindlberger  _et al._ (2019) L. Spindlberger, A. Csóré, G. Thiering, S. Putz, R. Karhu, J. Hassan, N. Son, T. Fromherz, A. Gali, and M. Trupke, Phys. Rev. Applied 12, 014015 (2019).
  * Csóré  _et al._ (2016) A. Csóré, A. Gällström, E. Janzén, and Á. Gali, in  _Silicon Carbide and Related Materials 2015_, Materials Science Forum, Vol. 858 (Trans Tech Publications, 2016) pp. 261–264.
  * Gällström  _et al._ (2009) A. Gällström, B. Magnusson, and E. Janzén, in  _Silicon Carbide and Related Materials 2008_, Materials Science Forum, Vol. 615 (Trans Tech Publications, 2009) pp. 405–408.
  * Bosma  _et al._ (2018) T. Bosma, G. J. J. Lof, C. M. Gilardoni, O. V. Zwier, F. Hendriks, B. Magnusson, A. Ellison, A. Gällström, I. G. Ivanov, N. T. Son, R. W. A. Havenith, and C. H. van der Wal, npj Quantum Information 4, 48 (2018), exported from https://app.dimensions.ai on 2019/03/20.
  * Baur  _et al._ (1997) J. Baur, M. Kunzer, and J. Schneider, physica status solidi (a) 162, 153 (1997).
  * Kaufmann  _et al._ (1995) B. Kaufmann, A. Dörnen, and F. S. Ham (Trans Tech Publications, 1995) pp. 707–712.
  * Dörnen  _et al._ (1992) A. Dörnen, Y. Latushko, W. Suttrop, G. Pensl, S. Leibenzeder, and R. Stein (Trans Tech Publications, 1992) pp. 1213–1218.
  * Maier  _et al._ (1992) K. Maier, H. D. Müller, and J. Schneider (Trans Tech Publications, 1992) pp. 1183–1194.
  * Reinke  _et al._ (1993) J. Reinke, H. Weihrich, S. Greulich-Weber, and J. M. Spaeth, Semiconductor Science and Technology 8, 1862 (1993).
  * Bersuker (2006) I. B. Bersuker,  _The Jahn-Teller effect_ (Cambridge University Press, 2006).
  * Ivády  _et al._ (2011) V. Ivády, A. Gällström, N. T. Son, E. Janzén, and A. Gali, Phys. Rev. Lett. 107, 195501 (2011).
  * (14) See Supplementary Information at [URL] for the results on molybdenum in 4H and 6H SiC, for the details of computational methodology, polaronic wavefunctions, spin-orbit coupling and the identification of vanadium defect configurations in 6H SiC.
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Steiner  _et al._ (2016) S. Steiner, S. Khmelevskyi, M. Marsmann, and G. Kresse, Phys. Rev. B 93, 224425 (2016).
  * Heyd  _et al._ (2003) J. Heyd, G. E. Scuseria, and M. Ernzerhof, The Journal of Chemical Physics 118, 8207 (2003).
  * Ivády  _et al._ (2014) V. Ivády, R. Armiento, K. Szász, E. Janzén, A. Gali, and I. A. Abrikosov, Phys. Rev. B 90, 035146 (2014).
  * Ham (1968) F. S. Ham, Phys. Rev. 166, 307 (1968).
  * Thiering and Gali (2018) G. Thiering and A. Gali, Phys. Rev. X 8, 021063 (2018).
  * Abragam and Bleaney (1970) A. Abragam and B. Bleaney,  _Electron paramagnetic resonance of transition ions_ (Oxford: Clarendon Press, 1970).
  * Chibotaru and Ungur (2012) L. F. Chibotaru and L. Ungur, The Journal of Chemical Physics 137, 064112 (2012), https://doi.org/10.1063/1.4739763 .
  * Stevens and Pryce (1953) K. W. H. Stevens and M. H. L. Pryce, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 219, 542 (1953).
  * Wolfowicz  _et al._ (2020) G. Wolfowicz, C. P. Anderson, B. Diler, O. G. Poluektov, F. J. Heremans, and D. D. Awschalom, Science Advances 6 (2020), 10.1126/sciadv.aaz1192, https://advances.sciencemag.org/content/6/18/eaaz1192.full.pdf .