# Perceptual Evaluation of 360 Audiovisual Quality and Machine Learning
Predictions

1st Randy Frans Fela  SenseLab
FORCE Technology
Hørsholm, Denmark
<EMAIL_ADDRESS>2nd Nick Zacharov  SenseLab
FORCE Technology
Hørsholm, Denmark
<EMAIL_ADDRESS>3rd Søren Forchhammer  Dept. Photonics Engineering
Technical University of Denmark
Kgs. Lyngby, Denmark
<EMAIL_ADDRESS>
###### Abstract

In an earlier study, we gathered perceptual evaluations of the audio, video,
and audiovisual quality for 360 audiovisual content. This paper investigates
perceived audiovisual quality prediction based on objective quality metrics
and subjective scores of 360 video and spatial audio content. Thirteen
objective video quality metrics and three objective audio quality metrics were
evaluated for five stimuli for each coding parameter. Four regression-based
machine learning models were trained and tested here, i.e., multiple linear
regression, decision tree, random forest, and support vector machine. Each
model was constructed using a combination of audio and video quality metrics
and two cross-validation methods (k-Fold and Leave-One-Out) were investigated
and produced 312 predictive models. The results indicate that the model based
on the evaluation of VMAF and AMBIQUAL is better than other combinations of
audio-video quality metric. In this study, support vector machine provides
higher performance using k-Fold (PCC = 0.909, SROCC = 0.914, and RMSE =
0.416). These results can provide insights for the design of multimedia
quality metrics and the development of predictive models for audiovisual
omnidirectional media.

###### Index Terms:

perceptual evaluation, 360 video, spatial audio, machine learning, multimedia
quality metrics, higher order ambisonsics.

††publicationid: pubid: 978-1-6654-3288-7/21/$31.00 ©2021 IEEE

##  I Introduction

In recent years, 360 video or omnidirectional video (ODV) has become popular
and increasingly developed to be playable in a more efficient way. The
features of ODV which allows users to explore a spherical image by rotating
their head offers the possibility to pair this type of video with spatial
audio. Several platforms, such as VLC, Youtube, and Facebook, allow users to
upload 360 audiovisual content and playback through the traditional flat
displays or head mounted displays (HMD). This technology raises the question
about how the users perceive the quality of 360 audiovisual content and how to
achieve a high-level user experience in 360 audiovisual.

As a common approach, both the affective testing and predictive measures are
employed together in perceptual quality assessment of audiovisual aiming to
provide the validation process of the results obtained from predictive
metrics. In the 2D video, eight audiovisual quality models were identified
from previous reports as comprehensively summarized in [1, 2]. In particular
with ODV, the perceptual quality has been studied in [3, 4] through subjective
and objective methods. Using these models we expect to be able to produce a
statistical-based perceptual quality model by utilizing techniques such as
curve fitting [3], metadata-based approached [5], and by incorporating
viewport information for adaptive streaming application [6]. In [7], a machine
learning-based QoE model have been proposed by converting continuous scores
into a dichotomous score in order to build a logistic regression model.

However, in order to obtain the overall impression of perceptual events, the
presence of auditory stimuli is required. In terms of spherical projection
that encourages users to look around in ODV, ambisonic spatial audio is
considered a highly compatible pair for ODV. It preserves the spatial
information of audio signals, allowing users to perceive sounds coming from
specific directions. After its first development in the late ’70s, Ambisonic,
which was firstly proposed in [8] has recently gained popularity with the
progress of virtual reality (VR) technology. In spatial audio quality
evaluation, overall quality, as well as attributes relating to spatial
qualities such as localization and timbral quality, are the central area of
interest [9]. Furthermore, there is a full-reference metrics available for
ambisonic as has been proposed in [10].

Although there are several perceptual evaluation studies of 360 video or
spatial audio individually [3, 9, 11, 4, 7], studies in immersive contents
that combine 360 video and spatial audio is relatively unexplored, thus will
be the main contribution of this paper. In this domain, our earlier work has
investigated the perceptual audio, video and audiovisual quality subjectively
[12]. In this paper, machine learning-based models, i.e., multiple linear
regression, decision tree, random forest, and support vector machine (SVM),
were investigated in order to evaluate the performance of these methods in
predicting perceptual quality models. According to this objective, we address
the following questions:

  * •

Based on the objective quality measures, which objective metric contributes to
multimodal quality in terms of correlation?

  * •

In terms of combined audiovisual quality metrics, which combination could
produce the highest correlation with subjective audiovisual quality?

  * •

Among the machine learning algorithms implemented, which algorithm could
provide the best prediction?

TABLE I: Quality specification of the audiovisual contents Audio | Video
---|---
| B-Format PCM First Order
---
Ambisonic (FOA) AmbiX
ERP 4K (3840x1920)
| 48 kHz, 16 bit, 3,072 Mbps
---
(768 kbps/channel)
| 29.97 fps, 8-bit depth
---
$\sim$30 Mbps, YUV 4:2:0
Figure 1: Temporal and spatial indexes of testing sequences.

##  II Content and Methods

###  II-A Content

The content was downloaded from Jump video dataset[13] and used with
permission from the creator. Five 360 videos in equirectangular projection
format containing first-order ambisonic (FOA) audio were carefully selected
with the internal expectation that the selected materials would be able to
elicit both the audio and visual quality responses. The specification of
source materials is listed in Table I. In order to ensure that there was
adequate visual information for subjective assessment, the temporal and
spatial index (TI & SI) of the source video were calculated as described in
[14] and the result is shown in Figure 1.

The overall workflow implemented in this study is illustrated in Figure 2.
First, audio clips were extracted from test items to be processed
independently (audio and video) for different encoding parameters. Each source
video was encoded in FFmpeg using H.264 (libx264) to create the processed
video sequences (PVS) with four quantization parameters (QP: 22, 27, 32, 37)
and four resolutions (3840x1920, 2560x1280, 1920x1080 and 1280x720).
Meanwhile, a low-bitrate codec (AAC-LC) and ambisonic decoding technique were
employed to create processed audio excerpts (PAE) in three bitrates (64 kbps,
128 kbps, and 256 kbps) and three types of channel playback (5.0, 11.0, 22.0).
Each test signal (PVS and PAE) was further processed with representative
objective quality metrics for video (OQMV) or audio (OQMA) and delivered to
the display device (head-mounted display or multichannel loudspeakers) for
laboratory testing.

Figure 2: Workflow of predictive models built in this study.

###  II-B Measures

Three full-reference of OQMA i.e., perceptual evaluation of audio quality
(PEAQ) [15], ViSQOLAudio111https://github.com/google/visqol [16, 17], and
AMBIQUAL [10] were computed to estimate overall listening quality. In terms of
video quality metrics, nine full-reference OQMV were computed by following the
common test condition and testing procedure for 360 video described in [18].
The OQMV was measured in the codec, cross-format (CF), and end-to-end (EE)
stages, including the basic peak signal-to-noise ratio (PSNR) and its
variants, i.e., weighted to spherically uniform PSNR (WS-PSNR) [19], Sphere-
based PSNR with interpolation (S-PSNR-I) [20], with the nearest neighbor
(S-PSNR-NN) [21], and craster parabolic projection PSNR (CPP-PSNR) [22]. In
addition, a VMAF (video multimethod assessment fusion) [23, 24] source code
originally developed by Netflix was also computed, generating four additional
metrics including structural similarity index metrics (SSIM) [25] and multi-
scale SSIM [26], PSNR, and VMAF. The inclusion of VMAF in this study was
motivated by the report that it is feasible to measure OQMV of 360 video by
using VMAF222https://github.com/Netflix/vmaf without any adjustment [11].

Three subjective tests (audio, video, audiovisual) with single stimulus
absolute category ratings with continuous quality scale (ACR-CQS) were
conducted and participated by twenty assessors who passed a basic audiovisual
screening test [14]. The user interface (UI) reliably displayed either on the
projection screen (for listening session) or as a pop-up interface in a
virtual environment (visual and audiovisual session). SenseLabOnline 4.0 was
used to integrate the entire test setup, define the experimental design, and
precisely run the tests in a double-blind random presentation order. Please
note that only audiovisual quality will be discusssed in this paper. The
readers interested to single modality results as well as detailed experimental
description are encouraged to refer to [12].

TABLE II: Computational setup for machine learning prediction Settings | LM | DT | RF | SVM
---|---|---|---|---
returnResamp | All | All | All | All
search | random | random | random | random
method | lm | rpart | rf | svmRadialSigma
tuning |  | minsplit: 20
---
maxdepth: 30
- |  | ntree: 500
---
tuneLength: 15
-

###  II-C Implementation and Evaluation

The machine learning models were implemented in the R programming language
with the caret package [27, 28]. Note that this is a basic benchmarking study
and four machine learning algorithms were selected based on the common
practice in perceptual quality studies. For each model, there are three
variable inputs which consist of OQMA ($n$ = 3), OQMV ($n$ = 13), and each
combination gives $n$ = 39. In contrast, the variable output is only a mean
opinion score of perceived audiovisual quality (MOSAV). This form is motivated
by the basic audiovisual quality model, which was originally expressed in
linear form, as shown in (1) [2],

$MOS_{AV}=\alpha_{1}+\beta_{1}Q_{A}+\gamma_{1}Q_{V}+\zeta_{1}Q_{A}Q_{V}$ (1)

where in this case, the symbol $Q$ could be the quality obtained from the
objective metrics of audio or video, respectively, and $MOS_{AV}$ denotes the
audiovisual quality score obtained from the subjective evaluation. The
relevant settings for the machine learning model are presented in Table II.

In order to maintain the accuracy of the model prediction due to the
relatively small dataset, we split the data into 80:20 ratio for the training
set and test set. Two types of cross validation (CV) were carried out in this
study i.e., k-Fold ($n$ = 10 splits) and leave-one-out (LOOCV) in order to
investigate how the results from content-based split will differ to random
split based on k-Fold. In LOOCV, each class of content was treated as a test
set. The remaining class (N = 5$-$1) was implemented as a training set
consecutively so that the prediction accuracy was the average of the overall
results.

