# Inference solves a boundary-value collision problem, with relevance to
neutrino flavor transformation

Eve Armstrong<EMAIL_ADDRESS>Department of Physics, New York
Institute of Technology, New York, NY 10023, USA  Department of Astrophysics,
American Museum of Natural History, New York, NY 10024, USA

###### Abstract

Understanding neutrino flavor transformation in dense environments such as
core-collapse supernovae (CCSN) is critical for inferring nucleosynthesis and
interpreting a detected neutrino signal. The role of direction-changing
collisions in shaping the neutrino flavor field in these environments is
important and poorly understood; it has not been treated self-consistently.
There has been progress, via numerical integration, to include the effects of
collisions in the dynamics of the neutrino flavor field. While this has led to
important insights, integration is limited by its requirement that full
initial conditions must be assumed known. On the contrary, feedback from
collisions to the neutrino field is a boundary value problem. Numerical
integration techniques are poorly equipped to handle that formulation. This
paper demonstrates that an inference formulation of the problem can solve a
simple collision-only model representing a CCSN core – without full knowledge
of initial conditions. Rather, the procedure solves a two-point boundary value
problem with partial information at the bounds. The model is sufficiently
simple that physical reasoning may be used as a confidence check on the
inference-based solution, and the procedure recovers the expected model
dynamics. This result demonstrates that inference can solve a problem that is
artificially hidden from integration techniques – a problem that is an
important feature of flavor evolution in dense environments. Thus, it is
worthwhile to explore means of augmenting the existing powerful integration
tools with inference-based approaches.

## I. Introduction

The physics of flavor evolution in the neutrino field in high-density
environments such as core collapse supernovae (CCSN) can significantly affect
the transport of energy, entropy, and lepton number, with implications
regarding nucleosynthesis, the mechanism of explosion, and the mass assembly
histories of galaxies Mirizzi  _et al._ (2016); Duan  _et al._ (2010);
Chakraborty  _et al._ (2016). Due to the fierce nonlinearity of various
features of this problem, we are far from a complete understanding. One
important and poorly studied problem is the role of direction-changing
collisions in shaping the neutrino density profile in these environments.

It is well established that coherent forward scattering of neutrinos with
matter can lead to resonant effects Mikheev and Smirnov (2007, 1985);
Wolfenstein (1978), and that forward scattering of neutrinos with each other
also affects flavor evolution Fuller  _et al._ (1987); Nötzold and Raffelt
(1988); Pantaleone (1992a, b). These potentials, which govern neutrino flavor
transformation, themselves depend on the flavor states of the neutrinos. This
nonlinearity has led to the development of highly sophisticated numerical
approaches, including the multi-angle Neutrino BULB code Duan  _et al._ (2006,
2008) and the IsotropicSQA code Richers  _et al._ (2019).

The ability of direction-changing scattering, in addition, to significantly
impact flavor evolution in the SN envelope was first shown by Ref Cherry  _et
al._ (2012). There have since been efforts to understand this "halo effect" on
flavor transformation Cherry  _et al._ (2013); Cirigliano  _et al._ (2018);
Zaizen  _et al._ (2020); Cherry  _et al._ (2020), although a self-consistent
solution has not been obtained. More recent studies have shown that fast
pairwise conversions can happen due to a crossing between the angular
distributions of neutrinos and anti-neutrinos (See Ref Tamborra and Shalgar
(2020) and citations therein) – a phenomenon that is sensitive to the exact
shapes of the angular distributions.

Most recently, direction-changing collisions (hereafter "collisions") have
been shown to enhance fast flavor conversion in high-density environments
Shalgar and Tamborra (2021); Johns (2021); Sasaki and Takiwaki (2021). This
result is counter to the expectation that collisions should dampen flavor
evolution via dynamical decoherence Martin  _et al._ (2021). Moreover, even in
small-scale CCSN models, accounting for collisions has yielded significant and
counter-intuitive changes in the angular distributions of the neutrino density
profiles. For a more comprehensive summary of the role of collisions in flavor
evolution, see Ref Shalgar and Tamborra (2021) and citations therein.

To date, efforts to include collisions in the neutrino quantum kinetic
equations have used numerical integration to solve an initial value problem
(e.g. Refs Cirigliano  _et al._ (2018); Shalgar and Tamborra (2019, 2021)). It
is not necessarily desirable to place assumptions on initial conditions. On
the contrary, when feedback between collisions and flavor evolution is
included in the calculation, the problem becomes a two-point boundary-value
problem with no guarantee of a unique solution Shalgar and Tamborra (2021).

With that motivation, in this paper we solve a simple collision model, not as
an initial-value problem, but rather as a two-point boundary value problem
with partial information at the boundaries. We do this via an inference – or
"inverse" – formulation Tarantola (2005) of the problem. Inference is a means
to optimize a model with data, to predict model evolution at locations outside
of those where the data have been provided. The inference formulation does not
require that initial conditions be known. Rather, it requires that some
condition(s) (or constraint(s)) be known at at one or more locations that
parameterize the model dynamics. These conditions may be placed at the
endpoints of model evolution, or more generally at any point(s) along the
model trajectory. Further, the formulation of inference employed in this paper
offers a mean to determine whether a unique solution exists, and – if it does
not – the degree of degeneracy present.

