# Learning to Dynamically Select Cost Optimal Schedulers in Cloud Computing
Environments

Shreshth Tuli1 Giuliano Casale1 Nicholas R. Jennings1,2
1Imperial College London  UK  2Loughborough University  UK

###### Abstract

The operational cost of a cloud computing platform is one of the most
significant Quality of Service (QoS) criteria for schedulers, crucial to keep
up with the growing computational demands. Several data-driven deep neural
network (DNN)-based schedulers have been proposed in recent years that
outperform alternative approaches by providing scalable and effective resource
management for dynamic workloads. However, state-of-the-art schedulers rely on
advanced DNNs with high computational requirements, implying high scheduling
costs. In non-stationary contexts, the most sophisticated schedulers may not
always be required, and it may be sufficient to rely on low-cost schedulers to
temporarily save operational costs. In this work, we propose MetaNet, a
surrogate model that predicts the operational costs and scheduling overheads
of a large number of DNN-based schedulers and chooses one on-the-fly to
jointly optimize job scheduling and execution costs. This facilitates
improvements in execution costs, energy usage and service level agreement
violations of up to 11%, 43% and 13% compared to the state-of-the-art methods.

##  1 INTRODUCTION

The onset of the Artificial Intelligence (AI) and Deep Learning (DL) era has
led to a recent shift in computation from hand-encoded algorithms to data-
driven solutions for resource management in cloud systems [2]. This has given
rise to efficiently harnessing the data processing capacities of multiple
devices and provides services at scale with high Quality of Service (QoS).
However, the increasing computational demands of modern Internet of Things
(IoT) applications makes it crucial to curtail the operational costs of cloud
machines. This calls for efficient resource management schemes, such as task
scheduling policies, to execute workloads on cloud systems within tight cost
budgets. AI and DL offer promising solutions that rely on accurate surrogate
models that are inexpensive to evaluate at run-time.

Background and Motivation. In recent years, state-of-the-art resource
management solutions, which particularly focus on optimal placement of tasks
on cloud virtual machines (VMs), have increasingly explored data-driven DL
methods [7, 9, 10, 11, 12]. Such methods typically rely on a trace of resource
utilization characteristics and QoS metrics of tasks and cloud hosts and are
referred to as trace-driven schedulers. They utilize such traces to train a
deep neural network (DNN) to estimate a set of Quality of Service (QoS)
parameters and run optimization strategies to find the best scheduling
decision for each incoming task. However, most prior work run inference using
such trace-driven DNN models on a broker node to generate scheduling
decisions, while the incoming tasks are executed on worker nodes [11]. In such
cases, having a dedicated broker is often cost inefficient due to the idle
times between decisions in discrete-time control settings, wherein the task
placement decisions are taken at fixed scheduling intervals [11]. To tackle
this, we resort to paradigms such as Function as a Service (FaaS) that allow
execution of DL schedulers as serverless functions, only incurring costs for
the inference run time of each DNN model.

Contributions. In this work, to leverage the recent advantages brought by DL
in task scheduling and build a policy agnostic solution, we choose a set of
state-of-the-art schedulers. This work aims to solve the meta-problem of on-
the-fly selection of scheduling policies for a cloud computing environment
wherein the incoming tasks are executed on worker nodes and scheduling
decisions are run as serverless functions. To solve this meta-problem, we
develop the proposed solution that we call MetaNet. It uses a DNN as a
surrogate model to predict the task execution costs and scheduling time for
each policy. We then select the most cost-efficient policy at each scheduling
interval using the online estimates generated by this surrogate. As we
dynamically update the policy to tradeoff between simple and sophisticated DL
schedulers, we reduce overall operational costs, energy consumption and
Service Level Agreement (SLA) deadline violation rates by up to 11%, 43% and
13% respectively compared to state-of-the-art schedulers.

##  2 BACKGROUND AND RELATED WORK

Recent work in scheduling for cloud computing environments has demonstrated
that AI-based solutions are not only faster, but can also scale efficiently
compared to traditional heuristic and optimization techniques [9, 10, 11, 12].

Evolutionary Optimization. This class of methods forecasts QoS using non-DL
solutions such as ARIMA [8] or DL based models such as Long-Short-Term-Memory
(LSTM) neural networks [6]. It then applies evolutionary search strategies
such as Ant Colony Optimization (ACO) [1], to converge to a locally optimal
scheduling decision.

Surrogate Optimization. This class of methods uses differentiable function
approximators, particularly neural networks, to act as surrogates of the QoS
for a future state of the system. For instance, GRAF [7] uses a graph neural
network (GNN) as a surrogate model to predict service latencies and
operational costs for a given scheduling decision and uses gradient-based
optimization to minimize the service latencies or execution costs. To do this,
it uses the concept of neural network inversion [13], wherein the method
evaluates gradients of the objective function with respect to inputs and runs
optimization in the input space. Other methods, such as Decision-NN [12], GOBI
[11] and its second-order generalization GOSH [10], combine the prediction and
optimization steps by modifying the loss function to train and optimize in
tandem.

##  3 PROPOSED METHOD

Figure 1: MetaNet System Model.

System Model. In this work, we target a standard heterogeneous cloud computing
environment where all nodes are in the same Wide Area Network (see Figure 1
for an overview). Tasks are container instances that need to be processed,
generated from an IoT layer and transferred to the compute nodes via gateway
devices. We do not have any cloud broker; instead, we run the scheduling
policies as FaaS functions on a serverless platform such as AWS Lambda or
Azure Serverless. As is common in prior work [1, 9, 12], we consider a
discrete-time control problem, i.e., we divide the timeline into fixed-size
execution intervals (of $\Delta$ seconds). We denote the overall cost of the
system amortized by the number of completed tasks in the $t$-th interval as
$\phi_{t}$. This includes the operational costs of the worker nodes as well as
that of running the schedulers as serverless functions.