####  II-C1 Decision Trees

The decision tree in R worked based on Gini impurity, which measures the
proportion of the incorrect labels of the randomly selected elements from the
set according to the label distribution in the subset[29]. We used the rpart
library for the decision tree model with the tuning parameters of a minimum
split _“minsplit”_ and maximum depth _“maxdepth”_ adjusted by default to 20
and 30, respectively. The _“minsplit”_ represents the minimum number of data
points needed to attempt a split before it is enforced to build a terminal
node, whereas the _“maxdepth”_ is the maximum number of internal nodes built
between the root and the terminal nodes.

####  II-C2 Random forest

Random forest [30] is considered more effective than decision tree when
working with large datasets and can retain consistency when missing data
exists. The rf method in the caret package was used to perform the random
forest work in R. The default tuning parameter is set to _“mtry”_ = 7 and
_“ntree”_ = 500. The parameter _“mtry”_ specifies the number of randomly
sampled variables as candidates at each split, while _“ntree”_ specifies the
number of trees to grow. Here, we configured the tuneLength parameter, which
allowed the system to adjust the algorithm automatically. It indicated the
number of different values to be tested for each adjustment parameter, for
example, _“mtry”_ as a random forest. Assuming the tuneLength = 5, which means
to try five different _“mtry”_ values and find the best _“mtry”_ value based
on these five.

####  II-C3 Support Vector Machine

The support vector machine (SVM)333This technique is originally named as SVM.
Later, a term of support vector regression (SVR) is frequently used in certain
community to distinguish the application in regression problem. is a machine
learning technique which aims to construct a hyperplane in the an $n$
dimensional space, where $n$ represents the number of features, to find a
decision boundary of two or more classes. The objective of SVM is to find the
hyperplane that could maximize the distance margin to the closest vectors
called support vectors. Hence, in the case where the goal is to predict the
label of the class is called classification problem. SVM can also be used for
regression problem where the goal is to predict the appropriate hyperplane
position relative to the support vectors. By maximizing the distance margin
between support vectors to the hyperplane could approximate the actual
function represented by data points. Support vector machines use a
hypothetical space of linear functions in a higher-dimensional feature space
which are trained using an optimization theory learning algorithm, which uses
learning biases derived from statistical learning theory. In this study, we
used support vector regression with a radial basis function (RBF) for kernel
parameters as reported in a similar study showing a greater accuracy than a
linear kernel [7]. The method chosen is svmRadialSigma as available in library
kernlab[31]. The radial basis function can be expressed in (2) as follows[32],

$K(x_{1},x_{2})=\text{exp}(-{\sigma}\|x_{1}-x_{2}\|^{2})$ (2)

In the RBF kernel, the value depends on the distance from the origin or a
certain point. Finally, the distance information of the vector in the original
space can be used to determine the dot product (similarity) of $x_{1}$ and
$x_{2}$. The tuning parameters are regularization parameter $c$ and kernel
parameter $\sigma$. Parameter $c$ controls error by adjusting the margin
distance. In the case where the value of $c$ increases and $\sigma$ decreases,
the model is overfitting. In a random search, tuneLength parameter is the
total number of ($c$, $\sigma$) pairs to evaluate.

##  III Results

TABLE III: Pearson’s correlation coefficient (PCC) of all metric combinations
for 1) k-Fold (left) and 2) LOOCV (right). TABLE IV: Spearman’s rank order
correlation coefficient (SROCC) of all metric combinations for 1) k-Fold
(left) and 2) LOOCV (right). TABLE V: Root mean squared error (RMSE) of all
metric combinations for 1) k-Fold (left) and 2) LOOCV (right).

The audio-video quality prediction metrics were computed in all possible
combinations between OQMV and OQMA as independent variables and audiovisual
subjective scores (MOSAV) as a dependent variable. In this section, we aim to
answer three questions addressed earlier in Section I. The correlation
coefficients (PCC and SROCC) and prediction errors RMSE are computed as
evaluation performance in this study between the predicted MOS and actual MOS
from the test set. In this section, all description in this section refers to
Tables III, IV and V. The column number in tables is represented as OQMA as
follows: 1) PEAQ, 2) ViSQOLAudio, and 3) AMBIQUAL.

###  III-A Perceptual Metrics

In OQMA, AMBIQUAL significantly outperforms other metrics regardless of the
algorithms and the type of cross-validation used, proving its robustness
(range of PCC: 0.731$-$0.909, SROCC: 0.725$-$0.915, RMSE: 0.697$-$0.426). In
contrast, PEAQ and ViSQOL results vary depending on the algorithm and
validation method. In a linear correlation, ViSQOL produces comparable or
higher performance over PEAQ in a linear model (PCC: 0.653$-$0.816) and
support vector machine (PCC: 0.700$-$0.831) for which the highest scores were
produced in k-Fold CV. In monotonic correlation, ViSQOL yields better
performance in LOOCV for all algorithms but random forest, ranging from 0.562
(DT2) to 0.835 (SVM2).

Meanwhile, VMAF generally performs the best scores over OQMV (range of PCC:
0.729$-$0.909, SROCC: 0.633$-$0.915, RMSE: 0.752$-$0.426). SSIM metrics
perform relatively lower than PSNR metrics specifically in linear model (LM1 &
LM2). Nevertheless, in particular paired with OQMA and mainly for LOOCV,
several metrics i.e., SSIM (PCC-RF2: 0.905, RMSE-RF2: 0.428), MS-SSIM (SROCC-
RF2 & SVM2: 0.910 & 0.914) and CF-SPSNR-NN could present a slightly comparable
(CF-SPSNR-NN) or higher score (SSIM & MS-SSIM) than VMAF.

###  III-B Machine Learning Predictions

The results indicated that regardless of the cross-validation method used,
VMAF or AMBIQUAL can improve the prediction performance both in terms of
linear and monotonic correlation, thus reducing the prediction error.
Furthermore, the VMAF$-$AMBIQUAL pair can generally achieve the highest
performance of all machine learning prediction models (PCC $\geq$0.794, SROCC
$\geq$0.804, RMSE: 0.605$-$0.426.). However, it is observed that a number of
OQMV, including CF-SPSNR-NN, SSIM, and MS-SSIM paired with AMBIQUAL produces
slightly better performance (SSIM PCC: 0.905, SROCC: 0.908 & MS-SSIM PCC:
0.900, SROCC: 0.910 both in RF2) or close to VMAF$-$AMBIQUAL (MS-SSIM PCC:
0.904, SRCC: 0.914 in SVM2). The root-mean-squared error (RMSE) values are
determined as small as 0.426 and 0.428 respectively in SVM1: VMAF$-$AMBIQUAL
and RF2: SSIM$-$AMBIQUAL. According to the comparison of the machine learning
models, the results imply that overall, the support vector machine could
produce the highest performance with PCC up to 0.909 (VMAF$-$AMBIQUAL) and
SROCC up to 0.914 (MS-SSIM$-$AMBIQUAL). Although some related studies have
found that random forest-based prediction models could better be suited among
their respective models [33], it can be said that a limitation of random
forest is that it cannot be extrapolated and the prediction is resulted only
from the average of previous data observed in a training set. Therefore, in
the regression problem, the prediction range of the random forest is bound by
the highest and lowest labels in the training data. It can be problematic when
the range or data distribution vary for the training and test sets. However,
most studies reported that random forest can perform closely, the same, or
better than SVM but repeatedly perform better than the remaining models such
as multiple linear regression and decision tree.

In cross-validation method, it is shown that the type of CV could diversify
the distribution of prediction results yielded among various pairs of metrics
and machine learning approaches. It is noticeable that k-Fold could produce
wider range than LOOCV, for instance by comparing SVM1 (PCC: 0.802$-$0.909,
SROCC: 0.779$-$0.914 & RMSE: 0.593$-$0.426) and SVM2 (PCC: 0.864$-$0.905,
SROCC: 0.852$-$0.914 & RMSE: 0.520$-$0.432). A k-Fold also produces more
consistent results in terms of different pairs of audiovisual metrics and
machine learning approaches, depicted by the VMAF$-$AMBIQUAL that outperforms
in all other models whereas the SSIM-based results can be found as superior in
LOOCV. In addition, it should be noted that the cross-validation method also
depends on the data size and type of data. This is a bias-variance trade-off
when choosing a cross-validation method for building a model. In LOOCV, since
each training set contains only $n-1$ examples, the estimation of the test
error has a lower bias. It can produce a higher variance, which means that it
will use the training set for each iteration. Since there is a large overlap
between the training sets, the effect is a larger variance, which means that
the average of the test estimates of the test error will have a more
significant variance. In contrast, the overlap between the training sets in
k-Fold CV is relatively small, so the correlation of the test error estimates
is small. As a result, the average test error value will not have as large a
variance as LOOCV. At the end of this results, based on Tables III, IV, and V,
k-Fold shows better results than LOOCV.

##  IV Conclusions

Perceptual evaluation of audiovisual quality was carried out through
subjective experiment for 360 video over the head-mounted display and low-
bitrate ambisonic over loudspeaker playback. The perceptual subjective
results, as well as perceptual quality models based on subjective data, have
been described earlier in [12]. This study evaluates a combination of
objective quality metrics computed for ambisonic audio and 360 video to build
the perceptual audiovisual quality models by utilizing machine learning
prediction. A test methodology, as well as the test results, have been
described, and the main conclusions can be summarized as follows:

  * •

In general, VMAF and AMBIQUAL, respectively, provide consistent performance
better than the other video/audio quality metrics in terms of the highest
correlation coefficient hence lowering prediction error (PCC $\geq$0.794,
SROCC $\geq$0.804, RMSE: 0.605$-$0.426.).

  * •