In this paper, we use a formulation of inference – statistical data
assimilation (SDA) – that is built to handle the case of extremely sparse
data. SDA was invented for numerical weather prediction Kimura (2002); Kalnay
(2003); Evensen (2009); Betts (2010); Whartenby  _et al._ (2013); An  _et al._
(2017) and has gained considerable traction in neurobiology Schiff (2009);
Toth  _et al._ (2011); Kostuk  _et al._ (2012); Hamilton  _et al._ (2013);
Meliza  _et al._ (2014); Nogaret  _et al._ (2016); Armstrong (2020). Within
astrophysics, inference has been used mainly for pattern recognition (e.g.
Djorgovski  _et al._ (2006)), while its utility for model completion is
gaining notice in the exoplanet community Madhusudhan (2018) and solar physics
Kitiashvili and Kosovichev (2008). Recently, SDA has been applied to small-
scale models of neutrino flavor evolution Armstrong  _et al._ (2017, 2020);
Rrapaj  _et al._ (2021). In those papers, the model contained coherent forward
scattering only, and thus numerical integration could be taken as a confidence
check on SDA solutions.

By contrast to Refs Armstrong  _et al._ (2017, 2020); Rrapaj  _et al._ (2021),
this work tackles a boundary-value problem formulation, which numerical
integration cannot access. For this reason, we shall examine a simple model
whose solution can be intuited via physical reasoning. The model represents
the high-density region of a supernova core, wherein neutrinos are effectively
trapped and distributed isotropically. In this context, we consider collisions
alone: absent neutrino-matter and neutrino-neutrino coupling potentials. The
model consists of four neutrino beams: two outgoing and two incoming, and four
angular bins. Absent flavor evolution, we will seek to examine the angular
distributions of the neutrino number densities as functions of radius, due to
collisions alone.

Given partial information at the boundaries, the SDA procedure is able to
recover the density profiles that follow from a simple physical argument,
thereby illuminating physics that is artificially hidden from numerical
integration techniques. Further, multiple independent trials converge to one
solution.

This paper proceeds as follows. Section II describes the model and imposed
boundary conditions, and presents a simple physical argument regarding the
expected model evolution. Section III describes the inference methodology and
the simulated experiments performed. Section IV shows how the SDA procedure
predicts the expected result. Finally, Section V comments on the implications
for tackling the full quantum kinetic equations that includes nonlinear
feedback between collisions and flavor evolution, and advocates for exploring
means to augment existing integration methodologies with inference-based
techniques.

## II. Stationary box model of a supernova core

The goal in crafting the model is twofold: 1) to create an opportunity to
observe the effect on the neutrino angular distributions due to collisions
alone, and 2) to retain a simplicity such that the solution can be intuited
via physical reasoning, so as to provide a confidence check on the SDA
solution. To this end, we employ a two-dimensional flat "box" model and four
neutrino beams: two incoming and two outgoing.

The model schematic of the supernova core is shown in Figure 1. Neutrinos
radiate from the center (the x-axis at $y=0$) out to some final radius $R$ (at
$y=y_{max}$). To examine the effects of collisions alone, we omit flavor
evolution and adopt a single-flavor model, with no anti-neutrinos. There is a
sink at $r=R$ where the outgoing beams may escape. There is also a source and
sink permitted at $y=0$, as particles are permitted to travel in the negative
$\hat{y}$ direction.

Neutrino flavor evolution can be expressed in terms of a density matrix $\rho$
(and $\overline{\rho}$ for anti-neutrinos). For each neutrino momentum mode
$\vec{p}$, we can write the following quantum-kinetic equation Volpe (2015);
Vlasenko  _et al._ (2014):

$\displaystyle\begin{split}\left(\frac{\partial}{\partial
t}+\vec{v}\cdot\nabla\right)\rho(\vec{x},\vec{p},t)=-i[H(\vec{x},\vec{p},t),\rho(\vec{x},\vec{p},t)]\\\
+C\\{\rho(\vec{x},\vec{p},t),\overline{\rho}(\vec{x},\vec{p},t)\\},\end{split}$
(1)

Figure 1: Two-dimensional "box" model of the supernova core. The core center
is the x-axis, at $r=0$. The positive $\hat{y}$ direction is outward toward
the edge. Two beams are directed outward, and two inward, each along some
angle $\theta$.

where $\vec{x}$ and $t$ are spatial and temporal coordinates, respectively. On
the right side, the Hamiltonian term $H$ describes neutrino mixing in vacuum,
neutrino interactions with the matter background, and neutrino-neutrino self
interactions; it contains coherent scattering. The second term $C$ is the
"collision term," which accounts for direction-changing scattering111Note that
the separation of coherent from direction-changing scattering is an artificial
one; see Section V..

In this paper, we neglect the (coherent) Hamiltonian term and focus on the
collision term $C$. We shall assume that a steady-state can be defined on a
timescale shorter than the leakage timescale, and so we eliminate time from
Equation 1. That is, we will solve Equation 1 depending on scattering angle
and not on momentum. In this context, the velocity in the second operator on
the left side shall be considered a spatial – rather than a temporal –
velocity. Then our equation of motion simplifies to:

$\textrm{cos}\theta_{i}\diff{\rho_{i}(r)}{r}=\sum_{j}\left(-C^{\textrm{loss}}\rho_{i}(r)+C^{\textrm{gain}}\rho_{j}(r)\right)\frac{\Delta\textrm{cos}\theta_{j}}{2}.$
(2)