MetaNet Methodology. We create a surrogate model $f_{\theta}$ that is a DNN
with parameters $\theta$. We consider a set of scheduling policies
$\mathcal{P}$ of size $q$. Given a system state at the start of the $t$-th
interval, each scheduler produces a scheduling decision for this interval. The
state is denoted by $C_{t-1}$ and includes the resource (CPU, RAM and disk)
utilization characteristics of the active workloads in the system $W_{t-1}$,
resource utilization of cloud hosts $H_{t-1}$ and scheduling decision of the
previous interval $S_{t-1}$. Thus, $C_{t-1}=[W_{t-1},H_{t-1},S_{t-1}]$. Now,
$f_{\theta}$ estimates the average execution cost (sum of the cost of task
execution and serverless run) for each scheduler $p^{k}\in\mathcal{P}$ in the
$t$-the interval as $\hat{\phi}^{k}_{t}$. The collection of
$\hat{\phi}^{k}_{t}$ for all $k$ in $I_{t}$ is denoted by $\hat{\phi}_{t}$. In
such a case, our problem can be formulated as

$\displaystyle\underset{\theta}{\text{minimize}}$
$\displaystyle\sum_{t=1}^{T}\phi^{\pi}_{t}$ (1) subject to $\displaystyle
S_{t}=p^{k}(t),\forall\ t$
$\displaystyle\pi=\operatorname*{arg\,min}_{k}\hat{\phi}^{k}_{t},\forall\ t$
$\displaystyle\hat{\phi}_{t}=f_{\theta}(W_{t-1},H_{t-1},S_{t-1}),\forall\ t,$

where $f_{\theta}$ is a surrogate of the average execution cost.

Figure 2: MetaNet Pipeline.

Surrogate Model and Training. We use a feed-forward neural model that takes as
input the state of the cloud system $C_{t-1}=[W_{t-1},H_{t-1},S_{t-1}]$ and
outputs a single scalar. We use a fully connected neural network with 3
layers, each with 64 nodes and $\mathrm{ReLU}$ activation, except the last
layer where we use the $\mathrm{Sigmoid}$ activation to bring the output in
the range $(0,1)$. Thus, for any input $(W_{t-1},H_{t-1},S_{t-1})$, the
complete neural network can be described as

$\hat{\phi}=f_{\theta}(W_{t-1},H_{t-1},S_{t-1}).$ (2)

To train the above described MetaNet neural network $f_{\theta}$, we collect
traces from a cloud computing environment. To collect data for training, we
execute the scheduling policies $\mathcal{P}$, each for $\Gamma$ scheduling
intervals, and collect a dataset as

$\Lambda=\\{(k,W_{t-1},H_{t-1},S_{t-1},\phi^{k}_{t},\omega^{k}_{t})\\}_{t=1}^{\Gamma\cdot
q},$ (3)

where $q$ is the number of scheduling policies in $\mathcal{P}$. We initialize
$W_{0}$ as $H_{0}$ zero-matrices and $S_{0}$ as an empty graph. We find the
maximum cost from the dataset as $\phi^{k}_{max}$ for each scheduler $p^{k}$.
This allows us to denormalize the neural network output and bring it to the
same range as the one in the dataset. We then train the model using the loss
function

$L=\lVert\phi^{k}-\phi^{k}_{max}\cdot\hat{\phi}^{k}\rVert_{2}.$ (4)

Dynamic Policy Selection Using MetaNet. We first store the saved surrogate
model $f_{\theta}$ into a central network-attached-storage (NAS) that is
accessible to all worker nodes. At start the scheduling interval $I_{t}$, we
get the workload and host characteristics $W_{t-1},H_{t-1}$ with the
scheduling decision $S_{t-1}$. As we do not have any broker node in our setup,
we run the MetaNet model in the worker node with the least CPU utilization. On
this worker node we use the surrogate model to predict cost and decide the
scheduling policy as

$p^{\pi},\text{ s.t. }\pi=\operatorname*{arg\,min}_{k}\hat{\phi}^{k}_{t}.$ (5)

Overall, MetaNet selects a scheduling policy on-the-fly as per the system
states. To do this, at the start of each scheduling interval, it predicts the
task execution cost and scheduling time of each policy. It then selects the
one with the minimum cost estimate. This policy optimizes the scheduling
decision, initialized as $S_{t-1}$, to generate $S_{t}$. The complete pipeline
is summarized in Figure 2. We also tune $f_{\theta}$ with each new datapoint
to adapt to non-stationary settings.

##  4 PERFORMANCE EVALUATION