The combination of audio-video metrics analyzed for 360 audiovisual contents
demonstrates that VMAF$-$AMBIQUAL outperforms other combinations and presents
a good agreement in any condition and tested machine learning algorithms (PCC:
0.794$-$0.909, SROCC: 0.804$-$0.914, RMSE: 0.605$-$0.426.).

  * •

Corresponding to cross-validation techniques, there is a slightly different
performance between k-Fold and leave-one-out cross-validation for all machine
learning algorithms but decision tree, where noticeable differences were
identified (e.g., SVM1: PCC 0.802$-$0.909, SROCC 0.779$-$0.914 & RMSE
0.593$-$0.426; SVM2: PCC 0.864$-$0.905, SROCC 0.852$-$0.914 & RMSE
0.520$-$0.432).

##  V Future Works

At the time of the writing of this paper, it could be stated that a
comprehensive 360 audiovisual quality database remains scarce. As expected,
our subjective assessment as described in [12] showed that the MOS score did
not fully span the scale range. Therefore, the availability of a 360
audiovisual database, which we are currently developing, is highly critical in
order to collect better quality data from retraining and building the models.

As the results shows in Tables III, IV, and V, combination of certain audio-
video quality metrics reveals a promising correlation with subjective data. It
is therefore exposing a new interest to consider a number of potential
objective metrics used in various application to be included. Although the
metric is typically application-dependent (speech, music,
annotation)444https://www.amtoolbox.org/models.php, there are at least 28
audio quality metrics exist as has been identified and proposed in [34] and
could be selectively applied to advance this study.

The data augmentation by means increasing the features through perceptual
measures, feature extraction, or by incorporating annotated data recorded from
physiological sensors for learning will allow us to explore the more optimal
model based on the combinatorial features and/or deep learning algorithm,
which is expected to provide more accurate results.

## Acknowledgment

The authors wish to thank Prof. Angelo Farina for the permission in using the
source materials and to the authors in [10] for providing the Matlab code of
AMBIQUAL. This research is funded by the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No.765911 RealVision.

## References

  * [1] Z. Akhtar and T. H. Falk, “Audio-visual multimedia quality assessment: A comprehensive survey,” _IEEE Access_ , vol. 5, pp. 21 090–21 117, 2017\.
  * [2] J. You, U. Reiter, M. M. Hannuksela, M. Gabbouj, and A. Perkis, “Perceptual-based quality assessment for audio–visual services: A survey,” _Signal Processing: Image Communication_ , vol. 25, no. 7, pp. 482–501, 2010\.
  * [3] H. T. Tran, C. T. Pham, N. P. Ngoc, A. T. Pham, and T. C. Thang, “A study on quality metrics for 360 video communications,” _IEICE Trans. Inf. Syst._ , vol. 101, no. 1, pp. 28–36, 2018.
  * [4] Z. Fei, F. Wang, J. Wang, and X. Xie, “QoE evaluation methods for 360-degree VR video transmission,” _IEEE J. Sel. Topics Signal Process_ , 2019.
  * [5] S. Fremerey, S. Göring, R. R. R. Rao, R. Huang, and A. Raake, “Subjective test dataset and meta-data-based models for 360° streaming video quality,” in _2020 IEEE 22nd International Workshop on Multimedia Signal Processing (MMSP)_. IEEE, 2020, pp. 1–6.
  * [6] S. Xie, Y. Xu, Q. Shen, Z. Ma, and W. Zhang, “Modeling the perceptual quality of viewport adaptive omnidirectional video streaming,” _IEEE Transactions on Circuits and Systems for Video Technology_ , 2019.
  * [7] M. S. Anwar, J. Wang, W. Khan, A. Ullah, S. Ahmad, and Z. Fei, “Subjective QoE of 360-Degree Virtual Reality Videos and Machine Learning Predictions,” _IEEE Access_ , vol. 8, pp. 148 084–148 099, 2020.
  * [8] M. A. Gerzon, “Periphony: With-height sound reproduction,” _J. Audio Eng. Soc._ , vol. 21, no. 1, pp. 2–10, 1973.
  * [9] T. Rudzki, I. Gomez-Lanzaco, J. Stubbs, J. Skoglund, D. T. Murphy, and G. Kearney, “Auditory localization in low-bitrate compressed ambisonic scenes,” _Appl. Sci._ , vol. 9, no. 13, p. 2618, 2019.
  * [10] M. Narbutt, J. Skoglund, A. Allen, M. Chinen, D. Barry, and A. Hines, “AMBIQUAL: Towards a quality metric for headphone rendered compressed ambisonic spatial audio,” _Appl. Sci._ , vol. 10, no. 9, p. 3188, 2020.
  * [11] M. Orduna, C. Díaz, L. Muñoz, P. Pérez, I. Benito, and N. García, “Video Multimethod Assessment Fusion (VMAF) on 360VR contents,” _IEEE Trans. Consum. Electron._ , vol. 66, no. 1, pp. 22–31, 2019.
  * [12] R. F. Fela, N. Zacharov, and S. Forchhammer, “Towards a perceived audiovisual quality model for immersive content,” in _Proceedings of the 2020 Twelfth International Conference on Quality of Multimedia Experience (QoMEX)_ , 2020, pp. 2472–7814.
  * [13] “Index of /Public/Jump-Videos,” http://pcfarina.eng.unipr.it/Public/Jump-Videos/, 2019.
  * [14] ITU-T Rec. P.910, _Subjective Video Quality Assessment Methods for Multimedia Applications_ , International Telecommunication Union Std., 2008.
  * [15] ITU-R Rec. BS.1387-1, _Method for Objective Measurements of Perceived Audio Quality_ , International Telecommunication Union Std., 2001.
  * [16] C. Sloan, N. Harte, D. Kelly, A. C. Kokaram, and A. Hines, “Objective assessment of perceptual audio quality using ViSQOLAudio,” _IEEE Trans. Broadcast._ , vol. 63, no. 4, pp. 693–705, 2017.
  * [17] M. Chinen, F. S. Lim, J. Skoglund, N. Gureev, F. O’Gorman, and A. Hines, “Visqol v3: An open source production ready objective speech and audio metric,” in _2020 Twelfth International Conference on Quality of Multimedia Experience (QoMEX)_. IEEE, 2020, pp. 1–6.
  * [18] E. Alshina, J. Boyce, A. Abbas, and Y. Ye, “JVET Common Test Conditions and Evaluation Procedures for 360 Degree Video,” _JVET document, JVET-G1030_ , 2017.
  * [19] Y. Sun, A. Lu, and L. Yu, “Weighted-to-spherically-uniform quality evaluation for omnidirectional video,” _IEEE Signal Process. Lett._ , vol. 24, no. 9, pp. 1408–1412, 2017.
  * [20] M. Yu, H. Lakshman, and B. Girod, “A framework to evaluate omnidirectional video coding schemes,” in _2015 IEEE International Symposium on Mixed and Augmented Reality_. IEEE, 2015, pp. 31–36.
  * [21] H. Lin, C. Huang, C. Li _et al._ , “AHG8: inter Digital’s Projection Format Conversion Tool Joint Video Exploration Team of ITU-T SG16 WP3 and ISO.” IEC JTC1/SC29/WG11, JVET-D0021. In: 5th Meeting, Geneva, 2017.
  * [22] V. Zakharchenko, K. P. Choi, and J. H. Park, “Quality metric for spherical panoramic video,” in _Optics and Photonics for Information Processing X_ , vol. 9970. International Society for Optics and Photonics, 2016, p. 99700C.
  * [23] Z. Li, A. Aaron, I. Katsavounidis, A. Moorthy, and M. Manohara, “Toward a practical perceptual video quality metric,” 2017. [Online]. Available: https://netflixtechblog.com/toward-a-practical-perceptual-video-quality-metric-653f208b9652
  * [24] Z. Li, C. Bampis, J. Novak, A. Aaron, K. Swanson, M. Anush, and J. D. Cock, “VMAF: The journey continues,” 2018. [Online]. Available: https://netflixtechblog.com/vmaf-the-journey-continues-44b51ee9ed12
  * [25] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” _IEEE Trans. Image Process._ , vol. 13, no. 4, pp. 600–612, 2004.
  * [26] Z. Wang, E. P. Simoncelli, and A. C. Bovik, “Multiscale structural similarity for image quality assessment,” in _The Thirty-Seventh Asilomar Conference on Signals, Systems & Computers, 2003_, vol. 2. IEEE, 2003, pp. 1398–1402.
  * [27] M. Kuhn _et al._ , “Building predictive models in r using the caret package,” _J Stat Softw_ , vol. 28, no. 5, pp. 1–26, 2008.
  * [28] M. Kuhn, J. Wing, S. Weston, A. Williams, C. Keefer, A. Engelhardt, T. Cooper, Z. Mayer, B. Kenkel, R. C. Team _et al._ , “Package ‘caret’,” _The R Journal_ , p. 223, 2020.
  * [29] T. M. Therneau, E. J. Atkinson _et al._ , “An introduction to recursive partitioning using the rpart routines,” Technical report Mayo Foundation, Tech. Rep., 1997.
  * [30] L. Breiman, “Random forests,” _Machine learning_ , vol. 45, no. 1, pp. 5–32, 2001.
  * [31] A. Karatzoglou, A. Smola, K. Hornik, and A. Zeileis, “kernlab-an s4 package for kernel methods in r,” _Journal of statistical software_ , vol. 11, no. 9, pp. 1–20, 2004.
  * [32] A. Karatzoglou, D. Meyer, and K. Hornik, “Support vector machines in r,” _Journal of statistical software_ , vol. 15, no. 9, pp. 1–28, 2006.
  * [33] E. Demirbilek and J.-C. Gregoire, “Machine learning–based parametric audiovisual quality prediction models for real-time communications,” _ACM Transactions on Multimedia Computing, Communications, and Applications (TOMM)_ , vol. 13, no. 2, pp. 1–25, 2017.
  * [34] M. Torcoli, T. Kastner, and J. Herre, “Objective measures of perceptual audio quality reviewed: An evaluation of their application domain dependence,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 29, pp. 1530–1541, 2021.