Here, subscripts ($i,j$) denote indices of angular bins,
$\Delta\textrm{cos}\theta_{j}$ is the width of the jth angular bin, and the
cosine term on the left side accounts for the dependence of path length on
zenith angle. $C^{\textrm{loss}}$ and $C^{\textrm{gain}}$ are the coefficients
for loss and gain for each bin; they represent the strength of interaction
between neutrinos and nucleons222Physically, the loss and gain coeffients
represent the product of the number of scatterers and the cross-section for
each particular interaction channel averaged over some distribution of
neutrino energy Shalgar and Tamborra (2019).. Thus for our box model, there
are four Equations 2, each governing the evolution of the neutrino density in
one of four angular bins.

To create a problem with a physically intuitable solution, we make further
simplifications. First, the coefficients $C^{\textrm{loss}}$ and
$C^{\textrm{gain}}$ are constant numbers throughout the medium. That is, the
matter density is constant within the box and drops discontinuously to zero at
$r=R$. Next we set $C^{\textrm{loss}}$ equal to $C^{\textrm{gain}}$, to
conserve particle number at each interaction333Note that setting
$C^{\textrm{loss}}$ equal to $C^{\textrm{gain}}$ renders Equation 2 trivial,
unless there are slight asymmetries in the values of the $\textrm{cos}\theta$
terms on the left side of the equations. For this reason, we introduce slight
asymmetries between the fluxes of outgoing versus incoming beams. The
$\textrm{cos}\theta$ values are 0.8 and 0.9 for the outgoing beams, and -0.6
and -0.7 for those incoming..

Finally, we choose a value of $C$ to approximate nuclear density. Near nuclear
density, neutrinos are effectively trapped, and so we expect an isotropic
distribution of scatterers: at each collision, a neutrino has equal
probability of scattering inward as outward. To choose $C$ appropriately, we
perform a dimensional analysis of Eq 2. The left side is a derivative, with
units of 1/length. The right side also has units of 1/length, as the only
dimensionful quantities on the right side are the coefficients $C$, with units
of inverse length. Physically, as $C$ represents the strength of the
interaction between neutrinos and matter, $C$ can be considered the inverse of
the neutrino mean free path, with units of $m^{-1}$. For example, taking a
mean free path on the order of meters to represent nuclear density, we choose
$C$ to be 0.2 $m^{-1}$.

### A. The boundary conditions

We place on this model partially-known boundary conditions. First, the values
of the outgoing beams at the core center ($r=0$) are known. At that location,
we give them both a (dimensionless) number density of 0.5. Second, the values
of incoming beams at the core edge ($r=R$) are zero. That is, by definition,
at $r=R$, there occurs no further back-scattering. This is the sole
information furnished to the inference procedure. The procedure is tasked with
taking this information together with the model equations of motion (Eq 2), to
predict the complete density profile at all other spatial locations for all
four angular bins – including the values at $r=0$ of the incoming beams and
the values at $r=R$ of the outgoing beams.

### B. Expected density profile, to be recovered by the SDA procedure

Given a value of $C$ within our box that corresponds to nuclear density, and
the fact that a system represented as ordinary differential equations has no
memory, we expect to observe isotropy for the vast majority of the
trajectories of all four beams. That is, at each interaction the particles
possess no preferred direction. Thus, given the condition that the outgoing
beams begin their trajectory at $r=0$ with a value of 0.5, we expect the
densities of all beams to remain approximately 0.5 throughout the majority of
their their trajectories. Then, as we approach $r=R$, we should see a
precipitous drop. Near $R$ the densities of the incoming beams should drop to
zero, in accordance with the imposed boundary condition there, and the
inference procedure should predict that the densities of the outgoing beams
drop to some number below 0.5 – that number to scale inversely with the value
of $C$. Stated physically: as $C$ weakens, more particles should escape.

