# Changeover phenomenon in randomly colored Potts models

Nir Schreiber  Department of Mathematics, Bar Ilan University, Ramat Gan,
Israel 5290002  Reuven Cohen  Department of Mathematics, Bar Ilan University,
Ramat Gan, Israel 5290002  Gideon Amir  Department of Mathematics, Bar Ilan
University, Ramat Gan, Israel 5290002  Simi Haber  Department of Mathematics,
Bar Ilan University, Ramat Gan, Israel 5290002

###### Abstract

A hybrid Potts model where a random concentration $p$ of the spins assume
$q_{0}$ states and a random concentration $1-p$ of the spins assume $q>q_{0}$
states is introduced. It is known that when the system is homogeneous, with an
integer spin number $q_{0}$ or $q$, it undergoes a second or a first order
transition, respectively. It is argued that there is a concentration
$p^{\ast}$ such that the transition nature of the model is changed at
$p^{\ast}$. This idea is demonstrated analytically and by simulations for two
different types of interaction: the usual square lattice nearest neighboring
and mean field all-to-all. Exact expressions for the second order critical
line in concentration-temperature parameter space of the mean field model
together with some other related critical properties, are derived.

##  I Introduction

The Potts model Potts (1952); Wu (1982) has been extensively studied in the
past few decades and the amount of related literature is enormous. Though the
model has some experimental realizations, it is mostly useful as a laboratory
for analytical concepts and simulations methods Baxter (1973); Nienhuis  _et
al._ (1979); Nauenberg and Scalapino (1980); Cardy  _et al._ (1980); Mukamel
_et al._ (1976); Wilson (1989); Honmura  _et al._ (1984); Feldmann  _et al._
(1998); Graner and Glazier (1992); Bayong  _et al._ (1999); Binder (1981);
Swendsen and Wang (1987); Ferrenberg and Swendsen (1988); Wang and Landau
(2001a, b). In the context of temperature-driven phase transitions, the Potts
model can be realized as a generalization to the celebrated Ising model with
the typical “up-down” symmetry, as it has a multi-fold symmetry that is
spontaneously broken.

An interesting problem that has been first studied by Baxter Baxter (1973) is
the dependence of the transition nature (first versus second order) on the
spin number $q$. Baxter considered the square lattice ferromagnetic model with
nearest neighbors interaction (henceforth referred to as the standard model)
and obtained an exact expression for the latent heat. That expression remained
finite for $q>4$ and vanished for $q=4$. Using an equivalence of the Potts
model to a six-vertex model Temperley and Lieb (1971) Baxter has shown that
the free energy of the latter near the critical temperature was negligibly
small. Based on these findings, Baxter conjectured that the model in subject
exhibits a continuous transition for $q\leq q_{c}$ and discontinuous
transition for $q>q_{c}$ where $q_{c}=4$ is the changeover integer or the
maximal integer for which the transition is continuous. Recently, Duminil-
Copin et al Duminil-Copin  _et al._ (2016a, b) have rigorously proven Baxter’s
conjecture using the random cluster or the Fortuin-Kasteleyn Fortuin and
Kasteleyn (1972) representation of the Potts model. Another interesting
manifestation of $q_{c}=4$ being a changeover integer is the presence of a
multiplicative logarithmic correction term to the leading order power law
divergence at criticality of, e.g., the specific heat and the magnetic
susceptibility Cardy  _et al._ (1980); Salas and Sokal (1997). This is in
contrast to $q<4$ where this term is absent.

Though the changeover integer $q_{c}=4$ has been found for the square lattice,
it is expected, due to universality, to hold for other translation-invariant
lattices Nienhuis  _et al._ (1979). More precisely, two dimensional systems
accompanying the spontaneous breaking of the $q$-fold Potts symmetry, are
expected to maintain the changeover behavior of the standard model for a wide
range of local interaction patterns. There are, however, a few counterexamples
Schreiber  _et al._ (2019); Tamura  _et al._ (2010); Tanaka  _et al._ (2011).

In Tamura  _et al._ (2010); Tanaka  _et al._ (2011) the authors introduced the
so-called Potts model with invisible colors (IC) where the spins assume $q$
“visible” colors (Potts states) as in the local interaction standard model and
additional $r$ “invisible” colors that control the entropy of but do not
affect the energy of the system. The model has a growing interest in recent
years and problems such as its marginal dimensions Krasnytska  _et al._ (2016)
or its behavior on general Bethe lattices Ananikian  _et al._ (2013), thin
graphs Johnston and Ranasinghe (2013) or scale-free networks Sarkanych and
Krasnytska (2019), have been addressed.

It has been demonstrated Tamura  _et al._ (2010) using a mean field solution
to the Blume-Emery-Griffiths model Blume  _et al._ (1971) together with
simulations, and also rigorously proven van Enter  _et al._ (2011, 2012), that
the IC model exhibits a first order transition for any $q$ when $r$ is
sufficiently large. A non-trivial outcome of the latter result is the
occurrence of a first order transition for $q\leq 4$. The reason is that,
typically, in the continuous transition cases large fractal-like
configurations buffer between order and disorder, and, even though the number
of such fractals grows exponentially with their size, it does not compensate
the entropic cost of their construction when the total number of microstates
is large enough.

We present a hybrid Potts (HP) model that captures the concept of manipulating
the transition order, while keeping the interaction parameters fixed, by
introducing inhomogeneity in the number of colors. Strictly speaking, we
consider a Potts Hamiltonian ${\cal H}(\\{\sigma\\})$ 111Any Potts Hamiltonian
known to describe a system undergoing a ferromagnetic phase transition may
apply. where $\\{\sigma\\}$ is a configuration of $N$ spins, each can select
one out of the colors $\\{1,2,...,Q_{i}\\}$ where

$Q_{i}=-X_{i}(q-q_{0})+q,\ i=1,...,N\;,$ (1)

are i.i.d. random variables with $X_{i}\sim{\rm{Ber}}(p)$. In other words,
each spin can be colored in $q_{0}\leq q_{c}$ “strong” colors with probability
$p$ and in $q>q_{c}$ colors (containing the same $q_{0}$ colors and additional
$q-q_{0}$ “weak” colors), with probability $1-p$. We say that spins chosen
with probability $p$ or $1-p$ belong to strong or weak regions, respectively.

Since the HP model assumes varying concentrations of spins with different
types of colors, between limiting homogeneous instances that lead to a
transition which is either continuous or discontinuous, one might expect that
at some concentration, $0<p^{\ast}<1$, the model displays a changeover
phenomenon, that is, switches from one type of transition to another for every
combination of $q_{0}\leq q_{c}$ strong colors and a total number of $q>q_{c}$
colors.

The rest of the paper is organized as follows. In section II we find
$p^{\ast}$ for the standard model. We also present the results of Monte-Carlo
(MC) simulations. A detailed analysis of the mean field (MF) model equipped
with the HP prescription (1) is reported in section III. Finally, our
observations are summarized and discussed in section IV.

##  II The standard model

In the present section we detect $p^{\ast}$ using a simple energy-entropy
first principle argument. To be more specific, we compute the entropy loss
against the energy gain in generating long range order (LRO), where it is
known that in a second order transition, LRO is usually established by
monochromatically coloring fractal clusters that form an exponential family
with positive geometric entropy, and in a first order transition LRO is
associated with the formation of monochromatic clusters that are taken from a
sub-exponential family having zero geometric entropy. By comparing the changes
in the corresponding free energies, $p^{\ast}$ can be found.

We start with a system of $N$ spins situated on the vertices of the square
lattice and interacting via the usual Hamiltonian

${\cal H}_{B}=-J\sum_{\langle i,j\rangle}\delta_{\sigma_{i},\sigma_{j}}\;,$
(2)

where $J>0$ is the ferromagnetic coupling constant ($J$, together with
$k_{B}$, is henceforth set to one for convenience) and the summation is taken
over nearest neigboring spins.

We call a large cluster simple if its number of sites per bond is minimal,
i.e., $1/2+o(1)$. We term a snake, a fractal with a maximal number of $1+o(1)$
sites per bond. Any fractal therefore grows no faster than a snake and has
$1/2+\delta+o(1)$ sites per bond, where $\delta$ is a positive parameter no
larger than $1/2$.

Consider a fractal contour made of $l=O(N)$ bonds. The contour has energy
$-l$. Suppose the fractal has $n$ sites where $k$ of them are positioned in
strong regions and $n/l=1/2+\delta+o(1)$. Let $\\#(k,n)$ be the number of such
fractals. The expected number of these fractals is

$\langle\\#(k,n)\rangle=\mu^{n}{n\choose k}p^{k}(1-p)^{n-k}\;,$ (3)