# Reactor antineutrino experiments

Haoqi Lu  Institute of High Energy Physics, Chinese Academy of Sciences
Yuquan Road 19-B, Beijing 100049,P.R. China
Member of Daya Bay collaboration
<EMAIL_ADDRESS>
(Day Month Year; Day Month Year)

###### Abstract

Neutrinos are elementary particles in the standard model of particle physics.
There are 3 flavors of neutrinos that oscillate among themselves. Their
oscillation can be described by a 3$\times$3 unitary matrix, containing three
mixing angles $\theta_{12}$, $\theta_{23}$, $\theta_{13}$, and one CP phase.
Both $\theta_{12}$ and $\theta_{23}$ are known from previous experiments.
$\theta_{13}$ was unknown just two years ago. The Daya Bay experiment gave the
first definitive non-zero value in 2012. An improved measurement of the
oscillation amplitude $\sin^{2}2(\theta_{13})$ = $0.090^{+0.008}_{-0.009}$ and
the first direct measurement of the $\bar{\nu}_{e}$ mass-squared difference
$\mid$$\Delta m^{2}_{ee}$$\mid$ = $(2.59^{+0.19}_{-0.20})\times 10^{-3}\rm
eV^{2}$ were obtained recently. The large value of $\theta_{13}$ boosts the
next generation of reactor antineutrino experiments designed to determine the
neutrino mass hierarchy, such as JUNO and RENO-50 .

###### keywords:

antineutrino,Daya Bay, $\theta_{13}$, mass hierarchy

PACS numbers: 14.60.Pq, 29.40.Mc, 13.15.+g

##  1 Introduction

The neutrino is a fundamental particle and was first discovered in 1956 by
Cowan and Reines [1],[2]. In the last few decades, it has been proven that the
observed neutrino oscillations can be described in a 3-flavor neutrino
framework. A parameterization of the standard Pontecorvo-Maki-Nakagawa-Sakata
(PMNS) matrix describing the unitary transformation relating the mass and
flavor eigenstates, defines the three mixing angles ($\theta_{23}$,
$\theta_{12}$, and $\theta_{13}$) and one charge-parity(CP)-violating phase
[3],[4]. $\theta_{12}$ is about 34º and determined by solar and reactor
neutrino experiments. $\theta_{23}$ is about 45º and determined by atmospheric
and accelerator neutrino experiments. An upper limit of the last unknown angle
$\theta_{13}$ was given by CHOOZ $\sin^{2}(2\theta_{13})$$<$0.15 at 90%
confidence level (C.L.)[5]. It was hinted to be non-zero by recent results
from the T2K [6], MINOS [7] and Double Chooz experiments[8]. The value of
$\theta_{13}$ will guide the designs of future experiments for the
measurements of the mass hierarchy and CP-violation.

Long-baseline accelerator experiments have limited sensitivity to the
$\theta_{13}$ mixing angle due to dependencies on the yet unknown mass
hierarchy and CP-violating phase. Reactor antineutrino experiments can provide
a clear and accurate measurement of $\theta_{13}$, due to their pure
antineutrino source, clear signal and independence of the CP phase and matter
effects. For reactor-based antineutrino experiments, $\theta_{13}$ can be
determined in terms of the survival probability of $\bar{\nu_{e}}$ at certain
distances from the reactors,

$P(\bar{\nu}_{e}\rightarrow\bar{\nu}_{e})=1-\cos^{4}\theta_{13}\sin^{2}2\theta_{12}\sin^{2}\Delta_{21}-\sin^{2}2\theta_{13}(\cos^{2}\theta_{12}\sin^{2}\Delta_{31}+\sin^{2}\theta_{12}\sin^{2}\Delta_{32})$
(1)

where $\Delta_{ji}$ $\equiv$ 1.267$\Delta^{2}_{ji}(\rm eV^{2})$ L(m)/E(MeV), E
is the $\bar{\nu}_{e}$ energy in MeV and L is the distance in meters between
the $\bar{\nu}_{e}$ source and the detector. $\Delta m^{2}_{ji}$ is the
difference between the mass-squares of the mass eigenstates $\nu_{j}$ and
$\nu_{i}$( $\Delta m^{2}_{ji}$ = $\Delta m^{2}_{j}$ \- $\Delta m^{2}_{i}$).
Since $\Delta m^{2}_{21}$ $\ll$ $\mid$ $\Delta
m^{2}_{31}$$\mid$$\approx$$\mid$$\Delta$ $m^{2}_{32}$$\mid$[9], the short
distance ($\sim$ km) reactor $\bar{\nu}_{e}$ oscillation is mainly determined
by the $\Delta_{3i}$ terms.

Most reactor antineutrino oscillation experiments measure antineutrino events
via the inverse beta decay(IBD) reaction $\bar{\nu}+p\to e^{+}+n$. The IBD
reaction is characterized by two time correlated events, the prompt signal
coming from the production and subsequent annihilation of the positron, and
the delayed signal from the capture of the neutron in the liquid scintillator.
For the Daya Bay experiment, 0.1% gadolinium (Gd)-doped liquid scintillator is
used to increase the capture cross section of thermal neutrons on Gd and
reduce the capture time (about 30 $\rm\mu$s) to suppress accidental
coincidence backgrounds. The total energy of gammas from neutron capture on Gd
is about 8 MeV, which is much higher than the energy range of background
radioactivity. The observable antineutrino spectrum is shown in Fig.1.

Figure 1: The observable $\bar{\nu}$ spectrum, a product of the antineutrino
flux from reactors and cross section of IBD reaction.

Due to the small value of $\theta_{13}$, precise measurements are required to
reduce the experimental uncertainties. In previous reactor antineutrino
experiments, such as Palo Verde [10] and CHOOZ [5], single detectors were
utilized for antineutrino detection. Their measurements had large uncertainty
from errors related to detector response and antineutrino flux predictions.
There are three ongoing reactor experiments : Daya Bay, Double Chooz, and
RENO. They use multiple detectors at different baselines to reduce correlated
uncertainties by a relative measurement. A comparison of the three experiments
is shown in Table 1\. The baselines of their far detectors are 1.65, 1.05 and
1.44km, respectively. [11].

The following content of this paper mainly focuses on the Daya Bay experiment
and future prospects of reactor antineutrino experiments.

Comparison of three reactor experiments  Exp. Power Det.(t) Overburden 3y
Sens. (GW) (N/F) (m.w.e)N/F (90%CL) Double Chooz 8.5 8/8 120/300 0.03 RENO
16.5 16/16 120/450 0.02 DYB 17.4 80/80 250/860 0.008

Figure 2: Layout of the Daya Bay experiment. There are six reactor cores (D1,
D2, L1, L2, L3, L4) and 3 experimental halls (EH1, EH2, EH3).

##  2 Daya Bay experiment

The Daya Bay experiment is designed to explore the unknown value of
$\theta_{13}$ by measuring the survival probability of electron antineutrinos
from the nuclear reactors in Daya Bay, China. The Daya Bay Nuclear Power Plant
complex, one of the 5 most prolific sources of reactor neutrinos in the world,
consists of 6 reactors(see Fig.2) producing 17.4 GW of total thermal power.
The goal of the experiment is to measure $\theta_{13}$ with a sensitivity of
$\sin^{2}2\theta_{13}<0.01$ at a 90% C.L. In order to achieve such high
sensitivity, the experiment is optimized in several aspects. Multiple sites
(one far experimental hall (EH3) and two near experimental halls (EH1,EH2))
are used to effectively cancel the flux uncertainty by relative measurements.
The experiment employs 8 identically designed detectors to decrease detector
related errors. As shown in Fig.2, only 6 detectors were deployed before
August. 2012. All 8 detectors were installed by October, 2012\. The detectors
are installed underground to reduce the cosmic-ray muon flux. Each site has
redundant muon detectors (water Cerenkov detectors and resistive plate
chambers (RPCs)) for muon identification.

###  2.1 Detector

####  2.1.1 Antineutrino detector

Figure 3: Antineutrino detector of Daya Bay experiment.

The antineutrino detectors (ADs) are filled with Gd doped liquid scintillator
for antineutrino event detection. The experiment uses 8 functionally identical
detectors (2 at EH1, 2 at EH2 and 4 at EH3), which are cylindrical stainless
steel vessels(SSV) with a 5 m diameter and 5 meter height (Fig.3). Each
detector is instrumented with 192 8-inch photomultiplier tubes (PMTs). The
detector has a three-zone structure, including a Gd-doped liquid scintillator
(GdLS) zone, liquid scintillator (LS) zone and mineral oil (MO) zone. The
inner region is the primary target volume filled with 0.1% Gd-LS. The middle
layer is LS and act as a gamma catcher, and the outer layer is filled with MO
which reduces the impact of radioactivity from PMT glass and the SSV. Two
acrylic tanks are used to separate each layer. There are two reflectors at the
top and bottom of an AD, which are laminated with ESR (VikuitiTM Enhanced
Specular Reflector Film) film sealed between two 1-cm thick acrylic panels.The
reflectors improve light collection and uniformity of detector response.

Three automated calibration units (ACUs) are located on the top of each AD
(Fig.3). The ADs are calibrated periodically by deploying LEDs and radioactive
sources inside the detectors. LEDs are used for timing and PMT gain
calibration. Energy calibration is performed with radioactive sources. Based
on the first three months of data, the first pair of ADs in EH1 are shown to
be functionally identical [12].

####  2.1.2 Muon veto system

Figure 4: Veto system of Daya Bay experiment. It includes the inner water
shield, outer water shield and top RPC detector.