Note again that the model possesses a leak at $r=R$ (and at $r=0$, as
scatterers may pass beyond zero – in the negative $\hat{y}$ of Figure 1. We do
not expect that a steady-state solution exists for longer than a duration
corresponding to the leakage timescale. Rather, we expect that different
choices for the coefficient $C$ will yield different solutions, each
representing one freeze-frame of a steady state solution. Specifically, as $C$
decreases, the mean free path increases – and equivalently the total duration
increases, as increasingly more neutrinos have had time to escape the core.

## III. Inference methodology

Statistical data assimilation (SDA) is an inference procedure in which a
dynamical system is assumed to underlie any measured quantities. This model
$\bm{F}$ can be written as a set of D ordinary differential equations that
evolve in some parameterization $r$ as:

$\displaystyle\diff{x_{a}(r)}{r}$
$\displaystyle=F_{a}(\bm{x}(r),\bm{p}(r));\hskip 10.00002pta=1,2,\ldots,D,$

where the components $x_{a}$ of the vector x are the model state variables.
Any unknown parameters to be estimated are contained in $\bm{p}$, and may
themselves be variable. In this paper, model $\bm{F}$ is the set of four
equations (Eq 2) governing the neutrino density profile in each of four
angular bins.

A subset $L$ of the $D$ state variables is associated with measured
quantities. In this paper, the "measurements" are the four boundary conditions
placed on the outgoing beams at $r=0$ and incoming beams at $r=R$. One seeks
to estimate the evolution of all state variables that is consistent with the
measurements provided, to predict model evolution at parametrized locations
where measurements are not present.

We formulate the SDA procedure as an optimization wherein a cost function is
extremized. An optimization formulation of inference does not require
knowledge of initial conditions. Rather, it requires that some constraints be
placed on the model at some location(s) on the coordinate axis that
parameterizes the model equations of motion. These conditions may exist at the
endpoints of the problem, as is the case in this paper. Importantly, however,
the conditions need not be at the bounds. This flexibility is the key
advantage of inference, compared to integration, that we aim to demonstrate.

The cost function is written in two terms444Additional terms representing
equality constraints may be added to the cost function, depending on the aim
of a particular optimization procedure.. One term represents "measurement
error": the difference between state prediction and any measurements made. The
second term represents "model error": the difference between state prediction
and adherence to the model dynamics. It will be shown in Section IV that
treating the model error as finite offers a means to determine whether a
particular solution is consistent with both measurements and model dynamics,
as well as a means to assess uniqueness.

We search the surface of the cost function via the variational method. The
procedure in its entirety - that is: a variational approach to minimization
coupled with an annealing method to identify a lowest minimum of the cost
function (which will be described below in this Section) - is referred to as
variational annealing (VA). The procedure searches a
$(D\,(N+1)+p)$-dimensional state space, where N is the number of discretized
steps. One seeks the path
$\bm{X}^{0}={\bm{x}(0),...,\bm{x}(N),\bm{p}(0),...\bm{p}(N)}$ in state space
on which the cost function attains a minimum value. We refer to the cost
function $A_{0}$ as the action, because it can be derived from the concept of
a classical action of a particle on a path in a state space Abarbanel (2013).
Ref Armstrong  _et al._ (2020) demonstrated that the action formulation offers
a simple litmus test for identifying correct solutions: namely, they are
solutions that correspond to the path of least action.

The cost function $A_{0}$ used in this paper is written as:

$\begin{split}A_{0}=&R_{f}A_{model}+R_{m}A_{meas}\\\
A_{model}=&\frac{1}{{N}D}{\sum}_{n\in\\{\text{odd}\\}}^{N-2}{\sum}_{a=1}^{D}\left[\left\\{x_{a}(n+2)-x_{a}(n)-\frac{\delta
r}{6}[F_{a}(\bm{x}(n),\bm{p})+4F_{a}(\bm{x}(n+1),\bm{p})+F_{a}(\bm{x}(n+2),\bm{p})]\right\\}^{2}\right.\\\
&\hskip
100.0pt+\left.\left\\{x_{a}(n+1)-\frac{1}{2}\left(x_{a}(n)+x_{a}(n+2)\right)-\frac{\delta
r}{8}[F_{a}(\bm{x}(n),\bm{p})-F_{a}(\bm{x}(n+2),\bm{p})]\right\\}^{2}\right]\\\
A_{meas}=&\frac{1}{N_{\text{meas}}}{\sum}_{j}{\sum}_{l=1}^{L}(y_{l}(j)-x_{l}(j))^{2}.\end{split}$
(3)

The second (simpler) squared term of Equation 3 governs the transfer of
information from measurements $y_{l}$ to model states $x_{l}$555This term
derives from the mutual information of probability theory Abarbanel (2013)..
Here, the summation on j runs over all discretized locations $J$ at which
measurements are made, which may be some subset of all discretized steps of
the model. The summation on l is taken over all L measured quantities. In the
simulations of this paper, the measured quantities are the number densities in
the angular bins, the outgoing beams sampled at radius of $r=0$ (at $j=1$
only), and the incoming beams sampled only at $r=R$ (at the value of $j$ that
denotes the final radial location $R$).

The first squared term of Equation 3 incorporates the model evolution of all D
state variables $x_{a}$. Here, the outer sum on n is taken over all
discretized timepoints of the model equations of motion. The sum on a is taken
over all D state variables at all discretized locations.666This term can be
derived via consideration of Markov-chain transition probabilities Abarbanel
(2013). In our model, these measured quantities are all four beams, or
$L=D=4$. The first and second bracketed terms represent error in the first and
second derivative of the model, respectively. For further details on the
action formulation, see Appendix A of Ref Armstrong  _et al._ (2017).

Importantly, note that the SDA procedure is tasked with inferring from the
sparse boundary conditions the complete model evolution at all locations on
$r$.

Finally, if the model is nonlinear, the cost function surface will be non-
convex: multiple minima may exist. For this reason, we employ an annealing
procedure to identify the global minimum of the problem, and to ascertain
whether there exists a single global minimum – that is, a unique solution.
Details of this annealing procedure will be presented in Section IV.

### A. Details of the simulated experiments

We performed the simulated experiment for three distinct values of the
coefficient $C^{\textrm{loss}}=C^{\textrm{gain}}$: 2.0, 0.2, and 0.02
$m^{-1}$, where the choice of 0.2 $m^{-1}$ represents nuclear density (a mean
free path of 5 meters). At each discretized model location $r$, the procedure
was permitted to search the full dynamical range of each state variable of
[0:1]. For each of those three experiments, we took the size of the core to be
$R=50km$, and examined the robustness of the solution across three distinct
choices for step size $dr$: 1 m, 10 m, and 100 m777Within the dimensionless
mathematical framework of the SDA procedure, those choices of 1, 10, and 100
$m^{-1}$ corresponded to: a step size of 0.01 with 10101 steps, 0.1 with 1011
steps, and 1.0 with 101 steps, respectively – keeping to a total size of 101..

To examine the uniqueness of solutions, for each of the (nine) experiments
described, 20 independent paths were searched, each beginning at a randomly-
generated set of initial guesses for the four state variables at each
parameterized location $r$.

To perform the optimization, we used the open-source Interior Point Optimizer
(Ipopt) Wächter (2009). Ipopt uses a Simpson’s rule method of finite
differences to discretize the state space, a Newton’s method to search, and a
barrier method to impose user-defined bounds that are placed upon the
searches. The discretization of the state space, the calculations of the model
Jacobean and Hessian matrices, and the annealing procedure were performed via
an interface with Ipopt that was written in C and Python minAone interface
with Interior-point Optimizer . All simulations were run on a 720-core,
1440-GB, 64-bit CPU cluster.

## IV. Result

### A. General findings

Results are fourfold.

  * •

As expected for an isotropic distribution of scatterers, for a large value of
the coefficient $C^{\textrm{loss}}=C^{\textrm{gain}}$, the number densities of
the four angular bins remain constant at 0.5 for the majority of the
trajectory, with a precipitous drop near the bound at $r=R$. As the value of
$C$ decreases, the rate of depletion of the number densities as a function of
$r$ slows, and the predicted values for the outgoing beams at $r=R$ show that
the densities in those bins have been depleted to a lesser degree, as more
particles escape the core. In the limit of no collisions ($C=0$), the solution
to Equation 2 is trivial and all particles escape.

  * •

The results are invariant, to one part in $10^{6}$, across choices of
discretized step size that span three orders of magnitude.

  * •

When the procedure receives one additional constraint such that the rates of
depletion of the bin populations near $R$ are captured, and meanwhile the
value of $C$ is withheld, the procedure is able to correctly infer the value
of $C$ that corresponds to those rates of depletion.

  * •

For all simulations performed, the plot of the cost function over the course
of an annealing procedure (to be described below) shows that i) all results
converge; that is: solutions are consistent with both boundary conditions and
model dynamics, and ii) each solution is unique; that is: all paths searched
converge to one solution.