where $\sum_{k}\\#(k,n)\sim\mu^{n}$ and the growth number $\mu$ in general
depends on $\delta$ Schreiber  _et al._ (2019). In making such fractals
monochromatic, the entropy of the system is changed in the amount of
$\ln\left(\\#(k,n)q_{0}^{-k}q^{-(n-k)}\right)$. Assuming the quantity
$\\#(k,n)$ is narrowly distributed around the mean (NDAM) and using (3), the
change in the free energy per site with inverse temperature $\beta$, can be
written (to leading order in $N$)

$\displaystyle-\beta\Delta f_{\rm{frac}}=\frac{2\beta}{1+2\delta}+\ln\mu-\ln
q$ $\displaystyle+\sum_{k}\Bigg{(}\frac{1}{n}\ln\left({n\choose
k}p^{k}(1-p)^{n-k}\right)+\frac{k}{n}\ln\left(\frac{q}{q_{0}}\right)\Bigg{)}\;.$

Under the assumption of NDAM, the sum in (II) can be replaced with the maximum
of the summand obtained at $k$ satisfying

$\kappa=\frac{pq}{pq+(1-p)q_{0}}\;,$ (5)

with $\kappa=k/n$. This brings (II) into the form

$\displaystyle-\beta\Delta f_{\rm{frac}}$ $\displaystyle=$
$\displaystyle\frac{2\beta}{1+2\delta}+\ln\mu-\ln q-\kappa\ln\kappa$ (6)
$\displaystyle-$ $\displaystyle(1-\kappa)\ln(1-\kappa)+\kappa\ln p$
$\displaystyle+$
$\displaystyle(1-\kappa)\ln(1-p)+\kappa\ln\left(\frac{q}{q_{0}}\right)\;.$

Since, unlike fractals, large simple clusters grow sub-exponentially with
their size, they are randomly distributed across the lattice. Accordingly, the
change in the free energy per site due to simple clusters with the same energy
$-l$ is given by

$\displaystyle-\beta\Delta f_{\rm{sim}}=2\beta-p\ln q_{0}-(1-p)\ln{q}\;.$ (7)

It may be more constructive to form large monochromatic simple clusters rather
than fractals. Specifically, within the framework of (6),(7), if at $\beta$
solving $\Delta f_{\rm{sim}}=0$ we have $\Delta f_{\rm{frac}}\geq 0$, then it
is entropically disadvantageous for the system to possess large fractals at
that temperature. Instead, large simple monochromatic clusters are formed and
the system undergoes a first order transition Schreiber  _et al._ (2019) at

$\beta_{c}\approx\frac{1}{2}\left(p\ln q_{0}+(1-p)\ln{q}\right)\;.$ (8)

The marginal concentration $p^{\ast}$ below which large simple clusters are
entropically more favorable than fractals can be estimated by taking
$\kappa^{\ast}$ to satisfy (5) at $p^{\ast}$, plugging it together with the
RHS of (8) into (6) and solving

$\displaystyle\sup_{\delta}$ $\displaystyle\Bigg{(}$
$\displaystyle\frac{p^{\ast}\ln q_{0}+(1-p^{\ast})\ln q}{1+2\delta}+\ln\mu-\ln
q$ (9) $\displaystyle-$
$\displaystyle\kappa^{\ast}\ln\kappa^{\ast}-(1-\kappa^{\ast})\ln(1-\kappa^{\ast})+\kappa^{\ast}\ln
p^{\ast}$ $\displaystyle+$
$\displaystyle(1-\kappa^{\ast})\ln(1-p^{\ast})+\kappa^{\ast}\ln\left(\frac{q}{q_{0}}\right)\Bigg{)}=0$

for $p^{\ast}$. It may be useful to write (9) as

$\sup_{\delta}\varphi(\delta,p^{\ast},q_{0},q)=0\;.$ (10)

Apart from uncovering the dependence of $p^{\ast}$ on the spin numbers
$q_{0},q$, (10) tells that $p^{\ast}$ is indeed marginal, since, at a
temperature that is candidate to be the first order critical point, $\Delta
f_{\rm{frac}}=0$, even for fractals that are taken from the most probable
exponential family.

For a large total number of colors, a first order transition occurs when the
concentration of strong regions is small. The reason is that, for $q$ large,
the construction of monochromatic exponential families of fractals is
entropically too expensive when the system is too sparsely populated by
$q_{0}$-type spins. The critical behavior in this case is determined by the
simple clusters and the transition is of first order. In order to quantify how
sparse the strong regions are, it may be useful to find $p^{\ast}$ in the
large $q$ limit. To this end we expand $\varphi(\delta,p^{\ast},q_{0},q)$ up
to $o(p^{\ast})$ and then solve (10) for $p^{\ast}$, keeping the leading order
large $q$ term. This results in a monotonically decreasing behavior of
$p^{\ast}$. In particular, if for varying $q$ the sup in (10) is determined by
values of $\delta$ that are not arbitrarily small, i.e., satisfying
$\inf_{q}\delta>0$, we have

$p^{\ast}\sim\frac{\ln q}{q}\;,$ (11)

which means that a first order transition takes place for small
concentrations, typically below $\frac{\ln q}{q}$, of strong regions.

It should be noted that even when the concentration of strong spins is well
below the percolation threshold, the presence of these spins can affect the
transition nature. Specifically, when $q$ is large, sparsely populated strong
spins may lead to a second order transition since they considerably reduce the
entropic cost of forming a monochromatic cluster having one of the strong
colors.

Our next goal is to numerically support the theory presented up to this point.
For that purpose, we performed extensive MC simulations using the Wang-Landau
(WL) method Wang and Landau (2001a, b) to calculate the density of states with
energy $E$, $\Omega(E)$. Being equipped with the quantities $\Omega(E)$
enables us to compute the energy probability-distribution function (PDF) at
any desired temperature. In particular, we are interested in the PDF at
pseudo-criticality, e.g., at a temperature where the specific heat
$C_{L}=L^{-d}\beta^{2}(\langle E^{2}\rangle-\langle E\rangle^{2})$ with the
energy moments given by $\langle
E^{n}\rangle=\frac{\sum_{E}E^{n}\Omega(E)e^{-\beta
E}}{\sum_{E}\Omega(E)e^{-\beta E}}$, is maximal Schreiber  _et al._ (2018,
2019). The pseudo-critical (PC) temperatures, $T_{L}$, together with the
maximal values of the specific heat for different values of $p$, for a sample
with a linear size $L=40$ ($L=\sqrt{N}$) and spin numbers $q_{0}=2,\ q=50$,
are given in Table 1. In order to roughly estimate the error bars of $T_{L}$
we ran the WL simulations $10$ times for a concentration $p=1$ and $L=40$.
Taking the std of the PC temperatures obtained from these runs, where the
temperature resolution in the canonical ensemble is 0.0001, gives $0.0008$,
thus, a plausible estimate of the error bars would be $0.001$. The error bars
of the maxima are computed using the jackknife method.

Table 1:  PC temperatures, $T_{L}$, and maximal values of the specific heat, $\max C_{L}$, for different concentrations. The spin numbers are $q_{0}=2,\ q=50$ and the lattice linear size is $L=40$. The uniform error bar ($0.001$) of the temperatures is roughly estimated using $10$ WL runs. The jackknife method is employed to compute the error bars of $\max C_{L}$. $p$ | $T_{L}$ | $\max C_{L}$
---|---|---
0 | 0.4793 | 3989$\pm$ 16
0.05 | 0.4977 | 3405$\pm$ 9
0.1 | 0.5132 | 173.29$\pm$ 0.01
0.15 | 0.5316 | 61.498$\pm$ 0.001

In Fig. 1, a plot of the specific heat against temperature is given for
different concentrations. It can be verified that the peaks for the
concentrations $p=0$ and $p=0.05$ are similar in magnitude and width. On the
other hand, the peaks of the larger concentrations are (at least) an order of
magnitude smaller and also broader. This may qualitatively indicate different
scaling behaviors and thus different types of transitions.

Figure 1: Specific heat as a function of temperature, for a sample with a
linear size $L=40$, and spin numbers $q_{0}=2,\ q=50$, displayed on a
semilogarithmic scale. The peaks are positioned at temperatures given in Table
1. Note the differences in the peaks’ magnitude between the two lowest
concentrations and the other concentrations.

In Fig. 2 we plot the PC energy PDF against the energy (per site). The
temperatures at which the PDFs are computed are those presented in Table 1.
For $p=0$, the values of the ordered and disordered energies (the positions of
the two peaks) are $\varepsilon_{o}\simeq-1.9275$ and
$\varepsilon_{d}\simeq-0.3544$, respectively, yielding
$-(\varepsilon_{o}+\varepsilon_{d})\simeq 2.2819$, in excelent agreement with
$-(\varepsilon_{o}+\varepsilon_{d})=2(1+1/\sqrt{50})=2.282842...$ Janke
(1993). The excellent agreement between the PC temperature and the exact
$T_{c}=1/\ln(1+\sqrt{50})=0.478862...$ Hintermann  _et al._ (1978) should be
also noted. Furthermore, since $L=40$ is larger than the correlation length
Buddenoir and Wallon (1993), it is likely that the typical first order scaling
law $|T_{L}-T_{c}|=O(L^{-d})$ ($|T_{L}-T_{c}|\simeq 0.0004$, c.f.
$40^{-2}\simeq 0.0006$), is satisfied.

In a second order transition, the PC energy PDF is expected to display a
single peak centred around the PC energy. On the other hand, in a first order
transition, a double peaked distribution centred around the energies of the
coexisting ordered and disordered states is conventional Janke (1993).
Evidently, while the broad PDF for $p=0.1$ is a manifestation of either a weak
first order transition or a single second order peak that suffers from finite
size effects, two pronounced peaks are observed for concentrations lower than
$p=0.1$ and a clear single peak is present for $p=0.15$. The marginal
concentration for a first order transition is therefore expected to be
$0.05<p^{\ast}\leq 0.15$, in reasonable agreement with (11) yielding
$p^{\ast}\sim\frac{\ln 50}{50}\simeq 0.08$.

Figure 2: PC energy PDF for the bond-interaction model on the square lattice.
The “strong” Potts variable is $q_{0}=2$ and the total number of colors is
$q=50$. The lattice linear size is $L=40$.

##  III Mean field model

In the present section we adopt the Bragg-Williams approach to solve the MF HP
model. The MF Hamiltonian is given by

${\cal H}_{MF}=-\frac{1}{N}\sum_{i<j}\delta_{\sigma_{i},\sigma_{j}}\;,$ (12)

where the normalization factor $N^{-1}$ assures that the total energy is
extensive. Let $\xi_{i}$ and $\eta_{i}$ be the fraction of spins with color
$i\in\\{0,1,...,q_{0}-1,...,q-1\\}$, in strong and weak regions, respectively.
The total energy can be expressed as ${\cal
E}=-\frac{1}{2}N\sum_{i}(p\xi_{i}+(1-p)\eta_{i})^{2}$ and the number of states
with this energy is given by
$W={pN\choose{pN\xi_{0},...,pN\xi_{q-1}}}{(1-p)N\choose{(1-p)N\eta_{0},...,(1-p)N\eta_{q-1}}}$.
Thus, introducing the Lagrange multipliers $a,b$, the free energy per site
takes the form

$\displaystyle\beta f=\lim_{N\to\infty}\frac{1}{N}(\beta{\cal E}-\ln
W)+{\mathrm{constraints}}$ (13) $\displaystyle=$
$\displaystyle\sum_{i}\Big{(}p\xi_{i}\ln\xi_{i}+(1-p)\eta_{i}\ln\eta_{i}-\frac{1}{2}\beta(p\xi_{i}+(1-p)\eta_{i})^{2}\Big{)}$
$\displaystyle+$ $\displaystyle
a\Big{(}\sum_{i}\xi_{i}-1\Big{)}+b\Big{(}\sum_{i}\eta_{i}-1\Big{)}\;.$

Setting $\xi_{j}=0,\ j=q_{0},...q-1$, differentiating (13) with respect to
$\xi_{j},\eta_{j}$ and equating to zero give

$\displaystyle\ln\xi_{j}-\beta(p\xi_{j}+(1-p)\eta_{j})$
$\displaystyle+a/p+1=0\;$ $\displaystyle,j=0,...,q_{0}-1\;,$
$\displaystyle\ln\eta_{j}-\beta(p\xi_{j}+(1-p)\eta_{j})$
$\displaystyle+b/(1-p)+1=0\;$ $\displaystyle,j=0,...,q-1\;.$ (14)

The fractions $\xi_{j},\eta_{j}$ therefore satisfy the coupled equations

$\xi_{j}=\begin{cases}e^{\beta(p\xi_{j}+(1-p)\eta_{j})}/Z_{0}&,j=0,...,q_{0}-1\\\
0&,j=q_{0},...,q-1\;,\\\ \end{cases}$ (15)

and

$\eta_{j}=\begin{cases}e^{\beta(p\xi_{j}+(1-p)\eta_{j})}/Z&,j=0,...,q_{0}-1\\\
e^{\beta(1-p)\eta_{j}}/Z&,j=q_{0},...,q-1\;,\\\ \end{cases}$ (16)

where
$Z_{0}=e^{\frac{a}{p}+1}=\sum_{i=0}^{q_{0}-1}e^{\beta(p\xi_{i}+(1-p)\eta_{i})}$
and $Z=e^{\frac{b}{1-p}+1}=Z_{0}+\sum_{i=q_{0}}^{q-1}e^{\beta(1-p)\eta_{i}}$
are the colors partition functions in strong and weak regions, respectively.
The quantities in (15),(16) can take the form

$\xi_{j}=\begin{cases}\frac{1}{q_{0}}\left(1+(q_{0}-1)m_{0}\right)&,j=0\\\
\frac{1}{q_{0}}(1-m_{0})&,j=1,...,q_{0}-1\\\ 0&,j=q_{0},...,q-1\;,\\\
\end{cases}$ (17)

and

$\eta_{j}=\begin{cases}\frac{1}{q}(1+(q_{0}-1)m_{0})\\\
\times(1+(q/q_{0}-1)m_{1})&,j=0\\\
\frac{1}{q}(1-m_{0})(1+(q/q_{0}-1)m_{1})&,j=1,...,q_{0}-1\\\
\frac{1}{q}\left(1-m_{1}\right)&,j=q_{0},...,q-1\;,\end{cases}$ (18)

where $m_{0},m_{1}$ are the components of a two-fold magnetization
$\boldsymbol{m}$. Indeed, there are favored fractions $\xi_{0}\geq\xi_{j}$ and
$\eta_{0}\geq\eta_{j}$ for all $1\leq j\leq q-1$, assuring the existence of
LRO. Eqs. (17),(18) also make sure that the ratio between the favored $j=0$
color and the other $q_{0}-1$ strong colors in both regions is preserved.

Above the critical point, at the disordered state, strong colors are uniformly
distributed so that (17),(18) imply

$m_{0}=0\;.$ (19)

In the case that the transition is continuous, when plugging (17),(18) into
(13), the magnetization at the critical point, $\boldsymbol{m}^{\ast}$, must
satisfy $\nabla_{\boldsymbol{m}}f\equiv(g_{0},g_{1})=0$ where (19) holds. It
is known Wu (1982); Kihara  _et al._ (1954) that for the homogeneous model
$q_{c}=2$. Thus, in the following we set $q_{0}=2$. The gradient components
then read

$\displaystyle g_{0}$ $\displaystyle=$
$\displaystyle-\frac{1}{2q^{2}}\Bigg{(}(1-p)q\left(2+m_{1}(q-2)\right)\ln\left(\frac{1-m_{0}}{1+m_{0}}\right)+\beta
m_{0}(m_{1}(1-p)(q-2)+p(q-2)+2)^{2}\Bigg{)}+p\tanh^{-1}(m_{0})\;,$
$\displaystyle g_{1}$ $\displaystyle=$
$\displaystyle-\frac{(1-p)(q-2)}{2q^{2}}\Bigg{(}\beta\left(1+m_{0}^{2}\right)q(m_{1}(1-p)+p)+2\beta
m_{0}^{2}(1-m_{1})(1-p)+q\left(m_{0}\ln\left(\frac{1-m_{0}}{1+m_{0}}\right)\right.$
(21) $\displaystyle-$
$\displaystyle\left.2\ln\left(\frac{2+m_{1}(q-2)}{2q}\right)-\ln\left(1-m_{0}^{2}\right)\right)\Bigg{)}\;.$

Substituting (19) in (III), $g_{0}$ indeed vanishes. The critical temperature
is obtained by the further condition that the magnetization at criticality is
an unstable point of the free energy. This can be formulated by considering
the Hessian matrix $H(\boldsymbol{m},\beta,p,q)$ given by

$\left(\begin{array}[]{cc}-\frac{(m_{1}(1-p)(q-2)-p(q-2)-2)\left(\beta(1-m_{0}^{2})(qp-2p+m_{1}(1-p)(q-2)+2)-2q\right)}{2\left(1-m_{0}^{2}\right)q^{2}}&\frac{(1-p)(q-2)\left(2\beta
m_{0}(-qp+2p-m_{1}(1-p)(q-2)-2)-q\ln\left(\frac{1-m_{0}}{1+m_{0}}\right)\right)}{2q^{2}}\\\
\frac{(1-p)(q-2)\left(2\beta
m_{0}(-qp+2p-m_{1}(1-p)(q-2)-2)-q\ln\left(\frac{1-m_{0}}{1+m_{0}}\right)\right)}{2q^{2}}&-\frac{(1-p)(q-2)\left(\beta(1-p)((q-2)m_{0}^{2}+q)-\frac{2q^{2}}{(1-m_{1})(m_{1}(q-2)+2)}\right)}{2q^{2}}\\\
\end{array}\right)\;,$ (22)

computed at $\boldsymbol{m}^{\ast}$ satisfying $m_{0}^{\ast}=0,\
g_{1}(\boldsymbol{m}^{\ast},\beta_{c},p,q)=0$, where the Hessian has a
vanishing eigenvalue, or, equivalently, obeys

$\det H(\boldsymbol{m}^{\ast}(\beta_{c},p,q),\beta_{c},p,q)=0\;.$ (23)

Fixing $q$, (23) implicitly determines the second order critical line in
concentration-temperature plane by

$\frac{1}{2}\exp{\left(\frac{\beta_{c}-q}{q-2}\right)}=\frac{\beta_{c}-2}{(q-2)(2-\beta_{c}p)}\;.$
(24)

The critical magnetization reads

$\boldsymbol{m}^{\ast}=\left(0,\frac{q(2-\beta_{c}p)}{\beta_{c}(1-p)(q-2)}-\frac{2}{q-2}\right)\;.$
(25)

Plugging (25) into the energy part of (13) yields the following expression for
the critical energy

$\varepsilon_{c}=-\frac{(\beta_{c}-4)\beta_{c}+2q}{2\beta_{c}^{2}(q-2)}\;.$
(26)

It is expected that the critical temperature, $\beta_{c}$, behaves as a
continuous function of the concentration $p\in[0,1]$ 222Indeed, the implicit
functions theorem guarantees that, for $q$ fixed, (24) defines a continuously
differentiable function $\beta_{c}=\psi(p),\ p\in(0,1)$.. However, the
solution to (24) at $p=0$ gives $\beta_{c}=q$, whereas the homogeneous system
undergoes a first order transition at $\beta_{c}=\frac{2(q-1)\ln(q-1)}{q-2}$
Wu (1982). Thus, continuity implies that there exists a finite concentration
$p^{\ast}$ where for $p<p^{\ast}$ (24) does not hold, and the transition
becomes a first order transition. It should be noted that $p^{\ast}$ is
expected to vary with $q$.

A first order transition is associated with a discontinuity of the
magnetization at the critical point. In other words, there are two points
$\boldsymbol{m}^{\ast}$ and $\boldsymbol{m}^{\ast\ast}$ that simultaneously
minimize the free energy, that is, simultaneously solve
$\nabla_{\boldsymbol{m}}f=0,\
f(\boldsymbol{m}^{\ast})=f(\boldsymbol{m}^{\ast\ast})$.

The nature of the transition and, consequently, the presence of $p^{\ast}$,
can be phenomenologically detected by fixing the concentration and examining
the behavior of $\boldsymbol{m}$ that numerically minimizes (13) when the
temperature is varied. In Fig. 3, the “strong” component $m_{0}$ is plotted
against temperature, for $q=6$ and different concentrations. Apparently, for
$p\leq 0.2$ the magnetization has a discontinuity in the vicinity of $T_{0}$
(taken to be the point where the numerical derivative $\frac{\Delta
m_{0}}{\Delta T}$ is “large”), indicating the presence of
$m_{0}^{\ast}=0,m_{0}^{\ast\ast}>0$ 333In practice we measure, e.g., for
$p=0.1$, $m_{0}(0.2560)=0.9149$ associated with $m_{0}^{\ast\ast}$ and
$m_{0}(0.2561)=2.337\times 10^{-5}$ associated with $m_{0}^{\ast}$, such that
$f(0.2560)=-2.0711\simeq f(0.2561)=-2.0242$, whereas for $p\geq 0.4$ it is
continuous at $T_{0}$. It is therefore expected that $0.2<p^{\ast}\leq 0.4$
for this spin number.

Interestingly, as can be verified from the inset of Fig. 3, for $q=1000$ there
is a concentration ($p=0.3$) such that, when cooled, the system first exhibits
a continuous transition and then a discontinuous transition at a lower
temperature. A similar behavior has been found for other concentrations. For
$q=6$ this phenomenon has not been clearly observed. Thus, it is expected that
for small $q$, the dual transition is either a weak effect or does not exist.
For instance, the $p=0.2$ graph may obscure a small finite magnetization in
the vicinity of a second order critical point. If such a finite magnetization
exists, then $p^{\ast}<0.2$. This may exemplify the observation that in the
case where in addition to the second order line, a first order branch is
launched, $p^{\ast}$ being the concentration where the second order line
terminates, is well defined.

We next discuss the results of Metropolis Metropolis  _et al._ (1953) MC
simulations performed to capture numerically the concentration-dependent
changeover phenomenon for the MF model. We consider the simulated
magnetization (SM)

${\rm{SM}}=p\frac{q_{0}z_{0}-1}{q_{0}-1}+(1-p)\frac{qz-1}{q-1}\;,$ (27)

where $z_{0}$ and $z$ are the maximal fractions of monochromatic spins in
strong and weak locations, respectively, and $q_{0}=2$. A sample of $N=500$
spins is used in the simulations where the SM together with the energy (per
site) proportional to (12) are observed over a total number of $5\times
10^{5}$ MC sweeps. In Fig. 4(a) it is seen that the SM performs multiple
“flip-flops” between two values of ordered and disordered states, as expected
from a system undergoing a first order transition. On the other hand, in Fig.
4(b) the system exhibits a typical second order transition by means of
fluctuations around a single value.

Figure 3: Strong magnetization component $m_{0}$ as a function of temperature,
for different concentrations. The total number of colors is $q=6$. It should
be noted that the symbols composing the graphs for $p\geq 0.4$ are more dense
than for $p\leq 0.2$, in some neighborhood of $T_{0}$. A combination of
continuous and discontinuous transitions at critical temperatures near
$T_{0}=0.1518$ and $T_{0}=0.1021$, respectively, for $q=1000$ and $p=0.3$, is
presented in the inset.  Figure 4: MC time variation of the simulated MF
magnetization given by (27). The simulated system has $N=500$ spins and a
total spin number $q=6$. For clarity, the presented observables capture a
portion of $10\%$ of the total number of sweeps. Two different concentrations
corresponding to a first or a second order transition together with
temperatures in the vicinity of $T_{0}$ (to be precise
$T=T_{0}+\frac{5}{N}=T_{0}+0.01$), are chosen. (a) $p=0.1$ and $T=0.2661$. (b)
$p=0.6$ and $T=0.4199$. The magnetization is averaged over $10$ different
runs. Energy per site (also averaged over $10$ different runs) against MC time
is plotted in the inset. The energy fluctuates around
$\langle\mathrm{energy}\rangle\simeq-0.1822$ ($\mathrm{std(energy)}\simeq
0.0037$), in plausible agreement with $\varepsilon_{c}=-0.172011...$ due to
(26).

We conclude this section by delving into the continuous regime of Fig. 3,
first by solving (24) taking $q=6$ and $p=0.6$, to give
$T_{c}=\beta_{c}^{-1}=0.409809...$. Next, from Fig. 5, where a plot of
$\boldsymbol{m}$ for $q=6$ and $p=0.6$ against temperature is presented, it
can be deduced that since the value of $T_{0}$ agrees well with $T_{c}$,
$\boldsymbol{m}$ is continuous at $T_{c}$. It should be emphasized at this
point that, unlike $m_{0}$, $m_{1}$ is actually not a proper order parameter
as it does not vanish at any finite temperature. It has, however, a
computational value.

Figure 5: Magnetization components minimizing (13) as a function of
temperature, for $q=6$ and $p=0.6$. The pseudo-critical temperature is
$T_{0}=0.4099$. The pseudo-critical magnetization at $T_{0}$ is $m_{1}=0.3235$
(c.f. $m_{1}^{\ast}=0.323567...$ due to (25)).

For instance, it can be verified from Fig. 5 that $m_{1}$ is not
differentiable at $T_{c}$ which makes $T_{c}$ a unique temperature.
Furthermore, substituting $q=6,\ p=0.6,\ \beta_{c}=2.440161...$ in (25), gives
$m_{1}^{\ast}=0.323567...$ which is in excellent agreement with $m_{1}$
computed at $T_{0}$.

##  IV Summary and discussion

It is demonstrated on the standard Potts model that in a mixed system of spins
randomly colored in $q_{0}$ or $q>q_{c}\geq q_{0}$ colors, the transition
order depends on the concentration $p$ of $q_{0}$-type spins. For that system
there is a critical concentration $0<p^{\ast}<1$ such that the transition is
discontinuous when $p$ is below $p^{\ast}$. For a small concentration of
strong regions and a large total number of colors, a changeover phenomenon is
verified analytically and numerically. In the large $q$ limit, $p^{\ast}$
displays a monotonically decreasing beheavior. This means that in the
asymptotic regime, with increasing values of $q$, it is entropically
beneficial to construct sparser monochromatic fractals that may have a
declining growth constant, than simple large-scale clusters.

The Bragg-Williams approximation, applied to the MF model, uncovers the
presence of the marginal concentration for that model, where the strong
regions are occupied by Ising-like spins. Unlike in the standard model where
it is argued that a first order transition occurs for concentrations below
$p^{\ast}$, in the MF model we conclude the presence of $p^{\ast}$ from the
other direction, that is, the second order critical temperature, as a function
of $p$, is defined on some interval with a left endpoint which is equal to
$p^{\ast}$.

It is inspected that for $q$ large there are concentrations hosting both
continuous and discontinuous transitions. A similar phenomenon, namely, a
sequential continuous and discontinuous transition occuring at some values of
$r$, has been reported for the IC model Krasnytska  _et al._ (2016) where the
Bragg-Williams approximation in the context of the MF model was employed.

A reasonable explanation for the large $q$ dual transition is that, at the
higher temperature it is entropically favorable to typically leave the weak
regions disordered, hence, Ising-like long range order is established in the
strong regions. However, when the temperature is lowered, at some point it
becomes energetically beneficial to expand the predetermined LRO by invading
the weak regions. It is yet not clear whether the dual transition is
manifested for any $q>2$.

The simulations produce observables that have a different qualitative behavior
for different concentrations. In the standard model case the PC energy PDF
changes its shape from dual to single peaked, with increasing values of $p$
and in the MF case the magnetization has a (noisy) typical single versus non-
single valued temporal behavior, for different values of $p$. Thus, in both
cases, the essence of $p^{\ast}$ as a changeover point is revealed.

Our model provides a rather general framework of governing the transition
order of well known Potts systems, without interfering with the interaction
content of such systems when no heterogeneity in the number of colors is
introduced, in the sense that a changeover phenomenon is observed for any
setup of strong and weak colors. We believe that the HP machinery can be
applied to other Potts systems such that a concentration-dependent changeover
phenomenon is detected in those systems.

## References

  * Potts (1952) R. B. Potts, Mathematical Proceedings of the Cambridge Philosophical Society 48, 106 (1952).
  * Wu (1982) F. Y. Wu, Reviews of Modern Physics 54, 235 (1982).
  * Baxter (1973) R. J. Baxter, Journal of Physics C: Solid State Physics 6, L445 (1973).
  * Nienhuis  _et al._ (1979) B. Nienhuis, A. N. Berker, E. K. Riedel, and M. Schick, Physical Review Letters 43, 737 (1979).
  * Nauenberg and Scalapino (1980) M. Nauenberg and D. J. Scalapino, Physical Review Letters 44, 837 (1980).
  * Cardy  _et al._ (1980) J. L. Cardy, M. Nauenberg, and D. J. Scalapino, Physical Review B 22, 2560 (1980).
  * Mukamel  _et al._ (1976) D. Mukamel, M. E. Fisher, and E. Domany, Physical Review Letters 37, 565 (1976).
  * Wilson (1989) W. G. Wilson, Physics Letters A 137, 398 (1989).
  * Honmura  _et al._ (1984) R. Honmura, E. F. Sarmento, C. Tsallis, and I. P. Fittipaldi, Physical Review B 29, 2761 (1984).
  * Feldmann  _et al._ (1998) H. Feldmann, A. J. Guttmann, I. Jensen, R. Shrock, and S.-H. Tsai, Journal of Physics A: Mathematical and General 31, 2287 (1998).
  * Graner and Glazier (1992) F. Graner and J. A. Glazier, Physical Review Letters 69, 2013 (1992).
  * Bayong  _et al._ (1999) E. Bayong, H. T. Diep, and V. Dotsenko, Physical Review Letters 83, 14 (1999).
  * Binder (1981) K. Binder, Journal of Statistical Physics 24, 69 (1981).
  * Swendsen and Wang (1987) R. H. Swendsen and J.-S. Wang, Physical Review Letters 58, 86 (1987).
  * Ferrenberg and Swendsen (1988) A. M. Ferrenberg and R. H. Swendsen, Physical Review Letters 61, 2635 (1988).
  * Wang and Landau (2001a) F. Wang and D. P. Landau, Physical Review Letters 86, 2050 (2001a).
  * Wang and Landau (2001b) F. Wang and D. P. Landau, Physical Review E 64 (2001b), 10.1103/physreve.64.056101.
  * Temperley and Lieb (1971) H. N. V. Temperley and E. H. Lieb, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 322, 251 (1971).
  * Duminil-Copin  _et al._ (2016a) H. Duminil-Copin, V. Sidoravicius, and V. Tassion, Communications in Mathematical Physics 349, 47 (2016a).
  * Duminil-Copin  _et al._ (2016b) H. Duminil-Copin, M. Gagnebin, M. Harel, I. Manolescu, and V. Tassion, “Discontinuity of the phase transition for the planar random-cluster and potts models with $q>4$,” (2016b), arXiv:1611.09877 .
  * Fortuin and Kasteleyn (1972) C. Fortuin and P. Kasteleyn, Physica 57, 536 (1972).
  * Salas and Sokal (1997) J. Salas and A. D. Sokal, Journal of Statistical Physics 88, 567 (1997).
  * Schreiber  _et al._ (2019) N. Schreiber, R. Cohen, S. Haber, G. Amir, and B. Barzel, Physical Review E 100 (2019).
  * Tamura  _et al._ (2010) R. Tamura, S. Tanaka, and N. Kawashima, Progress of Theoretical Physics 124, 381 (2010).
  * Tanaka  _et al._ (2011) S. Tanaka, R. Tamura, and N. Kawashima, Journal of Physics: Conference Series 297, 012022 (2011).
  * Krasnytska  _et al._ (2016) M. Krasnytska, P. Sarkanych, B. Berche, Y. Holovatch, and R. Kenna, Journal of Physics A: Mathematical and Theoretical 49, 255001 (2016).
  * Ananikian  _et al._ (2013) N. Ananikian, N. S. Izmailyan, D. A. Johnston, R. Kenna, and R. P. K. C. M. Ranasinghe, Journal of Physics A: Mathematical and Theoretical 46, 385002 (2013).
  * Johnston and Ranasinghe (2013) D. A. Johnston and R. P. K. C. M. Ranasinghe, Journal of Physics A: Mathematical and Theoretical 46, 225001 (2013).
  * Sarkanych and Krasnytska (2019) P. Sarkanych and M. Krasnytska, Physics Letters A 383, 125844 (2019).
  * Blume  _et al._ (1971) M. Blume, V. J. Emery, and R. B. Griffiths, Physical Review A 4, 1071 (1971).
  * van Enter  _et al._ (2011) A. C. D. van Enter, G. Iacobelli, and S. Taati, Progress of Theoretical Physics 126, 983 (2011).
  * van Enter  _et al._ (2012) A. C. D. van Enter, G. Iacobelli, and S. Taati, Reviews in Mathematical Physics 24, 1250004 (2012).
  * Note (1) Any Potts Hamiltonian known to describe a system undergoing a ferromagnetic phase transition may apply.
  * Schreiber  _et al._ (2018) N. Schreiber, R. Cohen, and S. Haber, Physical Review E 97 (2018).
  * Janke (1993) W. Janke, Physical Review B 47, 14757 (1993).
  * Hintermann  _et al._ (1978) A. Hintermann, H. Kunz, and F. Y. Wu, Journal of Statistical Physics 19, 623 (1978).
  * Buddenoir and Wallon (1993) E. Buddenoir and S. Wallon, Journal of Physics A: Mathematical and General 26, 3045 (1993).
  * Kihara  _et al._ (1954) T. Kihara, Y. Midzuno, and T. Shizume, Journal of the Physical Society of Japan 9, 681 (1954).
  * Note (2) Indeed, the implicit functions theorem guarantees that, for $q$ fixed, (24) defines a continuously differentiable function $\beta_{c}=\psi(p),\ p\in(0,1)$.
  * Note (3) In practice we measure, e.g., for $p=0.1$, $m_{0}(0.2560)=0.9149$ associated with $m_{0}^{\ast\ast}$ and $m_{0}(0.2561)=2.337\times 10^{-5}$ associated with $m_{0}^{\ast}$, such that $f(0.2560)=-2.0711\simeq f(0.2561)=-2.0242$.
  * Metropolis  _et al._ (1953) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, The Journal of Chemical Physics 21, 1087 (1953).


Now at ]Department of Chemistry and Chemical Technology, Bronx Community
College, Bronx, New York 10453