The muon veto system of the Daya Bay experiment is shown in Fig.4 . The
antineutrino detectors are immersed in an octagonal pool of ultrapure water.
The pool is divided into outer and inner volumes by Tyvek sheets. Both volumes
are instrumented with PMTs as active muon detectors via muon Cherenkov light.
The Tyvek sheet with very high reflectivity($>$95%) can increase light
collection efficiency. The outer layer of the water pool is 1 m thick and the
thickness of the inner layer is $>$1.5 m. At least 2.5 m of water surrounds
each AD to shield against ambient radioactivity. 288 8-inch PMTs are installed
in each near hall pool and 384 in the Far Hall. There is a water circulation
and purification system in each hall to maintain water quality. The tops of
the water Cherenkov detectors are covered by 4 layers of RPCs. RPCs are
gaseous detectors with resistive electrodes operating in streamer mode at Daya
Bay[13]. RPC signals are read out by external strips, with a position
resolution of about 8 cm. There are 54 modules in each near hall and 81
modules in the Far Hall. The designed efficiency is $>$99.5% with uncertainty
less than 0.25% by combining the water Cherenkov and RPC detectors. From muon
data analysis, water Cherenkov detector efficiency is $>$99.7% for long track
muons[12], which is better than the design requirement.

###  2.2 Data analysis

####  2.2.1 Event selection

After energy calibration and event reconstruction, we select the IBD events
under the following criteria:

PMT light emission events will be rejected by the PMT flasher cut[12].

The energy of prompt and delayed candidates satisfy 0.7 MeV $<$ Ep $<$ 12 MeV
and 6.0 MeV $<$ Ed $<$ 12 MeV, where Ep is the prompt signal energy and Ed is
the delayed signal energy.

The delayed signal (neutron capture) time satisfies: 1$\mu$s $<$t $<$
200$\mu$s, where t is the time between the prompt and delayed signal.

If a muon goes through the water pool and fires $>$12 PMTs, it will be treated
as a ’water pool muon’. The rejection time window of water pool muons is
[-2$\mu$ s, 600$\mu$s].

If the visible energy in the AD is ($>$20 MeV), it will be treated as an ’AD
muon’. The rejection time window of AD muon is [-2$\mu$s, 1400$\mu$s].

If the visible energy in the AD is $>$2.5 GeV, the rejection time window is
[-2 $\mu$ s,0.4s]

There must be no additional prompt-like signals 400$\mu$s before the delayed
neutron signal and no additional delayed-like signal 400$\mu$s after the
delayed neutron signal. As Fig. 5 shows, the events within the dash lines are
the IBD candidate events after selection.

Figure 5: IBD candidate events are shown within the dashed lines

####  2.2.2 Background

Background control is crucial for a precise measurement in reactor
antineutrino experiments. Li9/He8 is the main background from cosmic-ray
muons. It is a correlated background that can mimic prompt signals with a
beta-decay signal and delayed signals with neutron capture. This background is
directly measured by fitting the distribution of IBD candidates versus time
since the last muon. The background to signal ratio (B/S) of Li9/He8 is
$\sim$0.3$\pm$0.1%. Another background from muons is the fast neutron
background. Energetic neutrons from muons can give a prompt signal from proton
recoil and a delayed signal from the neutron’s capture on Gd. By relaxing the
prompt energy range from 12 MeV to 50 MeV in the IBD events selection, the
fast neutron spectrum in the range of $<$12MeV can be extrapolated from the
higher energy distribution. The fast neutron spectrum and rate is cross-
checked with fast-neutron samples tagged by the muon veto system. The B/S
ratio of fast neutrons is 0.1% (0.06%) for the near(far) sites.

Two uncorrelated signals can accidentally mimic an IBD signal. This accident
background can be precisely measured from data. The B/S of accidentals is 1.5%
(4%) for the near(far) sites. The error in the accidental rate is 1% (4%).

The 238U, 232Th, 227Ac decay chains and 210Po produce a 13C($\alpha$ n)16O
backgrounds. This background is estimated by Monte Carlo(MC).The B/S of this
background is 0.01% (0.05%) for the near(far) sites. Am-C calibration sources
inside the ACUs on top of each AD constantly emit neutrons, which would
occasionally mimic IBD signals by scattering inelastically with nuclei in the
shielding material (emitting gamma rays) before being captured on a metal
nucleus, such as Fe, Cr, Mn or Ni (releasing more gamma rays). MC is also used
to estimate the rate of this background. The background rate and shape are
constrained by taking data with a temporary intense source on top of the AD.
The B/S ratio of Am-C sources is 0.04% (0.35%) for the near(far) sites.

####  2.2.3 Recent results

Figure 6: Allowed regions for the neutrino oscillation parameters
$\sin^{2}(2\theta_{13})$ and $\mid$$\Delta m^{2}_{ee}$$\mid$ at different
confidence levels(solid regions). The best estimates of the oscillation
parameters are shown by the black dot. The adjoining panels show the
dependence of $\Delta\chi^{2}$ on $\mid$$\Delta m^{2}_{ee}$$\mid$(right) and
$\sin^{2}(2\theta_{13})$ (top). The black square and dashed curve represent
the rate-only result. The dashed horizontal line represents the MINOS
$\mid$$\Delta m^{2}_{\mu\mu}$$\mid$ measurement17.

The Daya Bay experiment first presented the discovery of a non-zero
$\theta_{13}$(5.2$\sigma$) in 2012[14]. The most recent analysis[16] is based
on the full six-AD data sample from Dec. 24, 2011 to July 28, 2012, which
includes more than 300,000 IBD interactions. The rate-only analysis yields
$\sin^{2}$$(2\theta_{13})$ = 0.090 $\pm$ 0.010 with $\chiup^{2}$/NDF = 0.6/4.
Spectral information is also used after applying an energy nonlinearity
correction to the positron spectrum. The oscillation parameters can be
extracted from a fit taking into account both rate and spectral information.
In Eq. 1, the short baseline ($\sim$ km) reactor $\bar{\nu}_{e}$ oscillation
is mainly determined by the $\Delta_{3i}$ terms. We use the following
definition[15] of the effective mass-squared difference $\sin^{2}\Delta ee$
$\equiv$$\cos^{2}\theta_{12}\sin^{2}\Delta_{31}+\sin^{2}\theta_{12}\sin^{2}\Delta_{32}$.
The best-fit values are $\sin^{2}(2\theta_{13})$ = $0.090^{+0.008}_{-0.009}$
and $\mid$$\Delta m^{2}_{ee}$$\mid$ =
$(2.59^{+0.19}_{-0.20})$$\times$$10^{-3}\rm eV^{2}$ with $\chiup^{2}$/NDF =
163/153 (68.3% C.L.). The 68.3%, 95.5%, and 99.7% C.L. allowed regions in the
$\mid$$\Delta m^{2}_{ee}$$\mid$ vs. $\sin^{2}$$(2\theta_{13})$ plane are shown
in Fig.6.

The relative deficit and spectral distortion observed between far and near ADs
at Daya Bay give the first independent measurement of $\mid$$\Delta
m^{2}_{ee}$$\mid$ = $(2.59^{+0.19}_{-0.20})$$\times$$10^{-3}\rm eV^{2}$ and
the most precise measurement of $\sin^{2}(2\theta_{13})$ =
$0.090^{+0.008}_{-0.009}$ to date.

###  2.3 Future plan

The sensitivity of $\sin^{2}(2\theta_{13})$ versus time for Daya Bay is shown
in Fig. 7. The projected uncertainty of $\sin^{2}(2\theta_{13})$ is less than
4% after 3 years of data taking. Daya Bay is also expecting to measure
$\mid$$\Delta m^{2}_{ee}$$\mid$ in complementary precision to accelerator
neutrino experiments. Fig.8 shows the $\mid$$\Delta m^{2}_{ee}$$\mid$ error
versus data-taking time. The error of $\mid$$\Delta m^{2}_{ee}$$\mid$ will be
$<$ $0.10$$\times$$10^{-3}\rm eV^{2}$ after 3 years of data-taking.

Figure 7: Error of $\sin^{2}(2\theta_{13})$ vs time. The red circle is this
analysis result. The current data-taking time is above 600 days .

Figure 8: Error of $\mid$$\Delta m^{2}_{ee}$$\mid$ vs time. The red circle is
this analysis result. The current data-taking time is above 600 days. The
horizontal dashed line is the error of $\mid$$\Delta m^{2}_{\mu\mu}$$\mid$
from MINOS.

###  2.4 RENO and Double Chooz status

RENO updated their result in March, 2013[18]. They gave
$\sin^{2}(2\theta_{13})$ =0.100$\pm$0.010(stat)$\pm$0.015(syst) and observed a
clear rate deficit (7.1% reduction) between far and near sites. RENO aims to
achieve 7% uncertainty in $\sin^{2}(2\theta_{13})$ and suppress their systemic
error to 0.5%. Double Chooz[19] provided their new results from a combined fit
of Gd capture and Hydrogen capture with $\sin^{2}(2\theta_{13})$
=0.109$\pm$0.035. For Double Chooz, the near detector is under construction
and will start running in 2014.

##  3 Future prospects

The next generation of neutrino oscillation experiments will mainly focus on
the mass hierarchy and CP phase. The mass hierarchy can be determined by a
precisely energy spectrum measurement in reactor neutrino oscillations, which
is independent of the CP-violating phase. Fig.9 shows the L/E distribution of
reactor antineutrino oscillations. It is found that frequencies of oscillation
between the normal hierarchy and inverted hierarchy have some difference. The
two oscillation frequency components are driven by $\Delta m^{2}_{31}$ and
$\Delta m^{2}_{32}$ , respectively [20]-[23]. The difference is very sensitive
to detector energy resolution. According to the study in Ref. Zhanl2, for a
single detector at a baseline of 58 km and with a 35 GW reactor, the
probability to determine the sign of the hierarchy has a significant
difference between 3% and 4% energy resolution as shown in Fig.10. Another big
effect on mass hierarchy sensitivity is potentially from multiple reactor
baselines, which are required to be within 500 m[25]. A large detector with
good energy resolution (3%) and equal baselines from powerful reactor cores is
required for a mass hierarchy determination.