Setup. We consider the complete set of hosts $\mathcal{H}$ to be static with
time as is common in prior work for a fixed cloud platform [9]. We use
different VM types from the Microsoft Azure Platform, i.e., 60 Azure B2s with
a dual-core CPU and 4GB RAM (in UK-South), 20 Azure B4ms with a quad-core CPU
and 16GB RAM and 20 Azure B8ms with an octa-core CPU and 32 GB RAM (in East-
US). To save on costs, we define a host to be active when the CPU utilization
is $>0\%$. We utilize the Azure Automation
(https://azure.microsoft.com/services/automation/#overview) service to
hibernate and resume VMs based on their CPU utilization.

Figure 3: Predicted costs (line plots) and dynamic scheduler selection
(background color) in MetaNet.

Workloads. To generate the tasks in our system, we use the AIoTBench
applications [5]. AIoTBench is an AI based cloud computing benchmark suite
that consists of various real-world computer vision application instances. The
seven specific application types correspond to the neural networks they
utilize. These include three typical heavy-weight networks: ResNet18,
ResNet34, ResNext32x4d, as well as four light-weight networks: SqueezeNet,
GoogleNet, MobileNetV2, MnasNet. At the start of each scheduling interval, we
create new tasks, sampled uniformly from the seven applications, where the
number of new tasks comes from a Poisson distribution with rate $\lambda=1.2$.

Baselines. We compare MetaNet against seven baselines, which also form our
policy set $\mathcal{P}$. We integrate the ACO algorithm with two demand
forecasting methods: AutoARIMA and LSTM, and call these ARIMA+ACO and
LSTM+ACO. We also include Decision-NN, Semi-Direct and GRAF, GOBI and GOSH.
This makes the size of the policy set $q=7$. We also use a Multi-Armed Bandit
(MAB) model using Upper-Confidence-Bound exploration [3] and a Deep-Q-Network
(DQN) [4] that choose a policy based on pre-trained models using data
$\Lambda$.

Implementation and Training. We use implement MetaNet on the COSCO framework
[11]. We collect the dataset $\Lambda$ by executing all baselines for
$\Gamma=100$ intervals. Similarly, we execute all approaches for $T=1000$
scheduling intervals to generate QoS scores, with each interval being
$\Delta=$ 10 seconds long, giving a total experiment time of nearly 2 hours 46
minutes for each method.

Visualization of MetaNet. Figure 3 visualizes the MetaNet approach running on
the setup and workloads described above. The x-axis denotes the scheduling
interval and the y-axis denotes the cost estimate $\hat{\phi}^{k}_{t}$ for the
three most frequently selected policies in $\mathcal{P}$ for readability:
ARIMA+ACO, GOBI and GOSH. The highlighted bands indicate the selected
scheduling policy. The selected policy corresponds to the one with the least
estimated cost. The cost estimates are non-stationary, further corroborating
the need for dynamic selection of the scheduling policies in volatile workload
and host setups.

Comparison with Baselines. Table 1 shows the results of our experiments (cost
in USD, energy in KW-hr, resp. time in seconds). Across all metrics, MetaNet
outperforms the baselines. MetaNet improves the energy consumption by
allocating tasks, i.e., to the same hosts to minimize execution costs and
consequently the active hosts in the system. This is shown by the highest CPU
utilization of MetaNet, i.e., 87.4%. MetaNet also gives the lowest response
time and consequently the lowest SLA violation rates. Dynamic optimization
baselines perform poorly due to the stateless assumption in MAB that ignores
environment dynamism, and DQN being slow to adapt in volatile settings [11].

Table 1: Comparison with the baseline.

Model | Cost | Energy | Resp. T. | SLA V. | CPU %
---|---|---|---|---|---
ARIMA+ACO | 0.794 | 5.432 | 20.391 | 0.241 | 62.4
LSTM+ACO | 0.862 | 5.217 | 16.921 | 0.212 | 76.5
Decision-NN | 1.072 | 4.224 | 12.255 | 0.211 | 80.2
Semi-Direct | 1.021 | 4.095 | 17.092 | 0.131 | 72.2
GRAF | 0.993 | 4.228 | 14.722 | 0.127 | 76.1
GOBI | 0.752 | 3.827 | 14.293 | 0.122 | 80.2
GOSH | 0.911 | 3.267 | 13.292 | 0.118 | 84.3
MAB | 0.874 | 2.921 | 13.921 | 0.121 | 79.9
DQN | 0.907 | 3.121 | 13.877 | 0.119 | 80.4
MetaNet | 0.667 | 2.029 | 11.273 | 0.102 | 87.4

##  5 Conclusions

This work presents MetaNet, which dynamically selects scheduling policies for
cost-efficient task processing in cloud systems. Future work would explore
other DNNs as surrogates and extend MetaNet to other types of resource
management decisions such as dynamic VM provisioning and autoscaling.

## References

  * [1] Muhammad Aliyu, M Murali, Abdulsalam Y Gital and Souley Boukari  “Efficient metaheuristic population-based and deterministic algorithm for resource provisioning using ant colony optimization and spanning tree”  In _International Journal of Cloud Applications and Computing (IJCAC)_ 10.2 IGI Global, 2020, pp. 1–21
  * [2] Hongming Cai, Boyi Xu, Lihong Jiang and Athanasios V Vasilakos  “IoT-based big data storage systems in cloud computing: perspectives and challenges”  In _IEEE Internet of Things Journal_ 4.1 IEEE, 2016
  * [3] Volodymyr Kuleshov and Doina Precup  “Algorithms for multi-armed bandit problems”  In _arXiv:1402.6028_ , 2014
  * [4] Yuxi Li  “Deep reinforcement learning: An overview”  In _arXiv:1701.07274_ , 2017
  * [5] Chunjie Luo et al.  “AIoT bench: towards comprehensive benchmarking mobile and embedded device intelligence”  In _International Symposium on Benchmarking, Measuring and Optimization_ , 2018, pp. 31–35  Springer
  * [6] Soukaina Ouhame, Youssef Hadi and Arif Ullah  “An efficient forecasting approach for resource utilization in cloud data center using CNN-LSTM model”  In _Neural Computing and Applications_ Springer, 2021, pp. 1–13
  * [7] Jinwoo Park, Byungkwon Choi, Chunghan Lee and Dongsu Han  “GRAF: a graph neural network based proactive resource allocation framework for SLO-oriented microservices”  In _Proceedings of the 17th International Conference on emerging Networking EXperiments and Technologies_ , 2021, pp. 154–167
  * [8] Parminder Singh, Pooja Gupta and Kiran Jyoti  “TASM: Technocrat ARIMA and SVR model for workload prediction of web applications in cloud”  In _Cluster Computing_ 22.2 Springer, 2019, pp. 619–633
  * [9] Peter J Stuckey et al.  “Dynamic programming for predict+optimise”  In _AAAI_ 34.02, 2020
  * [10] Shreshth Tuli, Giuliano Casale and Nicholas R. Jennings  “GOSH: Task Scheduling using Deep Surrogate Models in Fog Computing Environments”  In _IEEE Transactions on Parallel and Distributed Systems_ , 2022
  * [11] Shreshth Tuli et al.  “COSCO: Container Orchestration Using Co-Simulation and Gradient Based Optimization for Fog Computing Environments”  In _IEEE Transactions on Parallel and Distributed Systems_ 33.1, 2021, pp. 101–116
  * [12] Bryan Wilder, Bistra Dilkina and Milind Tambe  “Melding the data-decisions pipeline: Decision-focused learning for combinatorial optimization”  In _AAAI_ 33.01, 2019
  * [13] Eric Wong and J Zico Kolter  “Neural network inversion beyond gradient descent”  In _Advances in Neural Information Processing Systems, Workshop on Optimization for Machine Learning_ , 2017


# Numerical and Physical Stability of Supernova Neutrino Flavor Evolution

B. D. Keister  Physics Division, National Science Foundation, Arlington, VA
22230

###### Abstract

This paper examines neutrino flavor evolution outside a supernova
neutrinosphere using a one-dimensional model that retains the non-linear
nature of neutrino-neutrino interactions as well as some aspects of the full
geometry. In some limiting cases analytic results can be obtained that display
different behavior from their counterparts in (linear) solar neutrino flavor
evolution. For more general cases, numerical solutions require extended
numerical algorithms to achieve stable solutions, and these solutions exhibit
standard chaotic behavior.

##  I Introduction

Neutrinos are trapped by multiple scattering from baryons at nuclear matter
densities in a supernova. However, outside the neutrinosphere, which
characterizes the boundary of the baryon-rich region, neutrinos can propagate
outward as coherent waves. There is an extensive literature aimed at
describing neutrino transport outside supernovae where the scattering medium
can include the neutrinos themselves as well as electrons and baryons Duan et
al. (2010); Pantaleone (1992); Notzold and Raffelt (1988); Samuel (1993,
1996); Pastor and Raffelt (2002); Hannestad et al. (2006); Dasgupta et al.
(2009); Pastor et al. (2002); Balantekin and Pehlivan (2007); Dasgupta et al.
(2012); Friedland (2010); Duan et al. (2011); Mirizzi and Serpico (2012). [For
a review, see Ref. Duan et al. (2010) and references therein.]

Recently concerns have been raised regarding numerical stability and/or
sensitivy of the numerical results to geometric symmetry Raffelt and Seixas
(2013); Raffelt et al. (2013); Duan and Shalgar (2014). The purpose of this
paper is to gain some insight into the results of three-dimensional
simulations by means of a one-dimensional model. The model retains the non-
linear nature of neutrino-neutrino scattering, as well as some of the density
variation of particle species found in the full three-dimensional picture. It
also permits a study of numerical algorithms beyond those typically used in
the three-dimensional case. In our analysis, we find that achieving numerical
stability requires a higher order algorithm than what is normally employed in
the three-dimensional case. Once numerical stability is achieved, the
solutions for all but the most symmetric cases exhibit the the standard
criteria for chaotic behavior.

The model also provides a second insight with respect to the asymptotic flavor
evolution of neutrinos far from a supernova where the neutrino density has
fallen to zero. It address the question whether neutrino flavors evolve in the
same manner as they do when traveling in a static medium electrons outside the
sun Wolfenstein (1978); Mikheev and Smirnov (1985). It will be seen that the
non-linear nature of supernova neutrino evolution leads to different
configurations at large distances than one might expect from the solar
neutrino case, which is linear in form.

##  II Background: Neutrino Medium Refractive Index

A refractive index describes coherent propagation of a particle in a medium of
scatterers, for which plane waves travel with an effective wavenumber:

$\kappa=kn,$ (1)

where the index $n$ depends upon the forward scattering amplitude of the
scatterer:

$n=1+{2\pi\rho\over k^{2}}f(0),$ (2)

and $\rho$ is the number density of the medium particles. van de Hulst (1949)

This leads to an effective wave equation: for any direction fixed by a
wavenumber ${\mathbf{k}}$, the spatial variation apart from the overall phase
$e^{ikx}$ is given by

${\partial\over\partial x}\psi=iK\psi,$ (3)

where

$K={2\pi\rho\over k}f(0).$ (4)

$K$ thus plays the role of a Hamiltonian.

For neutrinos, the index becomes a matrix in flavor space. The matrix
contributions are diagonal for electron and baryon scatterers. In the absense
of mixing, the neutrino medium also contributes a diagonal matrix index. For
simplicity, we consider only two neutrino flavors $e$ and $x$ (e.g. assuming
flavors $\nu$ and $\tau$ to be fully mixed). In the presence of neutrino
mixing, the matrix index contribution from each neutrino scatterer is diagonal
a basis consisting of its mixed state plus its orthogonal partner. The index
then becomes a density matrix Pantaleone (1992).

The mean-free path,

$L\sim{1\over\rho\sigma_{\rm tot}}.$ (5)

provides a figure of merit by which to judge a calculation using the index
alone: additional non-coherent contributions arising from inelastic scattering
or elastic scattering far from forward-going directions start to be come
significant only over distances comparable to $L$.

From Ref. Duan et al. (2006), we note that the maximum neutrino density
outside the neutrinosphere is roughly $10^{32}$ cm-3. The total elastic
$\nu-\nu$ cross section (without exchange) is

$\sigma={G^{2}s\over 4\pi}.$ (6)

For neutrino energies of 10 MeV, this corresponds to a mean-free path of order
50,000 km. This means that a refractive index picture in which neutrinos
travel through a medium of other neutrinos is an excellent approximation for
the length scales of interest, typically hundreds of kilometers or less. This
observation is important for what follows: however difficult it may be to
implement the equations of motion, they capture the dominant physics. Put
another way, under these density conditions, corrections to the refractive
limit are very unlikely to alter the conclusions.

In the presence of an additional medium of electrons, the complete expression
is

$K=K_{\nu}+K_{e}+K_{m},$ (7)

where

$K_{e}=\sqrt{2}G\rho_{e}\left[\begin{array}[]{cc}1&0\\\
0&0\end{array}\right],$ (8)

and

$K_{m}={\Delta m^{2}\over 2E}\frac{1}{2}\left[\begin{array}[]{cc}-\cos
2\theta&\sin 2\theta\\\ \sin 2\theta&\cos 2\theta\end{array}\right].$ (9)

The coupled neutrino term receives a contribution from each scatterer with
momentum ${\mathbf{p}}_{2}$:

$K_{\nu}=2\sqrt{2}G\int
d{\mathbf{p}}_{2}\,(1-{\mathbf{p}}_{1}\cdot{\mathbf{p}}_{2})\left[\begin{array}[]{cc}\rho_{\nu_{e}}({\mathbf{p}}_{2})&0\\\
0&\rho_{\nu_{x}}({\mathbf{p}}_{2})\end{array}\right],$ (10)

where ${\mathbf{p}}_{1}$ is the neutrino beam momentum. In full calculations,
the integral over ${\mathbf{p}}_{2}$ is weighted with a distribution function.

##  III Model Studies

The focus of our study is a model with a single space-time variable $x\sim t$.

${\partial\over\partial x}\psi(x)=iK\psi(x).$ (11)

We consider the evolution of two groups of neutrinos, labeled $A$ and $B$,
respectively, and and associated wave functions $\psi$ an $\chi$, evolving
under a single variable $x$:

$\displaystyle{\partial\over\partial x}\psi$ $\displaystyle=$ $\displaystyle
K^{A}\psi$ (12) $\displaystyle{\partial\over\partial x}\chi$ $\displaystyle=$
$\displaystyle K^{B}\chi,$

where $K^{A}$ and $K^{B}$ are identical except for the nonlinear neutrino-
neutrino interaction term:

$\displaystyle K^{A}_{\nu}$ $\displaystyle=$ $\displaystyle
G\rho_{\nu}\,\chi\otimes\chi^{\dagger};$ (13) $\displaystyle K^{B}_{\nu}$
$\displaystyle=$ $\displaystyle G\rho_{\nu}\,\psi\otimes\psi^{\dagger}.$

This model does not include the full effects of three-dimensional geometry,
but retains the non-linear nature of the optics.

Full-scale simulations employ a thermal distribution of neutrinos. In fact, a
frequently discussed outcome of these simulations is the idea an energy swap.
For this model, we retain a minimal aspect of a distribution by following the
evolution of two incoherent neutrinos in each group:
$\psi_{<},\psi_{>},\chi_{<},\chi_{>}$, corresponding to two distinct neutrino
energies $E_{<}$ and $E_{>}$.

At $x=0$, each of the four neutrinos begins propagation as a flavor
eigenstate, e.g.

$\psi_{e}(x=0)=\left[\begin{array}[]{c}1\\\ 0\end{array}\right]$ (14)

for pure $\nu_{e}$, and

$\psi_{x}(x=0)=\left[\begin{array}[]{c}0\\\ 1\end{array}\right]$ (15)

for pure $\nu_{x}$. We take $K_{m}$ as given in Eq. 9. For $K_{e}$, we write

$K_{e}=\sqrt{2}G\rho_{e}F_{e}(x){1\over 2}\left[\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right].$ (16)

In this study we have removed all matrix traces as they do not contribute to
flavor oscillations.

The form factor $F_{e}(x)$ accounts for the density falloff of electron
density in three-dimensional space:

$F_{e}(x)=({R\over R+x})^{2}$ (17)

where $R$ sets the scale of the neutrinosphere radius.

$K_{\nu}=\sqrt{2}G\rho_{e}F_{\nu}(x)\,\chi\otimes\chi^{\dagger},$ (18)

where there is an implied incoherent sum over the two energies $E_{<}$ and
$E_{>}$. The form factor $F_{\nu}(x)$ accounts both for the spatial falloff of
neutrino density (as for electrons) and the fact that a typical angle of
incidence of neutrino pairs in the full three-dimensional problem approaches
zero (see Eq. 19):

$(1-{\mathbf{p}}_{1}\cdot{\mathbf{p}}_{2})\sim O\left({1\over r^{2}}\right).$
(19)

We then have

$F_{\nu}(x)=({R\over R+x})^{4}$ (20)

The model has corresponding features to that reported in Ref. Raffelt and
Seixas (2013). That model also has a single variable, which can be considered
to be one spatial dimension. It has no electron or neutrino density variation,
unlike the framework described above. It also has fixed rather than variable
neutrino densities. Nevertheless, some conclusions are quite similar.

###  III.1 Analytic Results

We examine here the possibility of a flavor evolution effect due to neutrino-
neutrino interactions that is analagous to the MSW effect due to the neutrino-
electron interaction. In the MSW effect, solar electron neutrinos propagate
initially through a medium of electrons whose density is high enough that its
optical effect dominates over vacuum oscillations. As the electron density
decreases adiabatically with distance, the neutrino state gradually evolves
into a mass-2 eigenstate via the non-crossing rule. For supernova neutrinos,
one could imagine that, for a neutrino population restricted to electron
neutrinos only, these neutrinos would propagate in a medium dominated by
electron flavor; as the electron neutrino density falls off, these neutrinos
would evolve to mass-2 eigenstates. In fact, this does not happen: two pure
flavor eigenstates evolve as they pass through each other to a final state
containing the same mass eigenstate mixtures as their initial states.

We can demonstrate this by induction. In the mass basis, a pure $\nu_{e}$
initial state is

$\psi_{\nu_{e}}(x=0)=\left[\begin{array}[]{c}\cos\theta\\\
\sin\theta\end{array}\right].$ (21)

According to the claim, this state at finite $x$ has the form

$\psi_{\nu_{e}}(x)=\left[\begin{array}[]{c}e^{i\alpha}\cos\theta\\\
e^{i\beta}\sin\theta\end{array}\right].$ (22)

We now demonstrate that this form is preserved under a step $x\to x+h$.

$d\psi=ih(K_{\nu}+K_{m})\psi.$ (23)

In the mass basis,

$K_{m}\psi=\lambda\left[\begin{array}[]{cc}-1&0\\\
0&+1\end{array}\right]\left[\begin{array}[]{c}e^{i\alpha}\cos\theta\\\
e^{i\beta}\sin\theta\end{array}\right]=\lambda\left[\begin{array}[]{c}-e^{i\alpha}\cos\theta\\\
+e^{i\beta}\sin\theta\end{array}\right].$ (24)

The neutrino-neutrino matrix is proportional to $\chi\chi^{\dagger}$:

$K_{\nu}=\mu\chi\chi^{\dagger}.$ (25)

For this case, the symmetry requires that $\psi$ and $\chi$ evolve in the same
way. This means that the density matrix $K_{\nu}$ simply identifies $\psi$
with eigenvalue +1:

$K_{\nu}\psi=+\psi.$ (26)

Hence we get

$d\psi=ih\left[\begin{array}[]{c}(\mu-\lambda)e^{i\alpha}\cos\theta\\\
(\mu+\lambda)e^{i\beta}\sin\theta\end{array}\right].$ (27)

That is, the mass admixtures do not change in magnitude; only the relative
phases change. Unlike the case of MSW, the pure electron state never evolves
into a pure mass-2 eigenstate, even though the eigenvalues of the propagation
matrix never cross at any point in their evolution. The variation of $K$ is
adiabatic with respect to stellar dimensions; however, the non-linear term in
the propagation matrix varies at exactly the same rate as the wave function on
which it operates, so in that sense adiabaticity is lost. Since the neutrino
mass states eventually decouple, the asymptotic survival probability for a
state originally produced as $\nu_{e}$ is

$P_{e}\mathop{\longrightarrow}_{x\to\infty}(\cos^{4}\theta+\sin^{4}\theta)=1-2\cos^{2}\theta\sin^{2}\theta.$
(28)

The conclusions of this analysis can be readily extended for any case where
$\chi_{<}=\psi_{<};\,\chi_{>}=\psi_{>}$ at $x=0$.

###  III.2 Numerical Studies

For all but the most symmetric cases Eqs. 12 and 13 must be integrated
numerically. We use the following parameters:

$\Delta m^{2}=.003~{}{\rm eV}^{2};\,\theta=0.1;\,\rho_{\nu}=10^{32}~{}{\rm
cm}^{-3};\,\rho_{e}=10^{30}~{}{\rm cm}^{-3};\,R=10~{}{\rm km}.$ (29)

####  III.2.1 Integration Algorithms

The nonlinear dependence of $K_{\nu}$ upon neutrino wave functions means that
one must take extra care when integrating the equation of motion numerically,
so as to avoid serious loss of precision and/or loss of unitarity
(normalization of $\psi$ and $\chi$). The method used in Ref. Duan et al.
(2006) uses what we denote unitarized Euler (UE). For each step increment $h$,
the equation of motion has the form

$\Delta\psi=iKh\psi,$ (30)

where $\psi$ represents one of the two-component flavor spinors. If $K$ is
traceless, one can rewrite Eq. 30 as

$\Delta\psi=i\Theta{\bm{\sigma}}\cdot{\mathbf{n}}\psi,$ (31)

To the same order in $h$, we can also write (cf. Eq. 42 of Ref. Duan et al.
(2006))

$\psi+\Delta\psi=e^{i\Theta{\bm{\sigma}}\cdot{\mathbf{n}}}\psi=(\cos\Theta+i\sin\Theta{\bm{\sigma}}\cdot{\mathbf{n}})\psi,$
(32)

where

$\Theta=h\sqrt{|K_{ee}|};\,{\bm{\sigma}}\cdot{\mathbf{n}}={h\over\Theta}K.$
(33)

Equations 32 and 33 preserve unitarity at every step. The overall numerical
uncertainty is linear in $h$.

The UE method is very sensitive to small, likely unphysical, perturbations
near $x=0$. Accordingly, the authors of Ref. Duan et al. (2006) added an extra
contribution to $K_{e}$ that is very high near the neutrinosphere but falls
off exponentially. We include such a term for all numerical calculations:

$K_{e}^{\prime}=10^{6}\times K_{e}(x=0)\times e^{-50x/R}.$ (34)

We found for this model that we still could not always maintain numerical
stability, and therefore sought an alternative. One method that achieves a
higher level of precision while maintaining unitarity is fourth-order implicit
Gauss-Runge-Kutta (GRK) Dieci et al. (1994); Butcher (1964).

####  III.2.2 Convergence

The GRK method has a global convergence rate of $O(h^{4})$ for step size $h$.
To test convergence to $O(h^{4})$, we write the partial integration of one
component of the coupled, linear differential equation at a particular point
$x$ as

$F_{h}(x)=F_{e}(x)+Ah^{4},$ (35)

where $F_{e}(x)$ is the exact result and $Ah^{4}$ is the leading term in the
expression for the error. If we calculate $F(x)$ using various step sizes
$h_{i}$, we can form ratios involving different $h_{i}$:

$\frac{F_{h_{i}}(x)-F_{h_{j}}(x)}{F_{h_{j}}(x)-F_{h_{k}}(x)}=\frac{h_{i}^{4}-h_{j}^{4}}{h_{j}^{4}-h_{k}^{4}}.$
(36)

We can then take a ratio of ratios:

$\left(\frac{F_{h_{i}}(x)-F_{h_{j}}(x)}{F_{h_{j}}(x)-F_{h_{k}}(x)}\right)/\left(\frac{h_{i}^{4}-h_{j}^{4}}{h_{j}^{4}-h_{k}^{4}}\right)=1.$
(37)

For cases where convergence could be achieved, this ratio was tested for a
variety of numerical runs using a sequence of $h$ values, typically in the
1-cm range, at points where the result showed significant fluctuation. In each
case tested, the ratio in Eq. 37 was within one or two percent of unity. A
convergence at another power of $h$ (e.g. $h^{2}$) would not yield such a
ratio for this range of step sizes. In order to achieve this power-law
convergence, some runs required up to 25 iterations to compute the implicit
GRK integration coefficients at each step, together with an iteration
precision of one part in $10^{15}$.

The UE method should have a convergence rate of $O(h)$. A similar ratio can be
constructed that is linear in $h$ in order to test this property.. For runs
and regions where the results appear to be stable, the ratio is indeed close
to unity.

####  III.2.3 Test Runs

We discuss here a representative set of test runs for initial conditions that
have varying symmetry between groups $A$ and $B$.

Figure 1: Mass-1 emergence probability associated with each of four neutrinos
as a function of $x$, with parameters of Case 1. Left panel: GRK method; right
panel: UE method

Figure 2: Mass-1 emergence probability associated with each of four neutrinos
as a function of $x$, with parameters of Case 2. Left panel: GRK method; right
panel: UE method. In this particular case, $B_{<}=A_{<};\,B_{>}=A_{>}$.

Figure 3: Mass-1 emergence probability associated with each of four neutrinos
as a function of $x$, with parameters of Case 3. Left panel: GRK method; right
panel: UE method

Figure 4: Mass-1 emergence probability associated with each of four neutrinos
as a function of $x$, with parameters of Case 4. Left panel: GRK method; right
panel: UE method

  * •

Case 1: For our initial case, we take both energies to be the same:

$E_{<}=E_{>}=10~{}{\rm MeV}$ (38)

and

$|\psi_{<}^{e}|^{2}=0.5;\,|\psi_{<}^{x}|^{2}=0;\,|\psi_{>}^{e}|^{2}=0;\,|\psi_{<}^{x}|^{2}=0.5;\,|\chi_{<}^{e}|^{2}=0.45;\,|\chi_{<}^{x}|^{2}=0;\,|\chi_{>}^{e}|^{2}=0;\,|\chi_{<}^{x}|^{2}=0.55.$
(39)

That is, equal populations of $\nu_{e}$ and $\nu_{x}$ in one direction pass
through slightly differing populations of $\nu_{e}$ and $\nu_{x}$ in the other
direction. The results are shown in Fig. 1. In each case we plot the emergence
probability of mass-1 neutrinos (choosing mass eigenenstates since those
propagate over very large distances without further mixing). In this case,
both the GRK and UE integration methods converge, and yield the same results
to $O(h^{4})$ and $O(h)$, respectively.

  * •

Case 2: In this case, we use a two-energy spectrum:

$E_{<}=9~{}{\rm MeV};\,E_{>}=11~{}{\rm MeV}$ (40)

and equal neutrino populations $A$ and $B$:

$|\psi_{<}^{e}|^{2}=0.45;\,|\psi_{<}^{x}|^{2}=0;\,|\psi_{>}^{e}|^{2}=0;\,|\psi_{<}^{x}|^{2}=0.55;\,|\chi_{<}^{e}|^{2}=0.45;\,|\chi_{<}^{x}|^{2}=0;\,|\chi_{>}^{e}|^{2}=0;\,|\chi_{<}^{x}|^{2}=0.55.$
(41)

The symmetry of the initial conditions means that both $A$ and $B$ neutrino
populations will evolve in the same way. The results are shown in Fig. 2 In
this case, the GRK method converges to $O(h^{4})$. The UE method does not
converge at any value of $x$ after $x~{}1$ km.

  * •

Case 3: In this case, we gain use equal energies

$E_{<}=E_{>}=10~{}{\rm MeV}$ (42)

and a different population from that of Case 1.

$|\psi_{<}^{e}|^{2}=0.55;\,|\psi_{<}^{x}|^{2}=0;\,|\psi_{>}^{e}|^{2}=0;\,|\psi_{<}^{x}|^{2}=0.45;\,|\chi_{<}^{e}|^{2}=0.45;\,|\chi_{<}^{x}|^{2}=0;\,|\chi_{>}^{e}|^{2}=0;\,|\chi_{<}^{x}|^{2}=0.55.$
(43)

In this case the flavor populations are unequal in both directions. The
results are shown in Fig. 3 In this case, the GRK and UE integration method
converges to $O(h^{4})$ using step sizes of 0.5 cm or less, but the UE method
does not converge even at these small step sizes.

  * •

Case 4: In this case, we use a two-energy spectrum

$E_{<}=9~{}{\rm MeV};\,E_{>}=11~{}{\rm MeV}$ (44)

and an asymmetric initial population:

$|\psi_{<}^{e}|^{2}=0.55;\,|\psi_{<}^{x}|^{2}=0;\,|\psi_{>}^{e}|^{2}=0;\,|\psi_{<}^{x}|^{2}=0.45;\,|\chi_{<}^{e}|^{2}=0.45;\,|\chi_{<}^{x}|^{2}=0;\,|\chi_{>}^{e}|^{2}=0;\,|\chi_{<}^{x}|^{2}=0.55.$
(45)

This is the same as Case 2 except that the energies are unequal as well. The
results are shown in Fig. 4 In this case, the GRK method converges to
$O(h^{4})$ until $x\sim 50~{}{\rm km}$, for which the rapid variations in the
solution are confirmed to $O(h^{4})$. Beyond that range of $x$ it was
difficult to obtain a stable solution. In an attempt to improve convergence,
we used a sixth-order implicit Runge-Kutta Algorithm Dieci et al. (1994);
Butcher (1964). With this improved algorithm, we were able to main stability
with a precision of $O(h^{6})$ up to about $x\sim 75$ km, but could not move
much further, even with step sizes $h\sim~{}1$ mm. The UE method does not
converge at any value of $x$.

Figure 5: Growth in relative separation of solutions for small initial
displacements, where $\epsilon$ is defined in the text. Curves are normalized
to unity at $x=0$.

Further investigation of Case 4 suggests that the behavior is in fact chaotic.
To see this, we consider slightly modified parameter sets where

$|\psi_{<}^{e}|^{2}=0.55+\epsilon;\,|\psi_{<}^{x}|^{2}=0;\,|\psi_{>}^{e}|^{2}=0;\,|\psi_{<}^{x}|^{2}=0.45;\,|\chi_{<}^{e}|^{2}=0.45-\epsilon;\,|\chi_{<}^{x}|^{2}=0;\,|\chi_{>}^{e}|^{2}=0;\,|\chi_{<}^{x}|^{2}=0.55.$
(46)

The displacement of initial condition $\delta|\psi|^{2}$ leads to a
displacement of the outcome probability $\delta P_{1}$. For a chaotic system,
we expect to see Alligood et al. (1996)

$|\delta P_{1}(x)|\approx e^{\lambda x}|\delta P_{1}(x=0)|.$ (47)

In Fig. 5 we show a calculated normalized displacement:

$\delta p_{1}(x)=\left|\delta P_{1}(x)\over\delta P_{1}(0)\right|$ (48)

for three choices of $\epsilon$: $10^{-8}$, $10^{-7}$, $10^{-6}$. By
definition the curves begin at unity. To the extent that Eq. 47 is correct,
the curves will lie on top of one another. In practice, all three trajectories
rise to $p_{1}\sim 100$ and remain near that value until $x\sim 0$ km, after
which they rise at an exponential rate. For comparison, see Fig. 2.32 in Ref.
Alligood et al. (1996).

###  III.3 Discussion

From the examples described above, there are cases in the parameter space of
initial conditions that lead to stable solutions using either GRK or UE
methods, stable solutions using only the GRK method, and solutions that
ultimate become chaotic, as confirmed by very high precision GRK algorithms.
That said, the cases yielding, convergent, non-chaotic solution represent very
special instances in which some parameters are equal or satisfy a symmetry.
The more general case that one would encounter corresponds to Case 4, where
the neutrinos have unequal energies, and the $A$ and $B$ $\nu_{e}$ and
$\nu_{x}$ populations are also not the same. The transition to a chaotic
solution also occurs well before the neutrinos would be expected to reach
their asymptotic mass contents.

These results are similar to those found in Refs. Mirizzi and Serpico (2012);
Sarikas et al. (2012); Mirizzi (2013); Raffelt et al. (2013); Raffelt and
Seixas (2013), which examine sensitivities of solutions under very small
variations of initial conditions. In particular, the one-dimensional model of
Ref. Raffelt and Seixas (2013) produces very similar behavior to that found
here. That study found that results were stable for high neutrino densities.
Consistent with that study, we find that very high neutrino densities yield
stable results. However, we also find that as the neutrino density falls off,
the results become chaotic. Instabilities were also observed in a model study
in Ref. Sawyer (2005).

In this study we distinguish between numerical stability, which can eventually
be achieved to arbitrary precision using powerful integration algorithms, and
chaotic behavior, which is a measure of the physical sensitivity of the
results to initial condition. As noted in these references, it is not clear
whether a full three-dimensional approach using an algorithm that has proven
convergence as a function of step size can yield stable results when initial
conditions are asymmetric. Even if this is possible, there remains the second
issue of sensitivity, whose impact is physical rather than numerical.

##  IV Summary

With a one-dimensional model we have been able to explore several aspects of
non-linear neutrino evolution outside a supernova neutrinosphere. For certain
limiting cases (such as a single-flavor initial state), one can obtain exact
asymptotic results that differ from what one expects from the corresponding
linear case (such as solar neutrino evolution). For more general cases of
initial flavor distributions, neutrino flavor (or mass) evolution becomes
difficult to integrate numerically; one can achieve numerical convergence, but
the solution and eventually exhibits the standard signals for chaotic
behavior. We note here that the framework of refractive optics remains valid
even in the chaotic region. The large neutrino mean-free path means that this
system is in fact an excellent illustration of dynamical birefringence whose
microscopic origin can be derived, whether or not a stable numerical solution
exists.

###### Acknowledgements.

The author wishes to thank A. B. Balantekin, J. F. Cherry, G. M. Fuller, J. H.
Eberly, K. E. Keister, and P. Marronetti for helpful discussions and
correspondence. This work was supported by the National Science Foundation
through its employee Independent Research and Development Program. The
opinions and conclusions expressed herein are those of the author and do not
represent the National Science Foundation

## References

  * Duan et al. (2010) H. Duan, G. M. Fuller, and Y.-Z. Qian, Ann. Rev. Nucl. Part. Sci. 60, 569 (2010).
  * Pantaleone (1992) J. Pantaleone, Physics Letters B 287, 128 (1992).
  * Notzold and Raffelt (1988) D. Notzold and G. Raffelt, Nucl. Phys. B307, 924 (1988).
  * Samuel (1993) S. Samuel, Phys. Rev. D 48, 1462 (1993).
  * Samuel (1996) S. Samuel, Phys. Rev. D 53, 5382 (1996).
  * Pastor and Raffelt (2002) S. Pastor and G. Raffelt, Phys. Rev. Lett. 89, 191101 (2002).
  * Hannestad et al. (2006) S. Hannestad, G. G. Raffelt, G. Sigl, and Y. Y. Y. Wong, Phys. Rev. D 74, 105010 (2006).
  * Dasgupta et al. (2009) B. Dasgupta, A. Dighe, G. G. Raffelt, and A. Y. Smirnov, Phys. Rev. Lett. 103, 051105 (2009).
  * Pastor et al. (2002) S. Pastor, G. Raffelt, and D. V. Semikoz, Phys. Rev. D 65, 053011 (2002).
  * Balantekin and Pehlivan (2007) A. B. Balantekin and Y. Pehlivan, Journal of Physics G: Nuclear and Particle Physics 34, 47 (2007).
  * Dasgupta et al. (2012) B. Dasgupta, E. P. O’Connor, and C. D. Ott, Phys. Rev. D 85, 065008 (2012).
  * Friedland (2010) A. Friedland, Phys. Rev. Lett. 104, 191102 (2010).
  * Duan et al. (2011) H. Duan, A. Friedland, G. C. McLaughlin, and R. Surman, Journal of Physics G: Nuclear and Particle Physics 38, 035201 (2011).
  * Mirizzi and Serpico (2012) A. Mirizzi and P. D. Serpico, Phys. Rev. Lett. 108, 231102 (2012).
  * Raffelt and Seixas (2013) G. Raffelt and D. d. S. Seixas, Phys. Rev. D 88, 045031 (2013).
  * Raffelt et al. (2013) G. Raffelt, S. Sarikas, and D. d. S. Seixas, Phys. Rev. Lett. 111, 091101 (2013).
  * Duan and Shalgar (2014) H. Duan and S. Shalgar (2014), eprint arXiv:1412.7097.
  * Wolfenstein (1978) L. Wolfenstein, Phys. Rev. D17, 2369 (1978).
  * Mikheev and Smirnov (1985) S. P. Mikheev and A. Y. Smirnov, Sov. J. Nucl. Phys. 42, 913 (1985).
  * van de Hulst (1949) H. C. van de Hulst, _Light Scattering By Small Particles_ (Wiley, New York, 1949).
  * Duan et al. (2006) H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Phys.Rev. D74, 105014 (2006).
  * Dieci et al. (1994) L. Dieci, R. Russell, and E. Van Vleck, SIAM J. Num. Anal. 31, 261 (1994).
  * Butcher (1964) J. C. Butcher, Math. Comp. 18, 50 (1964).
  * Alligood et al. (1996) K. T. Alligood, T. Sauer, and J. A. Yorke, _Chaos: An Introduction To Dynamical Systems_ (Springer, New York, 1996).
  * Sarikas et al. (2012) S. Sarikas, D. de Sousa Seixas, and G. Raffelt, Phys. Rev. D 86, 125020 (2012).
  * Mirizzi (2013) A. Mirizzi, Phys. Rev. D 88, 073004 (2013).
  * Sawyer (2005) R. F. Sawyer, Phys. Rev. D 72, 045003 (2005).