### B. Effect of collision strength $C$

We examined the solution for three distinct values of collision strength
$C^{\textrm{loss}}=C^{\textrm{gain}}$, to determine whether the optimization
would recover our expectation that the rate of depletion of the number
densities in the angular bins should slow as $C$ weakens. Figure 2 shows the
result. The left, middle, and right columns correspond to values of $C$ of
2.0, C=0.2, and 0.02 $m^{-1}$, respectively. The top and bottom rows show the
density profile as a function of $r$ for one outgoing beam $\rho_{1}$ and one
incoming beam $\rho_{2}$, respectively. (The trajectories of the other two
beams in the four-beam model are identical; not shown.) The circles denote the
locations of the "measurements", or boundary conditions, that were given as
information to the SDA procedure. The remainder of the trajectories are the
SDA prediction, given knowledge of the model dynamics.

From left to right in Figure 2, note three effects. First, the number
densities near $r=0$ appear roughly the same, which follows from the argument
regarding isotropy presented in Section II.

Second, as $C$ weakens, the number densities fall off near $r=R$ increasingly
gradually. Specifically, the predicted value for the outgoing beam $\rho_{1}$
at $r=R$ is 0.114 (left), 0.115 (middle), and 0.158 (right). This follows from
the expectation that as the mean free path lengthens, more particles escape.

Third, as $C$ weakens, fewer incoming particles are predicted to reach $r=0$.
Specifically, the value for the incoming beam $\rho_{2}$ at $r=0$ is 5.0
(left), 4.9 (middle), and 4.4 (right). These numbers are summarized in Table
1.

Recalling that the model contains a leak at $r=R$, we can consider the left,
middle, and right panels of Figure 2 to be a succession of solutions in time,
each approximating a steady state over a timescale shorter than the leakage
timescale. As $C$ weakens from left to right, the time increases. In the limit
where $C=0$, the solution to Equation 2 is trivial, and all particles escape
as time goes to infinity (not shown).

Finally, we make several notes regarding the robustness of the results. The
procedure imposes the boundary conditions ($\rho_{1}=0.5$ at $r=0$ and
$\rho_{2}=0.0$ at $r=R$) to one part in $10^{7}$. The results are invariant,
to one part in $10^{6}$, across the three choices of discretized step size
$dr$ corresponding to 1, 10, and 100 m888For the reader interested in mapping
these choices to a realistic physical scenario: a choice of $C=0.2m^{-1}$,
which corresponds to a mean free path of 5m, together with a step size $dr$ of
1 m, would adequately sample a particle’s path.. The limit on that accuracy is
likely due to a discretization error internal to Ipopt, the optimizer used in
the procedure.

Figure 2: SDA result for the density profiles as functions of radius $r$ of one outgoing (top) and one incoming beam (bottom), for three values of coefficients $C$. For left, middle, and right columns: $C=$ 2.0, 0.2, and 0.02 $m^{-1}$, respectively. Circles denote locations of "measurements", or boundary conditions, provided to the procedure. The remaining trajectories are the SDA prediction. See the text for an interpretation in light of the expectation set forth in Section II. C ($m^{-1}$) | 2.0 | 0.2 | 0.02
---|---|---|---
$\rho_{1}$ at $r=R$ (outgoing) | 0.114 | 0.115 | 0.158
$\rho_{2}$ at $r=0$ (incoming) | 0.5 | 0.49999 | 0.44
Table 1: Predicted values of outgoing and incoming beams at the endpoints,
across values of $C$. As $C$ dilutes from left to right, top: more outgoing
neutrinos escape to $r=R$; bottom: fewer incoming neutrinos return to $r=0$.