Figure 9: Neutrino oscillation difference between normal hierarchy and
inverted hierarchy.

Figure 10: The probability to determine the sign of the hierarchy vs. the
number of detected events. Different energy resolutions are used to estimate
the effect on the probability.

###  3.1 The Jiangmen Underground Neutrino Observatory

The Jiangmen Underground Neutrino Observatory (JUNO) is located at Kaiping,
Jiangmen City, in Guangdong Province, China, as shown in Fig.11, 53 km from
the Yangjiang and Taishan nuclear power plants. The total thermal power of the
reactors is 36 GW. The detector will be constructed deep underground to reduce
the cosmic-ray muon flux, with an overburden of about 700 m of rock. The
detector will be filled with 20 ktons of liquid scintillator with 3% energy
resolution (at 1MeV). The major goal of JUNO is to determine the mass
hierarchy by precisely measuring the energy spectrum of reactor antineutrinos.
However, it is a multi-purpose experiment that can also measure neutrino
oscillation parameters, study atmospheric, solar and geo-neutrinos, and
perform other exotic searches.

Fig.12 is the detector concept of JUNO. There are around 15,000 20-inch PMTs
installed in the central detector for the 3% energy resolution. 6 ktons of
mineral oil is used to shield the scintillator from PMT glass radioactivity.
The muon veto system includes an outer water Cherenkov detector and a top
tracking system. As for the physics prospects of JUNO[25], if we take into
account the spread of reactor cores, uncertainties from energy non-linearity,
etc, the mass hierarchy sensitivity with 6 years of data-taking can reach
$\Delta$$\chi^{2}$$>$ 9 with a relative measurement and $\Delta$$\chi^{2}$
$>$16 with an absolute $\Delta m_{\mu\mu}^{2}$ measurement from accelerator
neutrino experiments. Civil construction has begun and will complete in 2019.
Liquid scintillator filling and data taking will begin around 2020.

Figure 11: Experiment site of JUNO. The Taishan and Yangjiang reactor
complexes are used for JUNO. The previous site is not considered in that the
third reactor complex (Lufeng) is being planed.

Figure 12: Detector concept of JUNO, including the central liquid scintillator
detector and outer muon veto systems.

###  3.2 RENO 50

RENO-50 was proposed in South Korea[26]. The detector will be constructed
underground and consist of 18 ktons of ultra-low-radioactivity liquid
scintillator and 15,000 20-inch PMTs, about 50 km away from the Hanbit nuclear
power plant. The scientific goals of the experiment include a high precision
measurement of $\theta_{12}$ and $\mid$$\Delta m^{2}_{21}$$\mid$,
determination of the mass hierarchy, the observation of neutrinos from
reactors, the Sun, the Earth, supernovae, and any possible stellar objects.
The total budget is 100 million dollars for 6 years of construction. Facility
and detector construction will be started from 2013 and finished in 2018\. The
experiment data taking will start in 2019.

##  4 Summary

Reactor-based antineutrino experiments have obtained many excellent
achievements in recent years. Daya Bay measured the non-zero $\theta_{13}$
with great precision, together with other experiments. The most recent results
from Daya Bay provide the first independent measurement of $\mid$$\Delta
m^{2}_{ee}$$\mid$ and the most precise measurement of $\sin^{2}(2\theta_{13})$
to date. The precision on $\sin^{2}(2\theta_{13})$ will be improved to 4% in
the future. The large value of $\theta_{13}$ boosts the next generation of
neutrino oscillation experiments to determine the neutrino mass hierarchy and
measure the CP violation phase. Reactor-based antineutrino experiments will
continue to play an important role in the mass hierarchy determination and
precise measurements of oscillation parameters.

## Acknowledgments

The article is supported by the National Natural Science Foundation of China
(Y3118G005C). I would like to acknowledge my Daya Bay collaborators for useful
suggestion and comments, especially Dr. Yufeng Li, Logan Lebanowski and Viktor
Pec who helped me a lot to improve this article.

## References

  * [1] C. L. Cowan et al. Science  124, 103 (1956);
  * [2] F. Reines and C. L. Cowan, Jr., Nature 178, 446 (1956).
  * [3] B. Pontecorvo, Sov. Phys. JETP 6, 429 (1957) and 26, 984 (1968);
  * [4] Z. Maki, M. Nakagawa, and S. Sakata, Prog. Theor. Phys. 28, 870 (1962).
  * [5] CHOOZ Collab. ( M. Apollonio et al.), Phys. Lett. B 466, 415 (1999).
  * [6] T2K Collab. (K. Abe et al.), Phys. Rev. Lett. 107, 041801 (2011).
  * [7] MINOS Collab. (P. Adamson et al.), Phys. Rev. Lett. 107, 181802 (2011).
  * [8] Double Chooz Collab. ( Y. Abe et al.), Phys. Rev. Lett. 108, 131801 (2012).
  * [9] J. Beringer et al. (Particle Data Group), Phys. Rev. D 86, 010001 (2012), Section 13.
  * [10] F. Boehm et al., Phys. Rev. Lett. 84, 400 3764(2000).
  * [11] J. Cao, plenary talk at the 36th International Conference on High Energy Physics, Melbourne, Jul. 2012.
  * [12] DayaBay Collab. ( F. An et al.), NIM A 685, 78-97 (2012).
  * [13] L. Ma et al.NIM A  659, 154(2011); J. Xu et al., Chin. Phys. C 35,844 (2011).
  * [14] Daya Bay Collab. ( F. Anet al.), Phys. Rev. Lett. 108, 171803 (2012).
  * [15] The Daya Bay $\mid$$\Delta m^{2}_{ee}$$\mid$ definition is consistent with H. Minakata et al., Phys. Rev. D bf 74, 053008 (2006) and H. Minakata et al., Phys. Rev. D 76, 053005 (2007).
  * [16] Daya Bay Collab.( F. Anet al.), http://arxiv.org/abs/1310.6732
  * [17] (MINOS Collaboration),(P. Adamsonet al.) Phys. Rev. Lett. 110, 51801 (2013).
  * [18] Seon-Hee Seo, talk on NuTel2013, March 2013
  * [19] Ralitsa Sharankova, talk on NuFact2013, Aug,2013, Beijing
  * [20] S. T. Petcov and M. Piai Phys. Lett. B  533, 94 (2002);
  * [21] S. Choubey, S. T. Petcov, and M. Piai Phys. Rev. D  68, 113006 (2003).
  * [22] J. Learned et al., hep-ex/0612022
  * [23] L. Zhan, Y. Wang, J. Cao, L. Wen,Phys. Rev. D  78, 111103 (2008); Phys. Rev. D  79, 073007 2009
  * [24] L. Zhan, Neutrino Oscillation Workshop 2012,September 9-16, 2012 Italy
  * [25] Yu-Feng Li, Jun Cao, Yifang Wang, and Liang Zhan,Phys. Rev. D  88, 013008 (2013).
  * [26] Soo-Bong Kim, talk at the international workshop on RENO-50, June 13-14, 2013, Seoul.


# Magnetism of competing high-spin/low-spin states in Ba2NiO2(AgSe)2 and
related two-orbital two-electron systems

Rajiv R. P. Singh  Department of Physics, University of California Davis, CA
95616, USA

###### Abstract

We discuss the thermodynamic and magnetic properties of a competing high-
spin/low-spin two-orbital two-electron model on a square-lattice possibly
relevant to the nickelates such as Ba2NiO2(AgSe)2 (BNOAS). We focus on
parameter regimes where a high-spin ($S=1$) and a low-spin ($S=0$) state are
energetically close to each other and discuss various exchange processes in
such a system. The model we study is a variant of, but different from the
Kondo-necklace model proposed for the system by Jin et al jin . Although there
are similarities between the two models in terms of the ground-state phases
and their symmetries, the detailed properties of the different phases and the
phase transitions are entirely different and should be easy to distinguish
from an experimental point of view.

##  I I. Introduction

Following the discovery of high-temperature superconductivity in the cuprates,
many new families of materials have been synthesized with competing magnetic
and superconducting phases. In particular, several iron based superconductors
are known to have competing magnetic and superconducting phases and, at
stoichiometry, are often modelled as multi-orbital spin-one systems qimiao .

The discovery of superconductivity in the nickelates has brought renewed
interest to these systems as well hwang . Although much recent excitement
comes from doping the Ndn+1NinO2n+2 family of materials where nickel ions have
a Ni1+ ionic state with $d^{9}$ electron configuration providing an analog of
the cuprates, there are other nickelate families with the more usual $Ni^{2+}$
ionic state and $d^{8}$ electron configuration klee2 ; botana . Our work is
focused on competing high-spin and low-spin states in these latter systems.

Jin et al jin ; klee  have recently studied the material Ba2NiO2(AgSe)2
(BNOAS) by Density Functional Theory (DFT) and found a rather striking result.
Ni atoms have $d^{8}$ electronic configuration, with competing spin-zero and
spin-one states, where the lowest energy state has spin-zero. Depending on the
value of the Hubbard repulsion $U$, this lowest energy singlet state can be
made of two electrons occupying two different orbitals ($d_{x^{2}-y^{2}}$ and
$d_{z^{2}}$). A lower energy of this singlet compared with the triplet amounts
to an antiferromagnetic intra-atomic exchange, which goes against the Hund’s
rules. Jin et al have proposed an effective Kondo-necklace model, where the
two electrons on an atom are exchange coupled by an antiferromagnetic intra-
atomic exchange. In addition, one of the electrons (in the $d_{x^{2}-y^{2}}$
orbital) has an exchange interaction with electrons on neighboring atoms in
the plane, whereas the electron in the other state (in the $d_{z}^{2}$
orbital) provides an analog of the local spin in the Kondo models. To leading
order, its only coupling is to the other spin in the same atom. This leads to
a Kondo-necklace model with two spin-half operators $\vec{S}_{1i}$ and
$\vec{S}_{2i}$ at each site of the square-lattice and Hamiltonian:

$H=J\sum_{\langle
i,j\rangle}\vec{S}_{1i}\cdot\vec{S}_{1j}+J_{l}\sum_{i}\vec{S}_{1i}\cdot\vec{S}_{2i},$
(1)

where the first sum runs over nearest-neighbor pairs. At small $J/J_{l}$ the
system is in a non-magnetic singlet phase, where the two spins on each atom
combine to form a singlet. At large $J/J_{l}$, the system is in the Neel phase
with antiferromagnetic order. The two phases are separated by a quantum
critical point in the universality class of the 3-dimensional classical
Heisenberg model brenig .

Here we will focus on an alternative model for this insulating system with
$d^{8}$ electron configuration and two d-holes in two orbitals with competing
high-spin and low-spin states. In a two-orbital system, two electron (or two-
hole) configuration leads to one triplet and 3 singlets. Two of the singlets
correspond to the two electrons occupying the same orbital, while a singlet
and a triplet correspond to electrons occupying different orbitals. We will
assume that the intra-atomic Hund’s coupling is the normal ferromagnetic one
so that the triplet state is lower in energy when the two electrons occupy
different orbitals. However, Hubbard repulsion $U$, crystal field splittings
and ligand interactions are such that one of the singlets formed from holes
occupying the same orbital (plus ligands) is lower in energy than the triplet
state. Since intra-atomic Hund’s exchange is typically of order $1$ eV, we
will assume that the other two singlet states are significantly higher in
energy and not relevant to the low energy properties. A similar competing low-
spin/high-spin scenario has been discussed recently for Ni $d^{8}$ ions in the
material NiO2 by Jiang et al sawatzky .

##  II II. Exchange processes and the Model Hamiltonian

The overlap between orbitals on neighboring atoms leads to inter-atomic
exchange processes. A general two-band Hubbard model for the nickelates was
studied by Hu and Wu wu , who used perturbation theory to argue for a Kugel-
Khomskii type spin-orbital model khomskii  with many possible terms.

Here we take as our starting point for deriving the effective low energy
Hamiltonian two key assumptions motivated by the work of Jin et al jin  and
Jiang et al sawatzky : (i) One spin-zero and one spin-one state forms the low
energy subspace of our system and all others states are at significantly
higher energy. Thus, all inter-atomic exchanges must be projected on to these
low energy states of each atom. (ii) Inter-atomic hopping is much larger for
spins in one of the orbitals (e.g. $d_{x^{2}-y^{2}}$) than the other (e.g.
$d_{z^{2}}$). These considerations lead to the following inter-atomic
exchanges in order of decreasing magnitude:

1\. Spin-one Heisenberg exchange J: The largest inter-atomic exchange is a
Heisenberg coupling J between neighboring atoms when both are occupied by
spin-one states. For this exchange, overlap between $d_{x^{2}-y^{2}}$ orbitals
is sufficient. The difference from the Kondo necklace model of Jin et al is
that the projection to the low energy subspace means that the Heisenberg
coupling is between spin-one states on neighboring atoms. It is not a coupling
for spin-half or just one of the two spins. That is forbidden by the large
intra-site Hund’s coupling.

2\. Triplet hopping process K: When there is a spin-zero on one site and a
spin-one on a neighboring site, the two sites can exchange their spin states.
This is a much weaker process than the Heisenberg exchange as it requires
overlap of both (e.g. $d_{x^{2}-y^{2}}$ and $d_{z^{2}}$ ) orbitals on
neighboring atoms and thus may require either a lattice distortion or a
higher-order process. In our model we will ignore this term in the numerical
study. Qualitatively, it is important in the non-magnetic phase where even a
small value is relevant to the dispersion of spin excitations.

3\. Triplet-pair creation/destruction process W: A pair of triplets with total
spin-zero can turn into a pair of spin-zero states on neighboring sites and
vice versa. This process will require cross overlap between neighboring
$d_{x^{2}-y^{2}}$ and $d_{z^{2}}$ orbitals. We will assume this process is the
weakest and will not discuss it any further.

Thus, our low energy states on an atom are characterized by $n_{i}=0$ or $1$,
with $n_{i}=0$ corresponding to spin-zero and $n_{i}=1$ to the spin-one state.
In addition, when $n_{i}=1$, there are three spin states corresponding to
$S^{z}_{i}=0,\pm 1$. Thus, there are 4-states per atom. The effective
Hamiltonian can be expressed in terms of $n_{i}$ and the quantum spin-one
operators $\vec{S}_{i}$ as:

$H=\Delta\sum_{i}n_{i}+J\sum_{\langle
i,j\rangle}n_{i}n_{j}\vec{S}_{i}\cdot\vec{S}_{j},.$ (2)

We will assume that the parameter $\Delta$ is greater than zero so that the
lowest energy state on each atom has spin-zero.

##  III III. Ground State Phases of the model

Ground state phases and properties of the Kondo-necklace model have been
studied previously brenig . Certain aspects of the ground state phase diagram
of the present model can be obtained using results from the literature singh ;
hamer ; hamer2 . Like the Kondo-necklace model, this model has two phases. For
small $J/\Delta$ the model is in a non-magnetic singlet phase, where each atom
is occupied by electrons in a spin-zero state. For large $J/\Delta$ the model
is in a spin-one Neel state with long-range antiferromagnetic order. However,
there are key differences in the properties of the two phases and the phase-
transition between the two models.

Neel phase: The Neel phase is characterized by long-range antiferromagnetic
order and magnon excitations. The key difference between the models lies in
the behavior of the uniform susceptibility. In the present model, deep in the
Neel phase, the susceptibility has the temperature dependence of an
antiferromagnet, that goes down below a temperature of order $J$, when short-
range antiferromagnetism develops. At low temperatures, as the spins are
locked antiferromagnetically, they have only a weak transverse susceptibility
singh ; hamer ; hamer2 . In contrast, in the Kondo-necklace model, as one gets
away from the critical point, the low temperature uniform susceptibility
becomes large. This is because the local spins become nearly free and thus
develop almost Curie-like susceptibility brenig .

Non-magnetic singlet phase: In the Kondo necklace model, the singlet phase has
a dispersive triplet excitation, where the dispersion band-width is set by the
inter-atomic exchange energy $J$. As the phase transition is approached, the
excitation gap will go to zero at the antiferromagnetic wave-vector and at
$k=0$, with neutron-scattering spectral weight strongly concentrated near the
anti-ferromagnetic wave-vector. In contrast, for the present model, the
triplet will be static and completely localized in the absence of the triplet-
hopping process discussed earlier. It is the weaker triplet hopping process
$K$ that will lead to a dispersion for the excitations. Thus, one would expect
the dispersion to be relatively weak and being unrelated to the exchange $J$
it should have a minimum at $k=0$. In addition, there will be a clustering of
the spinful states as when the spin-one states are clustered together their
energy will be lowered than when they are separated.

Ground state Phase Transition: In the Kondo necklace model there is a quantum
critical point separating the non-magnetic phase from the antiferromagnet. In
contrast, in the present model, there will be a first order phase transition.
Well before the gap to spin-one excitations closes, the exchange energy will
lower the energy of the Neel state below that of singlets. Thus the system
will undergo a phase separation and a first order ground-state phase-
transition between the two phases. Since the ground state energy of the spin-
one model is known accurately singh  to be approximately $-2.327J$ per site,
the phase transition will happen at $\Delta=2.327J$.

Figure 1: Sketch of the phase diagram for the model. A line of first order
phase transitions, shown by dashed line, separates the Renormalized Classical
phase (RC) from the non-magnetic spin-zero phase. The line of phase
transitions terminates at an Ising critical point.

Finite temperature phase transition: In the present model there will be a
first order liquid-gas type phase transition at finite temperatures ending at
a critical point. Even though the spin degrees of freedom cannot order at
finite temperatures, the Ising variables $n_{i}$ can have a true phase
transition. Since the transition will occur at temperatures well below the
mean-field Neel temperature of $8J/3$, by that temperature short-range order
will be well developed. Hence, we can study the transition to a good
approximation by replacing the Heisenberg coupling in the Hamiltonian by its
low temperature expectation value which will be only weakly temperature
dependent. Thus our model becomes:

$H=\Delta\sum_{i}n_{i}-{\tilde{J}}\sum_{\langle i,j\rangle}n_{i}n_{j},$ (3)

where ${\tilde{J}}$ is minus the expectation value of
$<\vec{S}_{i}\cdot\vec{S}_{j}>$ on a nearest neighbor bond. It is weakly
temperature dependent at low temperatures singh

${\tilde{J}}(T)={2.327\over 2}J-aT^{3}.$ (4)

The $T^{3}$ term comes from the reduction in energy due to the magnons.

By the standard mapping

$n_{i}=(s_{i}+1)/2,$ (5)

this can be expressed in terms of the Ising variables $s_{i}=\pm 1$ as:

$H=-J_{I}\sum_{\langle i,j\rangle}s_{i}s_{j}+h\sum_{i}s_{i}+C,.$ (6)

where

$J_{I}={1\over 4}{\tilde{J}},$

and,

$h={\Delta-2{\tilde{J}}\over 2}.$

For the Ising model in Eq. 6, the phase transition line is given by $h=0$,
that is $\Delta=2{\tilde{J}}$. At $T=0$, the transition is at $\Delta=2.327J$.
The line of first order transitions will terminate at a critical point, which
can be obtained from the solution of the 2D Ising model to be at
$T={2J_{I}\over\ln{1+\sqrt{2}}}$. In the $\Delta/J$ and $T/J$ plane the phase
boundary gradually bends down due to the temperature dependence of
${\tilde{J}}$ and the critical point is roughly at $\Delta/J=2.2$,
$T_{c}/J=0.6$. The phase diagram is sketched in Fig 1.

##  IV IV. Numerical Study of finite temperature properties

Figure 2: Antiferromagnetic structure factors as a function of temperature for
various values of $\Delta$. Solid lines are 6th order NLC, dotted lines are
5th order NLC, dashed lines are ED results for a 10-site cluster and dash-
dotted lines are ED results for an 8-site cluster.