# First measurement of high-energy reactor antineutrinos at Daya Bay

F. P. An  Sun Yat-Sen (Zhongshan) University, Guangzhou  W. D. Bai  Sun Yat-
Sen (Zhongshan) University, Guangzhou  A. B. Balantekin  University of
Wisconsin, Madison, Wisconsin 53706  M. Bishai  Brookhaven National
Laboratory, Upton, New York 11973  S. Blyth  Department of Physics, National
Taiwan University, Taipei  G. F. Cao  Institute of High Energy Physics,
Beijing  J. Cao  Institute of High Energy Physics, Beijing  J. F. Chang
Institute of High Energy Physics, Beijing  Y. Chang  National United
University, Miao-Li  H. S. Chen  Institute of High Energy Physics, Beijing  H.
Y. Chen  Department of Engineering Physics, Tsinghua University, Beijing  S.
M. Chen  Department of Engineering Physics, Tsinghua University, Beijing  Y.
Chen  Shenzhen University, Shenzhen  Sun Yat-Sen (Zhongshan) University,
Guangzhou  Y. X. Chen  North China Electric Power University, Beijing  J.
Cheng  North China Electric Power University, Beijing  Z. K. Cheng  Sun Yat-
Sen (Zhongshan) University, Guangzhou  J. J. Cherwinka  University of
Wisconsin, Madison, Wisconsin 53706  M. C. Chu  Chinese University of Hong
Kong, Hong Kong  J. P. Cummings  Siena College, Loudonville, New York 12211
O. Dalager  Department of Physics and Astronomy, University of California,
Irvine, California 92697  F. S. Deng  University of Science and Technology of
China, Hefei  Y. Y. Ding  Institute of High Energy Physics, Beijing  M. V.
Diwan  Brookhaven National Laboratory, Upton, New York 11973  T. Dohnal
Charles University, Faculty of Mathematics and Physics, Prague  D. Dolzhikov
Joint Institute for Nuclear Research, Dubna, Moscow Region  J. Dove
Department of Physics, University of Illinois at Urbana-Champaign, Urbana,
Illinois 61801  D. A. Dwyer  Lawrence Berkeley National Laboratory, Berkeley,
California 94720  J. P. Gallo  Department of Physics, Illinois Institute of
Technology, Chicago, Illinois 60616  M. Gonchar  Joint Institute for Nuclear
Research, Dubna, Moscow Region  G. H. Gong  Department of Engineering Physics,
Tsinghua University, Beijing  H. Gong  Department of Engineering Physics,
Tsinghua University, Beijing  W. Q. Gu  Brookhaven National Laboratory, Upton,
New York 11973  J. Y. Guo  Sun Yat-Sen (Zhongshan) University, Guangzhou  L.
Guo  Department of Engineering Physics, Tsinghua University, Beijing  X. H.
Guo  Beijing Normal University, Beijing  Y. H. Guo  Department of Nuclear
Science and Technology, School of Energy and Power Engineering, Xi’an Jiaotong
University, Xi’an  Z. Guo  Department of Engineering Physics, Tsinghua
University, Beijing  R. W. Hackenburg  Brookhaven National Laboratory, Upton,
New York 11973  S. Hans  [  Brookhaven National Laboratory, Upton, New York
11973  M. He  Institute of High Energy Physics, Beijing  K. M. Heeger  Wright
Laboratory and Department of Physics, Yale University, New Haven, Connecticut
06520  Y. K. Heng  Institute of High Energy Physics, Beijing  Y. K. Hor  Sun
Yat-Sen (Zhongshan) University, Guangzhou  Y. B. Hsiung  Department of
Physics, National Taiwan University, Taipei  B. Z. Hu  Department of Physics,
National Taiwan University, Taipei  J. R. Hu  Institute of High Energy
Physics, Beijing  T. Hu  Institute of High Energy Physics, Beijing  Z. J. Hu
Sun Yat-Sen (Zhongshan) University, Guangzhou  H. X. Huang  China Institute of
Atomic Energy, Beijing  J. H. Huang  Institute of High Energy Physics, Beijing
X. T. Huang  Shandong University, Jinan  Y. B. Huang  Guangxi University,
No.100 Daxue East Road, Nanning  P. Huber  Center for Neutrino Physics,
Virginia Tech, Blacksburg, Virginia 24061  D. E. Jaffe  Brookhaven National
Laboratory, Upton, New York 11973  K. L. Jen  Institute of Physics, National
Chiao-Tung University, Hsinchu  X. L. Ji  Institute of High Energy Physics,
Beijing  X. P. Ji  Brookhaven National Laboratory, Upton, New York 11973  R.
A. Johnson  Department of Physics, University of Cincinnati, Cincinnati, Ohio
45221  D. Jones  Department of Physics, College of Science and Technology,
Temple University, Philadelphia, Pennsylvania 19122  L. Kang  Dongguan
University of Technology, Dongguan  S. H. Kettell  Brookhaven National
Laboratory, Upton, New York 11973  S. Kohn  Department of Physics, University
of California, Berkeley, California 94720  M. Kramer  Lawrence Berkeley
National Laboratory, Berkeley, California 94720  Department of Physics,
University of California, Berkeley, California 94720  T. J. Langford  Wright
Laboratory and Department of Physics, Yale University, New Haven, Connecticut
06520  J. Lee  Lawrence Berkeley National Laboratory, Berkeley, California
94720  J. H. C. Lee  Department of Physics, The University of Hong Kong,
Pokfulam, Hong Kong  R. T. Lei  Dongguan University of Technology, Dongguan
R. Leitner  Charles University, Faculty of Mathematics and Physics, Prague  J.
K. C. Leung  Department of Physics, The University of Hong Kong, Pokfulam,
Hong Kong  F. Li  Institute of High Energy Physics, Beijing  H. L. Li
Institute of High Energy Physics, Beijing  J. J. Li  Department of Engineering
Physics, Tsinghua University, Beijing  Q. J. Li  Institute of High Energy
Physics, Beijing  R. H. Li  Institute of High Energy Physics, Beijing  S. Li
Dongguan University of Technology, Dongguan  S. C. Li  Center for Neutrino
Physics, Virginia Tech, Blacksburg, Virginia 24061  W. D. Li  Institute of
High Energy Physics, Beijing  X. N. Li  Institute of High Energy Physics,
Beijing  X. Q. Li  School of Physics, Nankai University, Tianjin  Y. F. Li
Institute of High Energy Physics, Beijing  Z. B. Li  Sun Yat-Sen (Zhongshan)
University, Guangzhou  H. Liang  University of Science and Technology of
China, Hefei  C. J. Lin  Lawrence Berkeley National Laboratory, Berkeley,
California 94720  G. L. Lin  Institute of Physics, National Chiao-Tung
University, Hsinchu  S. Lin  Dongguan University of Technology, Dongguan  J.
J. Ling  Sun Yat-Sen (Zhongshan) University, Guangzhou  J. M. Link  Center for
Neutrino Physics, Virginia Tech, Blacksburg, Virginia 24061  L. Littenberg
Brookhaven National Laboratory, Upton, New York 11973  B. R. Littlejohn
Department of Physics, Illinois Institute of Technology, Chicago, Illinois
60616  J. C. Liu  Institute of High Energy Physics, Beijing  J. L. Liu
Department of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai
Laboratory for Particle Physics and Cosmology, Shanghai  J. X. Liu  Institute
of High Energy Physics, Beijing  C. Lu  Joseph Henry Laboratories, Princeton
University, Princeton, New Jersey 08544  H. Q. Lu  Institute of High Energy
Physics, Beijing  K. B. Luk  Department of Physics, University of California,
Berkeley, California 94720  Lawrence Berkeley National Laboratory, Berkeley,
California 94720  B. Z. Ma  Shandong University, Jinan  X. B. Ma  North China
Electric Power University, Beijing  X. Y. Ma  Institute of High Energy
Physics, Beijing  Y. Q. Ma  Institute of High Energy Physics, Beijing  R. C.
Mandujano  Department of Physics and Astronomy, University of California,
Irvine, California 92697  C. Marshall  Lawrence Berkeley National Laboratory,
Berkeley, California 94720  K. T. McDonald  Joseph Henry Laboratories,
Princeton University, Princeton, New Jersey 08544  R. D. McKeown  California
Institute of Technology, Pasadena, California 91125  College of William and
Mary, Williamsburg, Virginia 23187  Y. Meng  Department of Physics and
Astronomy, Shanghai Jiao Tong University, Shanghai Laboratory for Particle
Physics and Cosmology, Shanghai  J. Napolitano  Department of Physics, College
of Science and Technology, Temple University, Philadelphia, Pennsylvania 19122
D. Naumov  Joint Institute for Nuclear Research, Dubna, Moscow Region  E.
Naumova  Joint Institute for Nuclear Research, Dubna, Moscow Region  T. M. T.
Nguyen  Institute of Physics, National Chiao-Tung University, Hsinchu  J. P.
Ochoa-Ricoux  Department of Physics and Astronomy, University of California,
Irvine, California 92697  A. Olshevskiy  Joint Institute for Nuclear Research,
Dubna, Moscow Region  H.-R. Pan  Department of Physics, National Taiwan
University, Taipei  J. Park  Center for Neutrino Physics, Virginia Tech,
Blacksburg, Virginia 24061  S. Patton  Lawrence Berkeley National Laboratory,
Berkeley, California 94720  J. C. Peng  Department of Physics, University of
Illinois at Urbana-Champaign, Urbana, Illinois 61801  C. S. J. Pun  Department
of Physics, The University of Hong Kong, Pokfulam, Hong Kong  F. Z. Qi
Institute of High Energy Physics, Beijing  M. Qi  Nanjing University, Nanjing
X. Qian  Brookhaven National Laboratory, Upton, New York 11973  N. Raper  Sun
Yat-Sen (Zhongshan) University, Guangzhou  J. Ren  China Institute of Atomic
Energy, Beijing  C. Morales Reveco  Department of Physics and Astronomy,
University of California, Irvine, California 92697  R. Rosero  Brookhaven
National Laboratory, Upton, New York 11973  B. Roskovec  Charles University,
Faculty of Mathematics and Physics, Prague  X. C. Ruan  China Institute of
Atomic Energy, Beijing  H. Steiner  Department of Physics, University of
California, Berkeley, California 94720  Lawrence Berkeley National Laboratory,
Berkeley, California 94720  J. L. Sun  China General Nuclear Power Group,
Shenzhen  T. Tmej  Charles University, Faculty of Mathematics and Physics,
Prague  K. Treskov  Joint Institute for Nuclear Research, Dubna, Moscow Region
W.-H. Tse  Chinese University of Hong Kong, Hong Kong  C. E. Tull  Lawrence
Berkeley National Laboratory, Berkeley, California 94720  B. Viren  Brookhaven
National Laboratory, Upton, New York 11973  V. Vorobel  Charles University,
Faculty of Mathematics and Physics, Prague  C. H. Wang  National United
University, Miao-Li  J. Wang  Sun Yat-Sen (Zhongshan) University, Guangzhou
M. Wang  Shandong University, Jinan  N. Y. Wang  Beijing Normal University,
Beijing  R. G. Wang  Institute of High Energy Physics, Beijing  W. Wang  Sun
Yat-Sen (Zhongshan) University, Guangzhou  College of William and Mary,
Williamsburg, Virginia 23187  X. Wang  College of Electronic Science and
Engineering, National University of Defense Technology, Changsha  Y. Wang
Nanjing University, Nanjing  Y. F. Wang  Institute of High Energy Physics,
Beijing  Z. Wang  Institute of High Energy Physics, Beijing  Z. Wang
Department of Engineering Physics, Tsinghua University, Beijing  Z. M. Wang
Institute of High Energy Physics, Beijing  H. Y. Wei  Brookhaven National
Laboratory, Upton, New York 11973  L. H. Wei  Institute of High Energy
Physics, Beijing  L. J. Wen  Institute of High Energy Physics, Beijing  K.
Whisnant  Iowa State University, Ames, Iowa 50011  C. G. White  Department of
Physics, Illinois Institute of Technology, Chicago, Illinois 60616  H. L. H.
Wong  Department of Physics, University of California, Berkeley, California
94720  Lawrence Berkeley National Laboratory, Berkeley, California 94720  E.
Worcester  Brookhaven National Laboratory, Upton, New York 11973  D. R. Wu
Institute of High Energy Physics, Beijing  Q. Wu  Shandong University, Jinan
W. J. Wu  Institute of High Energy Physics, Beijing  D. M. Xia  Chongqing
University, Chongqing  Z. Q. Xie  Institute of High Energy Physics, Beijing
Z. Z. Xing  Institute of High Energy Physics, Beijing  H. K. Xu  Institute of
High Energy Physics, Beijing  J. L. Xu  Institute of High Energy Physics,
Beijing  T. Xu  Department of Engineering Physics, Tsinghua University,
Beijing  T. Xue  Department of Engineering Physics, Tsinghua University,
Beijing  C. G. Yang  Institute of High Energy Physics, Beijing  L. Yang
Dongguan University of Technology, Dongguan  Y. Z. Yang  Department of
Engineering Physics, Tsinghua University, Beijing  H. F. Yao  Institute of
High Energy Physics, Beijing  M. Ye  Institute of High Energy Physics, Beijing
M. Yeh  Brookhaven National Laboratory, Upton, New York 11973  B. L. Young
Iowa State University, Ames, Iowa 50011  H. Z. Yu  Sun Yat-Sen (Zhongshan)
University, Guangzhou  Z. Y. Yu  Institute of High Energy Physics, Beijing  B.
B. Yue  Sun Yat-Sen (Zhongshan) University, Guangzhou  V. Zavadskyi  Joint
Institute for Nuclear Research, Dubna, Moscow Region  S. Zeng  Institute of
High Energy Physics, Beijing  Y. Zeng  Sun Yat-Sen (Zhongshan) University,
Guangzhou  L. Zhan  Institute of High Energy Physics, Beijing  C. Zhang
Brookhaven National Laboratory, Upton, New York 11973  F. Y. Zhang  Department
of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai Laboratory
for Particle Physics and Cosmology, Shanghai  H. H. Zhang  Sun Yat-Sen
(Zhongshan) University, Guangzhou  J. L. Zhang  Nanjing University, Nanjing
J. W. Zhang  Institute of High Energy Physics, Beijing  Q. M. Zhang
Department of Nuclear Science and Technology, School of Energy and Power
Engineering, Xi’an Jiaotong University, Xi’an  S. Q. Zhang  Sun Yat-Sen
(Zhongshan) University, Guangzhou  X. T. Zhang  Institute of High Energy
Physics, Beijing  Y. M. Zhang  Sun Yat-Sen (Zhongshan) University, Guangzhou
Y. X. Zhang  China General Nuclear Power Group, Shenzhen  Y. Y. Zhang
Department of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai
Laboratory for Particle Physics and Cosmology, Shanghai  Z. J. Zhang  Dongguan
University of Technology, Dongguan  Z. P. Zhang  University of Science and
Technology of China, Hefei  Z. Y. Zhang  Institute of High Energy Physics,
Beijing  J. Zhao  Institute of High Energy Physics, Beijing  R. Z. Zhao
Institute of High Energy Physics, Beijing  L. Zhou  Institute of High Energy
Physics, Beijing  H. L. Zhuang  Institute of High Energy Physics, Beijing  J.
H. Zou  Institute of High Energy Physics, Beijing

###### Abstract

This Letter reports the first measurement of high-energy reactor antineutrinos
at Daya Bay, with nearly 9000 inverse beta decay candidates in the prompt
energy region of 8-12 MeV observed over 1958 days of data collection. A
multivariate analysis is used to separate 2500 signal events from background
statistically. The hypothesis of no reactor antineutrinos with neutrino energy
above 10 MeV is rejected with a significance of 6.2 standard deviations. A 29%
antineutrino flux deficit in the prompt energy region of 8-11 MeV is observed
compared to a recent model prediction. We provide the unfolded antineutrino
spectrum above 7 MeV as a data-based reference for other experiments. This
result provides the first direct observation of the production of
antineutrinos from several high-$Q_{\beta}$ isotopes in commercial reactors.

Recently, antineutrino measurements at commercial Mention  _et al._ (2011);
Adey  _et al._ (2019a); An  _et al._ (2016a, 2017a); Choi  _et al._ (2016);
Abe  _et al._ (2014) and research Ashenfelter  _et al._ (2019); Andriamirado
_et al._ (2021); Almazán  _et al._ (2020, 2021, 2022) reactors have challenged
both the flux and spectral shape of model predictions using either the
conversion method Schreckenbach  _et al._ (1985); Von Feilitzsch  _et al._
(1982); Hahn  _et al._ (1989); Mueller  _et al._ (2011); Huber (2011); Haag
_et al._ (2014) or the summation method Mueller  _et al._ (2011); Fallot  _et
al._ (2012); Dwyer and Langford (2015); Estienne  _et al._ (2019). These
discrepancies have spurred substantial advancement of the theory and
reexamination of experimental basis of the predictions , although these
improvements have not substantially altered the underlying data-model
disagreements Kopeikin  _et al._ (2021); Fallot  _et al._ (2012); Hayes  _et
al._ (2014); Estienne  _et al._ (2019); Berryman and Huber (2020); Bernstein
_et al._ (2019); Hayen  _et al._ (2019); Giunti  _et al._ (2022).

In these experiments, antineutrinos are detected via the inverse beta decay
(IBD) reaction ($\bar{\nu}_{e}+p\rightarrow e^{+}+n$) Vogel and Beacom (1999).
Energy of the antineutrino (Eν) is estimated using the energy of the detected
prompt signal (Ep) which is generated by the positron and its annihilation
gammas. An approximate relationship
${E}_{\bar{\nu}}\approx{E}_{p}+0.8~{}\text{MeV}$. All previous measurements
have focused on Ep$<$8 MeV because signals in the higher-energy region are
rare and often contaminated by cosmogenic backgrounds. Despite the rarity,
high-energy reactor antineutrinos represent a background for future
measurements, such as the diffuse supernova neutrino background that is
expected to permeate the Universe Beacom and Vagins (2004).

High-energy reactor antineutrinos are likely generated by only a handful of
short-lived $\beta$-decay nuclei with high end-point energies ($Q_{\beta}$),
such as 88,90Br and 94,96,98Rb. Daughters of these high-$Q_{\beta}$ nuclei
have complex deexcitation pathways that may have been incorrectly determined
by measurements suffering from low efficiency for associated $\gamma$ rays
(the ”pandemonium effect” Hardy  _et al._ (1977)). Recent measurements using
the total absorption spectroscopy revealed the substantial impact of
pandemonium effects in comparatively lower $Q_{\beta}$ isotopes on the
calculation of reactor antineutrino spectra Fallot  _et al._ (2012); Zakari-
Issoufou  _et al._ (2015); Valencia  _et al._ (2017); Rice  _et al._ (2017);
Guadilla  _et al._ (2019); Estienne  _et al._ (2019); Fijałkowska  _et al._
(2017), indicating the likelihood of similar issues at higher $Q_{\beta}$
Valencia  _et al._ (2017). In this sense, direct measurements of high-energy
antineutrinos can provide a valuable new perspective for nuclear data
validations relevant well beyond the bounds of neutrino physics Algora  _et
al._ (2021); Fallot  _et al._ (2019); Schmidt  _et al._ (2021).

This Letter reports the measurement of reactor antineutrinos with Ep$>$8 MeV
using a multivariate analysis at Daya Bay. This is the first detailed analysis
of reactor-produced antineutrinos in this energy region, yielding an unfolded
antineutrino energy spectrum which provides a more reliable data-based
prediction than existing theoretical methods.

The Daya Bay experiment studies antineutrinos from the Daya Bay nuclear power
complex, which hosts six commercial pressurized-water reactor cores (each with
2.9 GW maximum thermal power). Eight identically designed antineutrino
detectors (ADs) are deployed in two near experimental halls (EH1 and EH2, each
containing two ADs) and in a far experimental hall (EH3, containing four ADs).
Details about the experiment are described in Refs. An  _et al._ (2016b,
2017b). IBD interaction candidates are selected following the criteria similar
to “selection A” described in Ref. An  _et al._ (2017b). While in previous
analyses the veto of IBD candidates following cosmic muon signals with high
reconstructed energy ($E_{\mu}^{\rm rec}>2.5$ GeV) was set to 1 s to suppress
backgrounds from cosmogenic isotopes, the veto time is shortened to 1 ms to
include these background-rich samples. This change allows exploitation of
background properties to reliably disentangle signal from background. In the
latest dataset taken in 1958 calendar days, we collected about 4 million IBD
candidates at three experimental halls but only about 9000 candidates with Ep
between 8 and 12 MeV.

Figure 1 displays the IBD candidate rate $R$ as a function of the weighted
reactor power, ${P}_{\rm reactor}$, for four Ep bins. $R$ is the sum of all
detectors without background subtraction. For each detector, the ${P}_{\rm
reactor}$ is defined as

${P}_{\rm reactor}=\frac{R_{\rm RealPower}}{R_{\rm FullPower}}\times
17.4~{}{\rm GW_{th}},$ (1)

where $R_{\rm RealPower}$ is the predicted IBD rate using the reactor
information provided by the company and $R_{\rm FullPower}$ is the predicted
rate assuming that all reactors were working with full power. It is found that
the event rate drops by more than a factor of 20 for the 8-12 MeV energy
region compared to that of 6-8 MeV energy region, suggesting a much larger
background to signal ratio in the high-Ep region. By decomposing $R$ into two
parts, $R=R_{\rm IBD}{P}_{\rm reactor}+R_{\rm bkg}$, where $R_{\rm IBD}$
represents the IBD reaction rate per unit of reactor power, and $R_{\rm bkg}$
represents the background rate which is uncorrelated with ${P}_{\rm reactor}$,
a strong correlation between $R_{\rm IBD}$ and ${P}_{\rm reactor}$ is observed
for Ep in 6 to 8 MeV that is inconsistent with a background-only hypothesis at
over 30$\sigma$. By contrast, the significance of the correlation decreases to
$<$2.5$\sigma$ above 8 MeV, primarily due to the much larger background to
signal ratio. This figure illustrates significant challenges in extracting the
IBD events with Ep$>$ 8 MeV.

Figure 1:  Relationship between IBD rates without background subtract in four
Ep bins and the weighted reactor power. The vertical error bars are
statistical uncertainties and the horizontal ones are systematic
uncertainties. The relationship is fitted with a linear function. Scopes which
represent rates of IBD signals are also shown.

The main backgrounds with Ep in 8 to 12 MeV are from muon decays, cosmogenic
fast neutrons, and cosmogenic isotope decays. Other contributions, such as
atmospheric neutrinos, are negligible. A probability distribution function
(PDF) is constructed for each event:

$F(\bm{r};\bm{\Delta t},z,w)=\sum_{p}r_{p}f_{p}(\bm{\Delta
t})h_{p}(z)k_{p}(w),$ (2)

where $p$ represents the event types (IBD, cosmogenic isotopes, or fast
neutron), $r_{p}$ is the ratio of the number of type-$p$ events over the total
event number in each Ep bin, and $\bm{r}$ represents the vector of $r_{p}$.
$f(\bm{\Delta t})$ is the expected distribution of time difference $\bm{\Delta
t}$ between the IBD candidate and its preceding muon events in the same
detector (the vector of $\Delta t$ for each of eight muon categories described
below), $h(z)$ is the distribution of vertical vertex coordinates $z$ of the
prompt signal, and $k(w)$ is the event distribution of the weighted reactor
power $w$ when the event occurred. $\bm{r}$ is determined by minimizing a
$\chi^{2}$ constructed as

$\chi^{2}(\bm{r})=-2\sum\left[\log F(\bm{r};\bm{\Delta
t},z,w)\right]+g(\bm{\epsilon}),$ (3)

where $\Sigma$ denotes the sum over all events and $g(\bm{\epsilon})$
constrains the nuisance parameters describing $f_{p}(\bm{\Delta t})$,
$h_{p}(z)$ and $k_{p}(w)$ utilizing information described below.

If a cosmic-ray muon stops in the detector, its energy deposit and the
subsequent Michel electron could mimic an IBD event. Muons were tagged with an
efficiency larger than 99.7% by the two water Cherenkov detectors surrounding
ADs. However, the efficiency decreased after 2017 because of an increase of
nonfunctional photomultipliers in the Cherenkov detectors, resulting in the
appearance of muon decay background. The background was negligible for
previous neutrino oscillation analyses. In this analysis it was found that a
$z>1.7$ m cut applied on the prompt signal can reject more than 99.8% of this
background, because stopping muons with less than 12 MeV energy deposits
concentrated in the top of AD. This background is not considered in Eq. 2 due
to the negligible contribution after applying the cut.

Neutrons generated by cosmic muon interactions outside ADs could penetrate
into ADs and mimic IBD reactions. They are dubbed ”fast neutrons”. Because of
the production mechanism, the majority of fast neutrons have downward momentum
and deposit energies near top of ADs. The $z$ distribution fast neutrons for
8$<$Ep$<$12 MeV is inferred from an almost pure fast neutron sample with
12$<$Ep$<$20 MeV. The $z$ distribution of antineutrinos is uniform and
determined using the IBD-enriched sample with $2<$Ep$<8$ MeV that has less
than 0.1% fast neutron contamination. Rates of fast neutrons obtained using
the vertex information are consistent with an alternative method based on IBD-
like events coincident with muons identified using only the outer water
Cherenkov detector An  _et al._ (2017b).

The cosmogenic production of 9Li and 8He (referred as 9Li in later text) with
subsequent $\beta-n$ decay or a coincidence of two 12B $\beta$ decays can
generate a correlated pair of signals nearly identical to IBDs. Since
cosmogenic isotopes are produced by muons penetrating through the whole
detector, their $z$ vertex distributions are found to be consistent with the
uniform distribution of IBDs. The $\Delta t$ distribution is described by
$f(\Delta t)=\kappa\cdot e^{-\kappa\Delta t}$, where $\kappa=R_{\mu}$ for
muon-uncorrelated events such as IBDs and $\kappa=R_{\mu}+\frac{n}{\tau}$ for
muon-correlated events. Here, $R_{\mu}$ is the muon rate, and $\tau$ is the
life of the isotope. Details are found in Ref Wen  _et al._ (2006), Fig. 26 in
Ref. An  _et al._ (2017b), and Fig. 2 in Ref. Adey  _et al._ (2018). The
factor $n$ is 1 for 9Li and 2 for the accidental coincidence of two 12B
$\beta$ decays. The different $\Delta t$ dependencies facilitate an effective
determination of amounts of cosmogenic isotopes. In addition, production of
unstable isotopes is usually associated with energetic particle showers
attributed to muons Abe  _et al._ (2010), which often include spallation
neutrons. Thus, muons are divided into eight categories: four $E^{\rm
rec}_{\mu}$ bins, and each bin with or without accompanying neutrons.

The cosmogenic background rates can be effectively determined from the $\Delta
t$ distributions except in the near halls (EH1 and EH2) for the lowest muon-
energy bin ($E^{\rm rec}_{\mu}<1$ GeV) without accompanying neutrons. The high
muon rates in these two cases result in indistinguishable $\Delta t$
distributions for cosmogenic isotope backgrounds and muon-uncorrelated events.
Previous analyses in Daya Bay studied the relation between isotope and neutron
productions An  _et al._ (2017b); Adey  _et al._ (2018):
$\varepsilon_{hi}=Y_{hi}^{n}/Y_{hi}^{\rm all}$, where $Y_{hi}^{\rm all}$ is
the total isotope yield from muons in the $i^{th}$ $E_{\mu}^{\rm rec}$ bin in
experimental hall $h$, and $Y_{hi}^{n}$ is the yield with accompanying
neutrons. Though EH3 is situated deeper underground than EH1 or EH2,
$\varepsilon$ is similar in all three halls. Thus the fitted $Y^{n}$ for
$E^{\rm rec}_{\mu}<1$ GeV can be combined with the determination of
$\varepsilon$ in EH3 to estimate $Y^{\rm all}$ in EH1 and EH2.

An independent method was used to cross-check 9Li yields from low-energy muons
in EH1 and EH2. A full GEANT4-based Agostinelli  _et al._ (2003) simulation of
muons was performed. For each muon track through Gd-doped liquid scintillator,
energies from steps with large energy losses are summed to yield a ”shower
energy”. The shower energy distribution $\bm{S}_{hi}^{\rm sh}$ is obtained for
muons in the $i^{\rm th}$ $E_{\mu}^{\rm rec}$ range. The relationship between
the known 9Li yields from data ($Y^{\rm all}_{hi}$) and $\bm{S}^{\rm sh}_{hi}$
is built by: $Y^{\rm all}_{hi}=\sum_{j}\bm{S}^{\rm sh}_{hij}\times\bm{Q}_{j}$,
where $\bm{Q}$ represents the 9Li yield from each muon in the $j^{th}$ shower
energy bin. When extracting the $\bm{Q}$ based on $Y^{\rm all}_{hi}$ and
$\bm{S}^{\rm sh}_{hi}$, an unfolding technique is utilized to suppress the
amplified fluctuation Almazán  _et al._ (2021). Finally, with the determined
$\bm{Q}$, the 9Li yield from the lowest muon-energy bin is found to account
for 18% (13%) of the total yield in EH1 (EH2), with 50% (55%) relative
uncertainty. These predictions are consistent with the results calculated
based on $\varepsilon$ and are used as constraints in the pull term of Eq. 3.

After the minimization of Eq. 3, using the best-fit values of $\bm{r}$, the
probability of being an IBD signal is calculated for each event:

${P}_{\rm IBD}=\frac{r_{\rm IBD}f_{\rm IBD}(\bm{\Delta t})h_{\rm IBD}(z)k_{\rm
IBD}(w)}{F(\bm{r};\bm{\Delta t},z,w)}.$ (4)

The predictions of ${P}_{\rm IBD}$ for different event types are calculated
using toy Monte Carlo datasets generated from the PDFs in Eq. 2. As
illustrated in Fig. 2, the reactor antineutrinos can be statistically
separated and the differentiating capability decreases as the prompt energy
increases. The ${P}_{\rm IBD}$ distributions from data and predictions are
consistent within statistical uncertainty, which demonstrates good agreement
between data and model.

It is found that cosmogenic isotopes are the dominant background for Ep$<$9.5
MeV and fast neutrons for Ep$>$9.5 MeV. As a sanity check, removing reactor
power information from the fitter leads to a similar result because the
reactors are almost always operating at full power. Removing the $\Delta t$
information or reconstructed vertex information yields consistent results with
significantly larger uncertainties. Analysis with the data of only one EH also
shows consistent results within 1$\sigma$. The uncertainty of the IBD yields
for Ep$>$8 MeV is dominated by systematic uncertainties from the background
estimation. Uncertainty from cosmogenic isotopes takes up to $\sim$60% of the
total uncertainty for Ep in 8 to 9.5 MeV, while the uncertainty from fast
neutrons accounts for more than 70% for Ep$>$10 MeV.

Figure 2:  (a-c) Distributions of probabilities of being an IBD signal based
on the best-fit results for different event types in different Ep bins ((a):
8-9, (b): 9-10, and (c): 10-11 MeV). The filled histograms show the expected
distributions which are calculated using toy Monte Carlo datasets generated
from the PDFs of the fitting model. (d) Ratios of the distributions between
data and predictions in (a) and (b). The error bars are statistical only. The
data points for Ep in the 9-10 MeV bin are shifted slightly left for visual
clarity.

The above analysis has extracted the antineutrino yield for Ep$>$8 MeV. An
unfolding technique is needed to remove the detector response and convert the
Ep spectrum to the Eν spectrum. This requires the Ep spectrum in lower-energy
regions. The prompt energy spectrum below 8 MeV has been updated with the new
9Li background estimation in this analysis. The result is consistent with
those released in previous publications Adey  _et al._ (2019b); An  _et al._
(2021) within uncertainties.

The prompt energy spectrum per fission above 6 MeV and the comparison with
theoretical predictions are shown in Fig. 3. An extended neutrino spectrum for
Eν$>$10 MeV was kindly provided by the authors of a recent summation model
(SM2018) Estienne  _et al._ (2019) for the comparison. The measured IBD yield
is 3% larger than the SM2018 prediction for Ep in 6-8 MeV (similar to the
conclusion from Ref. An  _et al._ (2021)) but is 29% smaller than the
prediction for Ep in 8-11 MeV. Interestingly, for Ep$>$10 MeV, while a nonzero
signal rate of greater than 5$\sigma$ statistical significance is predicted by
the SM2018 model, the significance is only 1$\sigma$ for the data. The large
deficit in the observed high-energy IBD yields is consistent with the presence
of pandemonium-affected decay data for high-$Q_{\beta}$ isotopes such as
94,96,98Rb Valencia  _et al._ (2017); Dwyer and Langford (2015) in the SM2018
summation prediction, which would result in overstated feeding to low-lying
states and an overprediction of high-energy antineutrinos. In this case,
overprediction of high-energy antineutrinos would be accompanied by
underprediction of beta-delayed neutron release by these isotopes Roberts  _et
al._ (1939); Abriola  _et al._ (2011); Birch  _et al._ (2014).

The measured high-energy IBD yields were also compared to predictions derived
from the Huber-Mueller beta conversion model Huber (2011); Mueller  _et al._
(2011). A polynomial extrapolation was used to obtain the predictions for
$E_{\nu}>$8 MeV Huber (2011). The extrapolated result is larger than the
measurement by 30% or more for Ep$>$7.5 MeV, and gives worse agreement with
data than the SM2018 model. For this reason, we recommend not using this
extrapolation in future high-energy reactor antineutrino studies.

Figure 3:  (a) Measured prompt energy spectrum compared with the prediction
from the SM2018 model. The components from different isotopes (235U, 238U,
239Pu, and 241Pu) are shown for fission fractions (0.564, 0.076, 0.304, and
0.056 respectively). (b) A enlarged plot of (a) above 8 MeV with logarithmic
vertical scales. (c) Ratio of the measurement over the prediction from SM2018
and the extrapolated Huber-Mueller (HM) model. The HM model is not valid above
8 MeV as described in the text. The error bars in the data points are the
square root of the diagonal elements of the covariance matrix, including both
the statistical uncertainty and systematic uncertainty. The result above 11
MeV is not shown due to the larger than 100% uncertainty.

A data-based reactor antineutrino energy spectrum is determined with an
unfolding technique, similar to the earlier analysis in the low-energy region
An  _et al._ (2021). We constructed a fitter by removing the detector response
based on

$\chi^{2}=(\bm{P}-\bm{M})\textbf{Cov}^{-1}(\bm{P}-\bm{M})^{T}.$ (5)

Here the $\bm{M}$ and Cov are the measured prompt energy spectrum per fission
above 6 MeV and its covariance matrix, respectively. The prediction of the
prompt energy spectrum per fission ($\bm{P}$) is calculated by:
$\bm{P}=\bm{R}\bm{S}^{\mathrm{fit}}$, where $\bm{R}$ is the detector response
matrix at Daya Bay, which maps the predicted antineutrino energy spectrum
($\bm{S}^{\mathrm{fit}}$) to prompt energy. Here
$\bm{S}_{i}^{\mathrm{fit}}=\bm{S}_{i}^{\rm init}\times\bm{\eta}_{i}$, with
free parameter $\bm{\eta}_{i}$ on the initial values ($\bm{S}_{i}^{\rm init}$)
in $i^{\rm th}$ energy bin. The starting point of the antineutrino energy
spectrum in the unfolding is set at 6 MeV to ensure that resolution-induced
feedup of 6-7 MeV antineutrinos into the high-energy region is properly
accounted for.

The unfolded antineutrino energy spectrum is shown in Fig. 4. In the unfolding
process, the postfit prediction ($\bm{P}$) is the same as the measurement
($\bm{M}$), with the best-fit $\chi^{2}$ value equal to 0. Mathematically,
this method is equivalent to the matrix inversion method, but this procedure
has the advantage of enabling a variety of statistical tests. While this
method could suffer from amplified statistical fluctuations and big bin-to-bin
anticorrelation in the unfolding process An  _et al._ (2021), these problems
are mitigated by the large bin widths used in this analysis. Assuming no
reactor antineutrinos above specific energies (10, 10.5, or 11 MeV) in
$\bm{S}^{\mathrm{fit}}$, the best-fit $\chi^{2}$ values (38.3, 1.6, and 0.03,
respectively) are obtained. Therefore, the significance of our result in
rejecting the hypothesis of no reactor antineutrinos above 10 MeV is
6.2$\sigma$. Above 10.5 MeV, the ability to reject the background-only
hypothesis is marginal. The prompt energy spectrum, unfolded antineutrino
energy spectrum and their covariance matrices are included in Supplemental
Material.

In summary, the Daya Bay experiment has determined the prompt IBD energy
spectrum up to 11 MeV. A 29% difference in IBD rate in the prompt energy
region of 8-11 MeV is found compared with a recent summation model. An
antineutrino energy spectrum is then obtained from these data using an
unfolding procedure that removes the detector response effects. The hypothesis
of no reactor antineutrinos with energy above 10 MeV is rejected with a
significance of 6.2$\sigma$. For the first time, this work extends the energy
region of reactor antineutrinos above 10 MeV by direct measurement. The
combination of very high statistics and low cosmogenic backgrounds of Daya Bay
suggests that the precision of this measurement is unlikely to be surpassed in
the foreseeable future.

Figure 4:  (a) Reactor antineutrino energy spectrum above 7 MeV and its
comparison with the prediction from SM2018. The IBD cross section Vogel and
Beacom (1999) has been removed. The fission fractions of 235U, 238U, 239Pu,
and 241Pu are 0.564, 0.076, 0.304, and 0.056 respectively. The central values
have larger deviation from SM2018 compared with Fig. 3 due to the amplified
statistical fluctuations with bin-to-bin anticorrelation in the unfolding
process An  _et al._ (2021). (b) A enlarged plot of (a) above 8.5 MeV with a
logarithmic vertical scale. (c) Ratio of the unfolded antineutrino energy
spectrum over the prediction from SM2018.

We thank Magali Estienne and Muriel Fallot for their generosity to extend
their summation model in Ref. Estienne  _et al._ (2019) above 10 MeV and their
helpful discussions.

The Daya Bay experiment is supported in part by the Ministry of Science and
Technology of China, the U.S. Department of Energy, the Chinese Academy of
Sciences, the CAS Center for Excellence in Particle Physics, the National
Natural Science Foundation of China, the Guangdong provincial government, the
Shenzhen municipal government, the China General Nuclear Power Group, the
Research Grants Council of the Hong Kong Special Administrative Region of
China, the Ministry of Education in Taiwan, the U.S. National Science
Foundation, the Ministry of Education, Youth, and Sports of the Czech
Republic, the Charles University Research Centre UNCE, the Joint Institute of
Nuclear Research in Dubna, Russia, the National Commission of Scientific and
Technological Research of Chile, We acknowledge Yellow River Engineering
Consulting Co., Ltd., and China Railway 15th Bureau Group Co., Ltd., for
building the underground laboratory. We are grateful for the ongoing
cooperation from the China Guangdong Nuclear Power Group and China Light &
Power Company.

## References

  * Mention  _et al._ (2011) G. Mention, M. Fechner, T. Lasserre, T. Mueller, D. Lhuillier, M. Cribier, and A. Letourneau, Phys. Rev. D 83, 073006 (2011), arXiv:1101.2755 [hep-ex] .
  * Adey  _et al._ (2019a) D. Adey _et al._ (Daya Bay), Phys. Rev. D 100, 052004 (2019a), arXiv:1808.10836 [hep-ex] .
  * An  _et al._ (2016a) F. P. An _et al._ (Daya Bay), Phys. Rev. Lett. 116, 061801 (2016a), [Erratum: Phys.Rev.Lett. 118, 099902 (2017)], arXiv:1508.04233 [hep-ex] .
  * An  _et al._ (2017a) F. P. An _et al._ (Daya Bay), Chin. Phys. C 41, 013002 (2017a), arXiv:1607.05378 [hep-ex] .
  * Choi  _et al._ (2016) J. H. Choi _et al._ (RENO), Phys. Rev. Lett. 116, 211801 (2016), arXiv:1511.05849 [hep-ex] .
  * Abe  _et al._ (2014) Y. Abe _et al._ (Double Chooz), JHEP 10, 086 (2014), [Erratum: JHEP 02, 074 (2015)], arXiv:1406.7763 [hep-ex] .
  * Ashenfelter  _et al._ (2019) J. Ashenfelter _et al._ (PROSPECT), Phys. Rev. Lett. 122, 251801 (2019), arXiv:1812.10877 [nucl-ex] .
  * Andriamirado  _et al._ (2021) M. Andriamirado _et al._ (PROSPECT), Phys. Rev. D 103, 032001 (2021), arXiv:2006.11210 [hep-ex] .
  * Almazán  _et al._ (2020) H. Almazán _et al._ (STEREO), Phys. Rev. Lett. 125, 201801 (2020), arXiv:2004.04075 [hep-ex] .
  * Almazán  _et al._ (2021) H. Almazán _et al._ (STEREO), J. Phys. G 48, 075107 (2021), arXiv:2010.01876 [hep-ex] .
  * Almazán  _et al._ (2022) H. Almazán _et al._ (Stereo, Prospect), Phys. Rev. Lett. 128, 081802 (2022), arXiv:2107.03371 [nucl-ex] .
  * Schreckenbach  _et al._ (1985) K. Schreckenbach, G. Colvin, W. Gelletly, and F. Von Feilitzsch, Phys. Lett. B 160, 325 (1985).
  * Von Feilitzsch  _et al._ (1982) F. Von Feilitzsch, A. Hahn, and K. Schreckenbach, Phys. Lett. B 118, 162 (1982).
  * Hahn  _et al._ (1989) A. Hahn, K. Schreckenbach, G. Colvin, B. Krusche, W. Gelletly, and F. Von Feilitzsch, Phys. Lett. B 218, 365 (1989).
  * Mueller  _et al._ (2011) T. Mueller _et al._ , Phys. Rev. C 83, 054615 (2011), arXiv:1101.2663 [hep-ex] .
  * Huber (2011) P. Huber, Phys. Rev. C 84, 024617 (2011).
  * Haag  _et al._ (2014) N. Haag, A. Gütlein, M. Hofmann, L. Oberauer, W. Potzel, K. Schreckenbach, and F. M. Wagner, Phys. Rev. Lett. 112, 122501 (2014), arXiv:1312.5601 [nucl-ex] .
  * Fallot  _et al._ (2012) M. Fallot _et al._ , Phys. Rev. Lett. 109, 202504 (2012), arXiv:1208.3877 [nucl-ex] .
  * Dwyer and Langford (2015) D. A. Dwyer and T. J. Langford, Phys. Rev. Lett. 114, 012502 (2015), arXiv:1407.1281 [nucl-ex] .
  * Estienne  _et al._ (2019) M. Estienne _et al._ , Phys. Rev. Lett. 123, 022502 (2019), arXiv:1904.09358 [nucl-ex] .
  * Kopeikin  _et al._ (2021) V. Kopeikin, M. Skorokhvatov, and O. Titov, Phys. Rev. D 104, L071301 (2021), arXiv:2103.01684 [nucl-ex] .
  * Hayes  _et al._ (2014) A. Hayes, J. Friar, G. Garvey, G. Jungman, and G. Jonkmans, Phys. Rev. Lett. 112, 202501 (2014), arXiv:1309.4146 [nucl-th] .
  * Berryman and Huber (2020) J. M. Berryman and P. Huber, Phys. Rev. D 101, 015008 (2020), arXiv:1909.09267 [hep-ph] .
  * Bernstein  _et al._ (2019) L. A. Bernstein, D. A. Brown, A. J. Koning, B. T. Rearden, C. E. Romano, A. A. Sonzogni, A. S. Voyles, and W. Younes, Ann. Rev. Nucl. Part. Sci. 69, 109 (2019).
  * Hayen  _et al._ (2019) L. Hayen, J. Kostensalo, N. Severijns, and J. Suhonen, Phys. Rev. C 100, 054323 (2019), arXiv:1908.08302 [nucl-th] .
  * Giunti  _et al._ (2022) C. Giunti, Y. F. Li, C. A. Ternes, and Z. Xin, Phys. Lett. B 829, 137054 (2022), arXiv:2110.06820 [hep-ph] .
  * Vogel and Beacom (1999) P. Vogel and J. F. Beacom, Phys. Rev. D 60, 053003 (1999), arXiv:hep-ph/9903554 .
  * Beacom and Vagins (2004) J. F. Beacom and M. R. Vagins, Phys. Rev. Lett. 93, 171101 (2004), arXiv:hep-ph/0309300 .
  * Hardy  _et al._ (1977) J. C. Hardy, L. C. Carraz, B. Jonson, and P. G. Hansen, Phys. Lett. B 71, 307 (1977).
  * Zakari-Issoufou  _et al._ (2015) A. A. Zakari-Issoufou _et al._ (IGISOL), Phys. Rev. Lett. 115, 102503 (2015), arXiv:1504.05812 [nucl-ex] .
  * Valencia  _et al._ (2017) E. Valencia _et al._ , Phys. Rev. C 95, 024320 (2017), arXiv:1609.06128 [nucl-ex] .
  * Rice  _et al._ (2017) S. Rice _et al._ , Phys. Rev. C 96, 014320 (2017), arXiv:1704.01915 [nucl-ex] .
  * Guadilla  _et al._ (2019) V. Guadilla _et al._ , Phys. Rev. Lett. 122, 042502 (2019).
  * Fijałkowska  _et al._ (2017) A. Fijałkowska, M. Karny, K. P. Rykaczewski, B. C. Rasco, R. Grzywacz, C. J. Gross, M. Woli ńska Cichocka, K. C. Goetz, D. W. Stracener, W. Bielewski, R. Goans, J. H. Hamilton, J. W. Johnson, C. Jost, M. Madurga, K. Miernik, D. Miller, S. W. Padgett, S. V. Paulauskas, A. V. Ramayya, and E. F. Zganjar, Phys. Rev. Lett. 119, 052503 (2017).
  * Algora  _et al._ (2021) A. Algora, J. L. Tain, B. Rubio, M. Fallot, and W. Gelletly, Eur. Phys. J. A 57, 85 (2021), arXiv:2007.07918 [nucl-ex] .
  * Fallot  _et al._ (2019) M. Fallot, B. Littlejohn, and P. Dimitriou, “Antineutrino spectra and their applications,” IAEA Report INDC(NDS)-0786 (2019), https://www-nds.iaea.org/publications/indc/indc-nds-0786.pdf.
  * Schmidt  _et al._ (2021) K. H. Schmidt, M. Estienne, M. Fallot, S. Cormon, A. Cucoanes, T. Shiba, B. Jurado, K. Kern, and C. Schmitt, Nucl. Data Sheets 173, 54 (2021), arXiv:2012.08199 [nucl-ex] .
  * An  _et al._ (2016b) F. An _et al._ (Daya Bay), Nucl. Instrum. Meth. A 811, 133 (2016b), arXiv:1508.03943 [physics.ins-det] .
  * An  _et al._ (2017b) F. P. An _et al._ (Daya Bay), Phys. Rev. D 95, 072006 (2017b), arXiv:1610.04802 [hep-ex] .
  * Wen  _et al._ (2006) L.-j. Wen, J. Cao, K.-B. Luk, Y.-q. Ma, Y.-f. Wang, and C.-g. Yang, Nucl. Instrum. Meth. A 564, 471 (2006), arXiv:hep-ex/0604034 .
  * Adey  _et al._ (2018) D. Adey _et al._ (Daya Bay), Phys. Rev. Lett. 121, 241805 (2018), arXiv:1809.02261 [hep-ex] .
  * Abe  _et al._ (2010) S. Abe _et al._ (KamLAND), Phys. Rev. C 81, 025807 (2010), arXiv:0907.0066 [hep-ex] .
  * Agostinelli  _et al._ (2003) S. Agostinelli _et al._ (GEANT4), Nucl. Instrum. Meth. A 506, 250 (2003).
  * Adey  _et al._ (2019b) D. Adey _et al._ (Daya Bay), Phys. Rev. Lett. 123, 111801 (2019b), arXiv:1904.07812 [hep-ex] .
  * An  _et al._ (2021) F. P. An _et al._ (Daya Bay), Chin. Phys. C 45, 073001 (2021), arXiv:2102.04614 [hep-ex] .
  * Roberts  _et al._ (1939) R. B. Roberts, L. R. Hafstad, R. C. Meyer, and P. Wang, Phys. Rev. 55, 664 (1939).
  * Abriola  _et al._ (2011) D. Abriola, B. Singh, and I. Dillmann, “Beta-delayed neutron emission evaluation,” IAEA Report INDC(NDS)-0599 (2011), https://www-nds.iaea.org/publications/indc/indc-nds-0599.pdf.
  * Birch  _et al._ (2014) M. Birch, B. Singh, D. Abriola, I. Dillmann, T. D. Johnson, E. A. McCutchan, and A. A. Sonzogni, Nucl. Data Sheets 120, 66 (2014).