### C. Correct inference of collision strength $C$ as an unknown parameter

Having determined that the SDA procedure can predict the rate of depletion of
the angular bins given a particular strength of coefficients $C$, we sought to
determine whether it could perform the inverse: to infer $C$ given a rate of
depletion.

To this end, we repeated the experiment leaving $C$ as an unknown parameter to
be estimated. As before, we gave the SDA procedure the original four boundary
conditions. In addition, from the original solutions shown in Figure 2 we gave
to the procedure the values of densities in each angular bin at five
discretized locations prior to $R$, with the aim of capturing the rate of
depletion of the bins. For each of the three cases, the procedure correctly
estimated the unknown value of $C$ to sixth-decimal precision.

### D. Convergence and uniqueness

For a nonlinear model, the surface of the action, or cost function, will be
non-convex. The complete SDA procedure anneals in terms of the ratio of model
and measurement error, $R_{f}$ and $R_{m}$, respectively999More generally,
$R_{m}$ and $R_{f}$ are inverse covariance matrices for the measurement and
model errors, respectively. In this paper the measurements are taken to be
mutually independent, rendering these matrices diagonal., to gradually freeze
out a lowest-minimum of the Action Ye  _et al._ (2015). This iteration works
as follows.

We define the coefficient of measurement error $R_{m}$ to be 1.0, and write
the coefficient of model error $R_{f}$ as: $R_{f}=R_{f,0}\alpha^{\beta}$,
where $R_{f,0}=10^{-1}$, $\alpha=1.5$, and $\beta$ is initialized at zero.
Parameter $\beta$ is the annealing parameter. When $\beta=0$, relatively free
from model constraints the Action surface is smooth and there exists one
minimum of the variational problem that is consistent with the measurements.
We estimate that minimum. Then we increase the weight of the model term
slightly, via an integer increment in $\beta$, and recalculate the action. We
do this recursively toward the deterministic limit of

Figure 3: Action, or cost function, as function of annealing parameter
$\beta$, for a value of $C$ of 0.2; representative of all simulations. As
$\beta$ passes 7, the action flatlines, indicating convergence to a solution
that is consistent with both measurements (or boundary conditions) and model
dynamics.

$R_{f}\gg R_{m}$. The aim is to remain sufficiently near to the lowest minimum
so as not to become trapped in a local minimum as the surface acquires the
structure imposed by the model dynamics.

Ref Armstrong  _et al._ (2020) demonstrated that, for a particular path
searched, the plot of action as a function of annealing parameter $\beta$
indicates whether a particular solution is stable as the model dynamics are
imposed increasingly rigidly. Figure 3 shows this result for the case of
$C=0.2$, representative of all experiments performed in this paper. Note that
at a value of $\beta=7$, the action flatlines. This behavior indicates that a
solution has been found that is consistent with both imposed boundary
conditions and with model dynamics. For a detailed explanation, see Ref
Armstrong  _et al._ (2020).

Finally, the 20 paths initialized for each experiment in this paper each
converged to a single solution. This was demonstrated by identical state
variable evolution (e.g. Figure 2) and identical action($\beta$) plots (e.g.
Figure 3) for each path sampled. Of course, to ascertain whether a solution is
truly unique, one must sample an infinite number of paths. The specific number
chosen for any particular model should depend on that model’s complexity and
dimensionality.

## V. Discussion

We have considered a simple two-point boundary value collision model with
partially-known boundary conditions, and challenged an inference procedure to
solve it. The model is sufficiently simple that a physical argument can intuit
the expected result, thereby serving as a confidence check on the inference-
based solution. Specifically, an optimization formulation of SDA recovers the
expected angular distributions of the density profiles of particles in a
collision-only flat two-dimensional model, where the particles represent
neutrino-nucleon interactions near nuclear density. The behavior of the cost
function over an iterative reweighting of the terms that impose the boundary
conditions and the equations of motion, respectively, demonstrates that the
SDA solution is consistent with both the boundary conditions and the model
dynamics. This was the case over a range of choices for the collision strength
that spanned three orders of magnitude, and the results are robust to the step
size used for discretization. Multiple randomly-initiated searches of the
state-and-parameter space converged to one solution.