In this section we study properties of the model more quantitatively at high
temperatures using Numerical Linked Cluster Expansions (NLC) rigol  and Exact
Diagonalization (ED) of an 8 and a 10 site periodic cluster. In NLC, extensive
properties on a lattice ${\cal L}$ with $N$ sites are calculated, in the
thermodynamic limit $N\to\infty$, as a sum over all linked clusters $c$ as

$P({\cal L})/N=\sum_{c}L(c)\times W(c)$ (7)

Here $L(c)$, called the lattice constant, is the number of times the cluster
arises in the lattice, per lattice site. The quantity $W(c)$ is the weight of
the cluster. It is defined recursively through the property on the finite
cluster $c$ as

$W(c)=P(c)-\sum_{s}W(s),$ (8)

where the sum runs over all proper sub-clusters $s$ of the cluster $c$.

Figure 3: Heat capacity as a function of temperature for different values of
$\Delta$. Solid lines correspond to 6th order NLC and dotted lines to 5th
order NLC. Figure 4: Entropy as a function of temperature for different values
of $\Delta$. Solid lines correspond to 6th order NLC and dotted lines to 5th
order NLC.

We calculate the heat capacity per site $C$, the entropy per site $S$, the
uniform susceptibility per site $\chi$ and the antiferromagnetic structure
factor $S_{\pi,\pi}$. Calculations are done to 6th order, that is weights of
all linked clusters with 6 or less bonds are included in our study. NLC
converges at high temperatures but will break down when correlation length
exceeds a few lattice constants rigol .

The parameter regimes of the model can be divided into 3 regions as seen from
the behavior of the structure factor in Fig. 2. The first region is
$-\infty<\Delta<2$. Representative plots for the structure factors are shown
for $\Delta=0$ and $\Delta=1.5$. As $\Delta\to-\infty$ this model reduces to
spin-one Heisenberg model which is known to have a renormalized classical (RC)
phase at low temperatures. All the way up to $\Delta=2.0$ the fundamental
behavior remains unchanged. The only change in properties is that at high
temperature the existence of spin-zero states pushes down the onset of
antiferromagnetic correlations. Note that NLC and ED agree above $T/J=3$. At
lower temperature the finite size data must saturate, they are limited by
total number of spins. NLC shows a rapid rise in the structure factor. At the
lowest temperatures NLC must also break down as the correlation length
increases exponentially.

The second region is $\Delta>2.5$. Representative plots for the structure
factors are shown for $\Delta=3.0$ and $5.0$. This is a spin-zero phase in
which NLC converges well at all temperatures. The only finite size effects in
ED are at intermediate temperatures as seen at $\Delta=3.0$ and they are
relatively small.

The third region is the transition region $2.0<\Delta<2.5$. This is a
difficult region for numerical study at intermediate temperatures. As
illustration, structure factor for $\Delta=2.3$ is shown. One can see that ED
and NLC diverge below a temperature of $1.5J$. Finite sizes are too small to
capture the intermediate temperature behavior. Although NLC also breaks down
as we get close to the transition, it does give us a glimpse of the transition
behavior, which we now discuss.

In Fig. 3 we show the heat capacity as a function of temperature for $\Delta=$
$2.3$, $2.4$, $2.5$ and $3.0$. As $\Delta$ is reduced from large values, the
peak in the specific heat moves to lower temperature and it increases in
magnitude. However, at $\Delta=2.3$ the specific heat curve develops a two-
peak structure. The weak higher T peak signifies a development of spin-zero
dominance over spin-one (gas phase in the liquid-gas analogy). This gives way,
through a first order transition, to the antiferromagnetic phase (liquid or
condensed phase in the liquid-gas analogy) at still lower temperature. The
corresponding entropy curves are shown in Fig. 4. At $\Delta=2.3$, there
develops a shoulder in the entropy function followed by a sharp downturn. NLC
cannot capture the first order transition well, we expect the plateau-like
behavior to be followed by a jump in entropy at the transition.

In Fig. 5, we show the uniform susceptibility as a function of temperature.
For $\Delta=2.3$, there is a sudden rise in the susceptibility as the
transition to high-spin states occurs. This sudden rise continues even in the
RC phase as $\Delta$ is lowered further. This result is potentially relevant
to experiments on BNOAS BNOAS-expt  as we discuss in the next section. The NLC
results are not valid on the ordered side, the RC phase must have a finite
susceptibility at $T=0$ CHN .

Figure 5: Magnetic susceptibility as a function of temperature for various
values of $\Delta$. Solid lines are 6th order NLC and dotted lines are 5th
order NLC.

The above properties clearly distinguish our model from the Kondo-necklace,
where a zero temperature quantum critical point turns into a quantum-critical
fan at finite temperatures brenig ; CHN .

##  V V. Summary and Discussions

In this paper we have studied the magnetic and thermodynamic properties of a
two-electron two-orbital model system with competing high spin (S=1) and low
spin (S=0) states, possibly relevant to the nickelates with a $d^{8}$
electronic configuration. The model has some similarities to the Kondo-
necklace model proposed earlier by Jin et al jin . However, detailed
experimental properties are quite different and should be easy to distinguish
experimentally. It is the purpose of the paper to bring out these differences.

Being a two-dimensional model, there cannot be a Neel phase at finite
temperatures. Instead, there is a Renormalized Classical phase for the
antiferromagnet, which will develop Neel order at $T=0$, or at finite
temperatures in presence of the slightest anisotropy or 3-dimensional
coupling. This phase is separated by a first order phase transition from a
non-magnetic spin-zero phase. The line of first order transitions terminates
at a liquid-gas critical point.

The magnetic behavior in the vicinity of the transition is quite unusual as
the magnetic susceptibility can be non-monotonic as a function of temperature
and show a sudden rise as one transitions from a low-spin to high-spin phase.
There can be hysteresis and meta-stability of the low and high spin phases. At
the phase transition there should be a jump in entropy and a latent heat.

Jiang et al sawatzky  have discussed the importance of charge transfer energy
and the tuning (and crossing) of singlet-triplet gap as a function of the
charge transfer energy. Thus, it is possible that different materials, with
small differences in ligand-environment, can be in different parts of the
phase diagram. Some maybe magnetic, some non-magnetic and some may have
parameters close to the critical value.

The experimental results for the BNOAS material by Matsumoto et al BNOAS-expt
look very intriguing. There is a Neel transition at $130$ K. Just above the
Neel transition at around $150$ K the system has a sudden increase in the
magnetic susceptibility. The susceptibility increases by almost a factor of
$6$ between $150$ and $130$ K. The author’s speculate that this may be due to
anisotropies that cause canted ferromagnetism. It is also possible that it is
caused in part by the sudden emergence of high-spin from low-spin as seen in
the present model. If true, this is direct evidence that the system is close
to a degeneracy between high and low-spin states.

This manuscript is focused on the parent insulating materials. But, it is
interesting to speculate about doping and superconductivity. From this point
of view also this model may provide a sharp distinction to the proposal of Jin
et al. Upon doping, the Kondo-necklace model should lead to the usual t-J
model, whereas this model could lead to a spin-one or type II t-J model
recently proposed by Zhang and Vishwanath ashvin , with new mechanisms for
superconductivity.

##  VI Acknowledgments

We thank Warren Pickett for many useful discussions and for getting us
interested in this problem. This work is supported in part by National Science
Foundation grant number NSF-DMR 1855111.

## References

  * (1) H. S. Jin, W. E. Pickett and K. W. Lee, Phys. Rev. Research 2, 033197 (2020).
  * (2) Q. Si, R. Yu, E. Abrahams, Nature Rev. Mater. 1, 16017 (2016).
  * (3) D. Li, K. Lee, B. Y. Wang, M. Osada, S. Crossley, H. R. Lee, Y. Cui, Y. Hikita and H. Y. Hwang, Nature 572, 624–627 (2019).
  * (4) ] K.-W. Lee and W. E. Pickett, Phys. Rev. B 70, 165109 (2004).
  * (5) A.-S. Botana and M. R. Norman, Phys. Rev. X 10, 011024 (2020).
  * (6) W.-L. Tu, E.-G. Moon, K.-W. Lee, W. E. Pickett, and H.-Y. Lee, http://arxiv.org/abs/2107.11936.
  * (7) W. Brenig, Phys. Rev. B 73, 104450 (2006).
  * (8) M. Jiang, M. Berciu and G. A. Sawatzky, Phys. Rev. Lett. 124, 207004 (2020).
  * (9) L.-H. Hu and C. Wu, Phys. Rev. Research 1, 032046(R).
  * (10) K. I. Kugel and D. I. Khomskii, Soviet Physics Uspekhi. 25: 231 (1982).
  * (11) R. R. P. Singh, Phys. Rev. B 41, 4873(R) (1990).
  * (12) Z. Weihong, J. Oitmaa, and C. J. Hamer Phys. Rev. B 43, 8321 (1991).
  * (13) C. J. Hamer, Z. Weihong, and J. Oitmaa Phys. Rev. B 50, 6877 (1994).
  * (14) M. Rigol, T. Bryant, and R. R. P. Singh, Phys. Rev. Lett. 97, 187202 (2006); M. Rigol, T. Bryant, and R. R. P. Singh Phys. Rev. E 75, 061118 (2007).
  * (15) S. Chakravarty, B. I. Halperin, and D. R. Nelson Phys. Rev. B 39, 2344 (1989).
  * (16) Y. Matsumoto, T. Yamamoto, Kousuke Nakano, H. Takatsu, T. Murakami, K. Hongo, Ryo Maezono, H. Ogino, D. Song, C. M. Brown, Cédric Tassel, and H. Kageyama, Angew Chem Int Ed Engl. 58, 756–759 (2019).
  * (17) Y.-H. Zhang and A. Vishwanath, Phys. Rev. Research 2, 023112 (2020).