This finding has important implications regarding the potential of inference
to inform flavor transformation models governed by the full quantum kinetic
equation (Equation 1; that is, taking the Hamiltonian together with direction-
changing collisions. While the integeration-based approaches to the collision
problem have led to valuable insights, considering the full feedback between
collisions and flavor evolution, will require that the problem be recast as a
boundary problem with no guarantee of a unique solution Shalgar and Tamborra
(2021). An inference formulation avoids reliance on initial conditions, and
rather seeks a solution that is compatible with constraints placed only at
locations where we have high confidence in our understanding of the physics
there. In the simulations performed in this paper, those constraints were
placed at the endpoints, but they need not be. Further, within the
optimization formulation, uniqueness can be investigated via the
initialization of multiple independent searches. The ability of this SDA
procedure to probe model degeneracy has been demonstrated in significantly
more detail in Refs Armstrong  _et al._ (2017, 2020).

It is relevant to note here another potential advantage of inference for
augmenting numerical integration tools: efficiency. Including feedback from
collisions self-consistently in the quantum kinetic calculations dramatically
increases the computational complexity, as both the collisions and the flavor
evolution will simultaneously shape the neutrino angular distributions Shalgar
and Tamborra (2021). The simple experiment presented in this paper did not
offer the opportunity to showcase the high efficiency with which this SDA
formulation performs state-and-parameter estimation; it has been demonstrated
for a small-scale flavor evolution model in Ref Armstrong  _et al._ (2020).

The question of how to implement inference within the existing numerical
integration framework has yet to be explored. It is relevant here to note that
the separation of coherent versus direction-changing scattering in Equation 1
– the former into the Hamiltonian and the latter into the collision term $C$ –
is an artificial choice. It might be worthwhile to consider possible
advantages of rewriting the quantum kinetic formulation to treat coherent and
direction-changing scattering collectively. Of course, this is a daunting
suggestion, as the the formulation of Equation 1 forms the foundation of
extremely large and powerful codes. Nevertheless, it is worth asking to what
extent that formulation limits our ability to access the full dynamical range
of the flavor field in compact object environments.

In parallel with exploring means to fold inference into the existing codes, it
will be instructive to continue examining how inference performs on collision-
only models of increasing complexity. For example, in the model used in this
paper, neutrino number density was conserved at each collision: collisions
were completely elastic and did not absorb or emit neutrinos. Further, the
matter density was held constant, with a discontinuous fall to zero at the
core’s edge. A more realistic matter profile will contain some radial
dependence. In addition, model behavior should be examined within a
spherically-symmetric geometry. Because adding complexity will obscure a
simple expectation regarding the result, it will be instructive to identify
alternative comparisons to inference-based solutions, such as random-walk
formulations.

## VI. ACKNOWLEDGEMENTS

Thank you to Shashank Shalgar for invaluable conversations, without which this
manuscript would not exist. Thank you also to George Fuller and Amol
Patwardhan for important comments. E. A. acknowledges the National Science
Foundation (NSF grant 2139004), an Institutional Support for Research and
Creativity grant from New York Institute of Technology, and Gretel Yeager.

## References

  * Mirizzi  _et al._ (2016) A. Mirizzi, I. Tamborra, H.-T. Janka, N. Saviano, K. Scholberg, R. Bollig, L. Hüdepohl, and S. Chakraborty, La Rivista del Nuovo Cimento 39, 1 (2016).
  * Duan  _et al._ (2010) H. Duan, G. M. Fuller, and Y.-Z. Qian, Annual Review of Nuclear and Particle Science 60, 569 (2010).
  * Chakraborty  _et al._ (2016) S. Chakraborty, R. Hansen, I. Izaguirre, and G. Raffelt, Nuclear Physics B 908, 366 (2016).
  * Mikheev and Smirnov (2007) S. Mikheev and A. Y. Smirnov, arXiv preprint arXiv:0706.0454 (2007).
  * Mikheev and Smirnov (1985) S. Mikheev and A. Y. Smirnov, Soviet Journal of Nuclear Physics 42, 913 (1985).
  * Wolfenstein (1978) L. Wolfenstein, Physical Review D 17, 2369 (1978).
  * Fuller  _et al._ (1987) G. M. Fuller, R. W. Mayle, J. R. Wilson, and D. N. Schramm, Astrophysical Journal 322, 795 (1987).
  * Nötzold and Raffelt (1988) D. Nötzold and G. Raffelt, Nuclear Physics B 307, 924 (1988).
  * Pantaleone (1992a) J. Pantaleone, Physics Letters B 287, 128 (1992a).
  * Pantaleone (1992b) J. Pantaleone, Physical Review D 46, 510 (1992b).
  * Duan  _et al._ (2006) H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Physical Review D 74, 105014 (2006).
  * Duan  _et al._ (2008) H. Duan, G. M. Fuller, and J. Carlson, Computational Science & Discovery 1, 015007 (2008).
  * Richers  _et al._ (2019) S. A. Richers, G. C. McLaughlin, J. P. Kneller, and A. Vlasenko, Physical Review D 99, 123014 (2019).
  * Cherry  _et al._ (2012) J. F. Cherry, J. Carlson, A. Friedland, G. M. Fuller, and A. Vlasenko, Physical review letters 108, 261104 (2012).
  * Cherry  _et al._ (2013) J. F. Cherry, J. Carlson, A. Friedland, G. M. Fuller, and A. Vlasenko, Physical Review D 87, 085037 (2013).
  * Cirigliano  _et al._ (2018) V. Cirigliano, M. Paris, and S. Shalgar, Journal of Cosmology and Astroparticle Physics 2018, 019 (2018).
  * Zaizen  _et al._ (2020) M. Zaizen, J. F. Cherry, T. Takiwaki, S. Horiuchi, K. Kotake, H. Umeda, and T. Yoshida, Journal of Cosmology and Astroparticle Physics 2020, 011 (2020).
  * Cherry  _et al._ (2020) J. F. Cherry, G. M. Fuller, S. Horiuchi, K. Kotake, T. Takiwaki, and T. Fischer, Physical Review D 102, 023022 (2020).
  * Tamborra and Shalgar (2020) I. Tamborra and S. Shalgar, Annual Review of Nuclear and Particle Science 71 (2020).
  * Shalgar and Tamborra (2021) S. Shalgar and I. Tamborra, Physical Review D 103, 063002 (2021).
  * Johns (2021) L. Johns, arXiv preprint arXiv:2104.11369 (2021).
  * Sasaki and Takiwaki (2021) H. Sasaki and T. Takiwaki, arXiv preprint arXiv:2109.14011 (2021).
  * Martin  _et al._ (2021) J. D. Martin, J. Carlson, V. Cirigliano, and H. Duan, Physical Review D 103, 063001 (2021).
  * Shalgar and Tamborra (2019) S. Shalgar and I. Tamborra, The Astrophysical Journal 883, 80 (2019).
  * Tarantola (2005) A. Tarantola,  _Inverse problem theory and methods for model parameter estimation_ (SIAM, 2005).
  * Kimura (2002) R. Kimura, Journal of Wind Engineering and Industrial Aerodynamics 90, 1403 (2002).
  * Kalnay (2003) E. Kalnay,  _Atmospheric modeling, data assimilation and predictability_ (Cambridge university press, 2003).
  * Evensen (2009) G. Evensen,  _Data assimilation: the ensemble Kalman filter_ (Springer Science & Business Media, 2009).
  * Betts (2010) J. T. Betts,  _Practical methods for optimal control and estimation using nonlinear programming_ , Vol. 19 (Siam, 2010).
  * Whartenby  _et al._ (2013) W. G. Whartenby, J. C. Quinn, and H. D. Abarbanel, Monthly Weather Review 141, 2502 (2013).
  * An  _et al._ (2017) Z. An, D. Rey, J. Ye, and H. D. Abarbanel, Nonlinear Processes in Geophysics (Online) 24 (2017).
  * Schiff (2009) S. J. Schiff, in  _2009 Annual International Conference of the IEEE Engineering in Medicine and Biology Society_ (IEEE, 2009) pp. 3318–3321.
  * Toth  _et al._ (2011) B. A. Toth, M. Kostuk, C. D. Meliza, D. Margoliash, and H. D. Abarbanel, Biological cybernetics 105, 217 (2011).
  * Kostuk  _et al._ (2012) M. Kostuk, B. A. Toth, C. D. Meliza, D. Margoliash, and H. D. Abarbanel, Biological cybernetics 106, 155 (2012).
  * Hamilton  _et al._ (2013) F. Hamilton, T. Berry, N. Peixoto, and T. Sauer, Physical Review E 88, 052715 (2013).
  * Meliza  _et al._ (2014) C. D. Meliza, M. Kostuk, H. Huang, A. Nogaret, D. Margoliash, and H. D. Abarbanel, Biological cybernetics 108, 495 (2014).
  * Nogaret  _et al._ (2016) A. Nogaret, C. D. Meliza, D. Margoliash, and H. D. Abarbanel, Scientific reports 6, 1 (2016).
  * Armstrong (2020) E. Armstrong, Physical Review E 101, 012415 (2020).
  * Djorgovski  _et al._ (2006) S. G. Djorgovski, C. Donalek, A. Mahabal, R. Williams, A. J. Drake, M. J. Graham, and E. Glikman, in  _18th International Conference on Pattern Recognition (ICPR’06)_ , Vol. 1 (2006) pp. 856–863.
  * Madhusudhan (2018) N. Madhusudhan, “Atmospheric retrieval of exoplanets,” in  _Handbook of Exoplanets_, edited by H. J. Deeg and J. A. Belmonte (Springer International Publishing, Cham, 2018) pp. 2153–2182.
  * Kitiashvili and Kosovichev (2008) I. Kitiashvili and A. G. Kosovichev, The Astrophysical Journal 688, L49 (2008).
  * Armstrong  _et al._ (2017) E. Armstrong, A. V. Patwardhan, L. Johns, C. T. Kishimoto, H. D. Abarbanel, and G. M. Fuller, Physical Review D 96, 083008 (2017).
  * Armstrong  _et al._ (2020) E. Armstrong, A. V. Patwardhan, E. Rrapaj, S. F. Ardizi, and G. M. Fuller, Physical Review D 102, 043013 (2020).
  * Rrapaj  _et al._ (2021) E. Rrapaj, A. V. Patwardhan, E. Armstrong, and G. M. Fuller, Physical Review D 103, 043006 (2021).
  * Volpe (2015) C. Volpe, International Journal of Modern Physics E 24, 1541009 (2015).
  * Vlasenko  _et al._ (2014) A. Vlasenko, G. M. Fuller, and V. Cirigliano, Physical Review D 89, 105004 (2014).
  * Abarbanel (2013) H. Abarbanel, AMC 10, 12 (2013).
  * Wächter (2009) A. Wächter, in  _Dagstuhl Seminar Proceedings_ (Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2009).
  * (49) minAone interface with Interior-point Optimizer, https://github.com/yejingxin/minAone, accessed: 2021-10-13.
  * Ye  _et al._ (2015) J. Ye, D. Rey, N. Kadakia, M. Eldridge, U. I. Morone, P. Rozdeba, H. D. Abarbanel, and J. C. Quinn, Physical Review E 92, 052901 (2015).