An optimization-based approach to calculating neutrino flavor evolution
Eve Armstrong, Amol V. Patwardhan, Lucas Johns, Chad T. Kishimoto, Henry D.I. Abarbanel, George M. Fuller
AAn optimization-based approach to calculating neutrino flavor evolution
Eve Armstrong, ∗ Amol V. Patwardhan,
2, 3, † Lucas Johns,
2, 3, ‡ Chad T.Kishimoto,
4, 3, § Henry D. I. Abarbanel,
2, 5, ¶ and George M. Fuller
2, 3, ∗∗ BioCircuits Institute, University of California, San Diego, La Jolla, CA 92093-0328, USA Department of Physics, University of California, San Diego, La Jolla, CA 92093-0319, USA Center for Astrophysics and Space Sciences, University of California, San Diego, La Jolla, CA 92093-0424, USA Department of Physics and Biophysics, University of San Diego, San Diego, CA 92110, USA Marine Physical Laboratory (Scripps Institution of Oceanography),University of California, San Diego, La Jolla, CA 92093-0210, USA (Dated: March 8, 2018)We assess the utility of an optimization-based data assimilation (D.A.) technique for treating theproblem of nonlinear neutrino flavor transformation in core collapse supernovae. D.A. uses measure-ments obtained from a physical system to estimate the state variable evolution and parameter valuesof the associated model. Formulated as an optimization procedure, D.A. can offer an integration-blind approach to predicting model evolution, which offers an advantage for models that thwartsolution via traditional numerical integration techniques. Further, D.A. performs most optimallyfor models whose equations of motion are nonlinearly coupled. In this exploratory work, we considera simple steady-state model with two mono-energetic neutrino beams coherently interacting witheach other and a background medium. As this model can be solved via numerical integration, wehave an independent consistency check for D.A. solutions. We find that the procedure can capturekey features of flavor evolution over the entire trajectory, even given measurements of neutrino flavoronly at the endpoint, and with an assumed known initial flavor distribution. Further, the procedurepermits an examination of the sensitivity of flavor evolution to estimates of unknown model param-eters, locates degeneracies in parameter space, and can identify the specific measurements requiredto break those degeneracies.
PACS numbers: TBA
I. INTRODUCTION
We assess the efficacy of a data assimilation (D.A.)technique for constraining neutrino flavor evolution his-tories inside core collapse supernovae. The specific tech-nique of interest in this paper is an integration-blind pro-cedure, which offers advantages for problems that thwarttraditional numerical integration techniques. The proce-dure is crafted to efficiently find solutions for an extrem-ized action from a sparse set of measurements [1]. Im-portantly, the technique of transporting information fromthe measured quantities through the complete model tothe un measured quantities requires that the model’s dif-ferential equations be coupled. For this reason, D.A.lends itself particularly well to the exploration of non-linear systems. Notably, collective neutrino oscillationphenomena in astrophysical environments are essentiallynonlinear. This feature begs the question of whether aD.A. approach to solving these problems is feasible. Herewe investigate this issue of feasibility in the context of asimple toy model that captures nonlinearity.Neutrino flavor evolution in compact object environ- ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] ¶ [email protected] ∗∗ [email protected] ments is a vexing and unsolved problem [2–52]. It isinherently nonlinear and fraught with difficulties. Thesedifficulties are exacerbated by the limitations of our un-derstanding of key supernova physics, for example: theequation of state and weak interaction properties of nu-clei and nuclear matter in hot and dense conditions.Even accounting for the inherent uncertainties in su-pernova microphysics, obtaining convincing numericalsimulations of supernova neutrino flavor histories insidethe supernova remains problematic. In just one exam-ple, inelastic neutrino back-scattering could contributeto flavor evolution in some regimes [27, 34, 53], givingrise to the “neutrino halo” problem. The neutrino haloeffect changes the flavor evolution problem from an ini-tial value problem with neutrino fluxes and flavor contentspecified at the edge of the proto-neutron star (the “neu-trino sphere”), to something more akin to a boundaryvalue problem, where flavor phase information is prop-agating both outward and inward . The computationaldifficulties endemic to neutrino direction-changing scat-tering represent just one of the many challenges we facein transitioning from the usual coherent “index of re-fraction” treatment of neutrino flavor evolution to a fullquantum kinetic approach [54–70]. Other outstandingproblems in supernova neutrino flavor physics includefast flavor conversion [48–51] and spatial and temporalinstabilities [26, 36, 37, 40, 52]. For background on neu-trino physics in massive stars and supernovae, see Ap-pendix D.Against this backdrop of uncertainty in theoretical cal- a r X i v : . [ a s t r o - ph . H E ] A ug culations, there is an effort to configure water or iceˇCerenkov (for example, HyperK, IceCube), liquid argon(for example, DUNE), and liquid scintillator detectors(for example, JUNO), to be able to capture the neutrinosignal from a Galactic core collapse supernova [71–78].The potential for such a detection to provide a probeof beyond-standard-model physics in the neutrino sectorand insight into the nuclear equation of state, and a hostof astrophysical issues like neutrino heating and nucle-osynthesis, is alluring. Consequently, the stakes are highwhen it comes to gaining confidence in modeling nonlin-ear neutrino flavor conversion. Given this context, it isworth exploring new techniques.The extraction of information from measurements is ageneral procedure known in the geosciences as data as-similation (D.A.) [79–82]. D.A. is used commonly in fluiddynamics [83, 84] and more recently in neuroscience [85–88], where the available data from a physical system aresparse and the corresponding model consists of degrees offreedom coupled in a nonlinear manner. The aim of D.A.is to incorporate information contained in measurementsdirectly into a model, to estimate unknown parametersand the dynamics of the model state variables - bothmeasured and unmeasured. The test of a successful esti-mation is the ability of the completed model to predictthe system state outside the times or locations at whichthe measurements were obtained.D.A. has two key advantages over numerical integra-tion. First, when cast within the framework of opti-mization, it can be written as an integration-blind pro-cedure. An integration-blind approach may be amenableto boundary-value problems for which solution via for-ward integration is unrealistic. Second, the procedurecan systematically and efficiently identify the existenceof degenerate solutions, and the specific measurementsthat are required to break degeneracy. For these reasons,we considered D.A. worth exploring as a possible alter-native attack to the standard initial-value treatment ofneutrino flavor evolution employed in, for example, the“bulb” model [2].We examine the potential utility of D.A. in the contextof collective neutrino oscillations by applying it to a sim-ple model that can be solved via numerical integration—and thus where there exists an independent consistencycheck for D.A. solutions. The model describes the flavorevolution of two mono-energetic neutrino beams emanat-ing from a supernova event. The measurements used are:the flavor content of each neutrino beam at some final ra-dius r = R at which a detector might be placed, given anassumed known initial flavor distribution at the surfaceof the “neutrino sphere” ( r = 0). We seek to determinewhether the sparse measurements suffice to yield the fla-vor evolution history over the interim distance, and toestimate unknown model parameter values: namely, thestrength of neutrino coupling to matter and the strengthof neutrino-neutrino coupling.In this first exploratory paper, three main resultsemerge. First, over repeated trials we obtain realistic overall flavor evolution history for both neutrinos, evengiven the extreme sparsity of measurements. We shalldescribe exceptions to this finding in Section V. Second,we find that generally the given measurements are insuf-ficient to distinguish among multiple parameter solutionsets. Third, and as a consequence of the first and sec-ond findings, we gain insight regarding the sensitivity offlavor evolution to these parameter values.This paper proceeds as follows. • Section II,
Inverse Problems and Optimization , ex-plains the general framework for data assimilationvia optimization, and it states the specific objec-tive function and method of evaluation used in thispaper. • Section III,
Model , describes the specific model usedin this paper: a simplified version of the dynam-ics of neutrino flavor evolution that ensue from thesurface of a core collapse supernova event. • Section IV,
The Experiments , explains the full pro-cedure for simulated experiments given to an opti-mization algorithm, and our physical rationale forthe experimental designs. • Section V,
Results , describes the solutions. • In Sec. VI,
Discussion , we comment on the impli-cations of the results with respect to more realisticproblems in neutrino astrophysics, and we describeimmediate future work. • Section VII,
Summary , contains concluding re-marks. • Appendix A gives the path-integral derivation ofthe objective function used in this paper. • Appendix B gives details of the optimization pro-cedure. • Appendix C gives an overview of how to interpretour simplistic model in terms of a constant-entropyenvelope model. • Appendix D gives background on relevant neutrinoastrophysics.
II. INVERSE PROBLEMS ANDOPTIMIZATIONA. General framework
D.A. is a procedure whereby information in measure-ments is used to complete a model of the system fromwhich the measurements were obtained (see Ref. [89] foran introduction to this “inverse problem” formulation).The model F is written as a set of D ordinary differentialequations that evolve in affine parameter r as:d x a ( r )d r = F a ( x ( r ) , p ); a = 1 , , . . . , D, where the components x a of the vector x are the modelstate variables. The affine parameterization r may be,for example, time or distance. The unknown model pa-rameters to be estimated are contained in p ; note thatthe model evolution depends on p .A subset L of the D state variables is associated withmeasured quantities. One seeks to estimate the p un-known parameters and the evolution of all state variablesgiven the provided measurements, and to then use thoseestimates to predict the model evolution in regions wherethere exist no measurements. The prediction phase is thetest of estimation quality.In the simulated experiments described in this paper,we integrate forward the equations of motion from aninitial known state at r = 0 (the surface of the neu-trino sphere), to obtain a “measurement” of L vari-ables at some detector location R . The question forD.A. is whether sufficient information can be propagatedthrough the coupled equations - from measured to un-measured variables - such that the parameters that hadgenerated the measurements can be inferred. As noted,the test of a successful parameter estimation is its abil-ity to predict the evolution of all D state variables in the interim r ∈ (0 , R ) during which no measurements areobtained. B. Specific optimization formulation used in thispaper
One method commonly employed to solve an inverseproblem is optimization. Optimization is the process offinding the extremum of a function, called the “cost (orobjective or penalty) function”. The cost function usedin this paper is motivated from a path-integral-like for-mulation of D.A., and for this reason we nickname it an“action”.In constructing the action that will be used to yield pa-rameter estimates, we consider three factors that will dic-tate those estimates: 1) measurements obtained from thephysical system of interest, 2) the dynamics of the modeldescribing that system, and 3) additional equality con-straints that are specific to the model formulation. Themeans by which we incorporate information from mea-surements, and the manner in which each of the abovethree factors is considered, may be best understood viaan examination of this cost function’s specific formula-tion. We do this now, and then proceed to describe eachterm in turn.The action A used in this paper is written as: A = R f ( N − D N − (cid:88) n ∈{ odd } D (cid:88) a =1 (cid:34)(cid:26) x a ( n + 2) − x a ( n ) − δr F a ( x ( n ) , p ) + 4 F a ( x ( n + 1) , p ) + F a ( x ( n + 2) , p )] (cid:27) + (cid:26) x a ( n + 1) −
12 ( x a ( n ) + x a ( n + 2)) − δr F a ( x ( n ) , p ) − F a ( x ( n + 2) , p )] (cid:27) (cid:35) + R m N meas (cid:88) j L (cid:88) l =1 ( y l ( j ) − x l ( j )) + k N (cid:88) n | g ( x ( n )) − | + k N (cid:88) n | g ( x ( n )) − | . (1)We seek the path X = { x (0) , . . . , x ( N ) , p } in statespace on which A attains a minimum value.The two squared terms in the first double sum in Eq. 1incorporate the model evolution of all D state variables x a . Of these, the first term in curly braces represents er-ror in the first derivative (with respect to r ) of the statevariables, whereas the second term corresponds to errorin the second derivative. These terms can be derivedfrom a consideration of Markov-chain transition prob-abilities. Here, the outer summation in n is taken overall odd-numbered grid points—discretized steps in r thatparameterize the model equations of motion. The step-size δr is defined as the distance between alternate gridpoints: δr ≡ r n +2 − r n . The inner summation in a is taken over all D state variables.The squared term in the second double sum governs thetransfer of information from measurements y l to states x l . It derives from the concept of conditional mutualinformation of probability theory. The y l are the mea-surements, and x l are the model variables correspondingto the measurements. Here, summation on j runs overthe set of all N meas discretized grid points where the mea-surements are made, which may in general be some subsetof all the model grid points. The summation in l is takenover the L measured quantities. R m and R f are inverse covariance matrices for the mea-surement and model errors, respectively. In this paper wetake the measurements to be mutually independent andthe state variables to be independent, rendering thesematrices diagonal. Additionally, we constrain R f to beuniform across all state variables, and likewise R m for allmeasurements. For our purposes, R m and R f are relativeweighting terms; the utility of relative weighting will bedescribed below in this section. For a short derivation ofthe first two terms, see Appendix A; for a full treatment,see Ref. [1].The third and fourth terms (with coefficients k ) areequality constraints, which were added to increase theefficiency of the search algorithm. These will be writtenout explicitly in Section III.The optimization is performed at all locations along apath simultaneously, so as not to impart greater impor-tance to a measurement at any particular location overanother. An integration-blind technique may lend itselfwell to problems that cannot be solved in a straightfor-ward manner via forward integration. An example ofsuch a problem — the “neutrino halo” — is discussed inAppendix D.To minimize A we employ a variational approach, viathe open-source Interior-point Optimizer (Ipopt) [90]; seeSection IV. C. Iterative reweighting of measurement andmodel contributions (“Annealing”)
The complete D.A. procedure involves an iteration thatis aimed to identify the parameter set corresponding tothe global minimum of the action. Local minima will rep-resent degenerate sets of parameter estimates that maywell fit the measurements but which are poor predic-tors of state variable evolution outside of the locationsat which the measurements were obtained. The rem-edy consists of recursively calculating A as the ratio ofthe model and measurement coefficients - R f and R m ,respectively - is gradually increased. Specifically: we de-fine R f = R f, α β , where α is a small number greaterthan 1, and β is increased from 0 in uniform increments.For all simulations performed in this paper, R f takes auniform value over all state variables x a . R m is uniformover all measured variables x l , and it is held fixed acrossiterations while R f is incremented. For an explanation ofwhy this iterative procedure - which we call “annealing” -aids in identifying the global minimum, see Appendix A. III. MODELA. Toy model motivation and scheme
A complete description of neutrino flavor evolution inrealistic astrophysical environments, such as supernovae,involves complications imposed by geometry, multi-angleeffects, realistic emission spectra, and non-forward scat-tering, among other effects. Turning D.A. into a use-ful tool to simulate the dynamics in such a model, while matching the fidelity and sophistication of cur-rent forward-integration codes like “BULB”, is a daunt-ing task. We shall not attempt this here. Rather, ourmotivation in this paper is to take a first step toward as-sessing the efficacy of D.A. in treating the astrophysicalneutrino flavor transformation problem. This first steprequires the use of a vastly simplified model.The model we craft possesses two key features. First,it is nonlinear - a key aspect of the physics that gives riseto collective neutrino flavor evolution in these environ-ments. Notably, D.A. is particularly useful for estimat-ing model evolution and parameter values in nonlinearmodels where only a subset of the state variables canbe accessed experimentally. Such is the case for neu-trino flavor evolution. A second key point is that themodel is sufficiently simple to be solvable via traditionalforward-integration techniques. This feature enables anindependent consistency check for D.A. solutions.For our model, we consider a scenario in which two mo-noenergetic neutrino beams with different energies inter-act with each other and with a background consisting ofnuclei, free nucleons, and electrons. The densities of thebackground particles and of the neutrino beams them-selves are taken to dilute as some functions of a positioncoordinate r , which could be interpreted as the distancefrom the neutrino sphere in a supernova. In what fol-lows, we first discuss general two-flavor neutrino flavoroscillations, and we then adapt this formalism for ourparticular toy model. B. Two-flavor neutrino flavor evolution
Since ν µ and ν τ neutrino flavors experience identical in-teractions in the supernova environment, the three-flavorproblem can be reduced to a two-flavor mixing between ν e and a state ν x that is a particular superposition of ν µ and ν τ [91, 92]. The flavor state of neutrinos of energy E , as a function of position, can then be expressed usinga 2 × ρ E ( r ) = (cid:18) ρ ee,E ρ ex,E ρ xe,E ρ xx,E (cid:19) = (cid:18) | a e,E | a ∗ e,E a x,E a e,E a ∗ x,E | a x,E | (cid:19) . (2)Here the last matrix representation of the density op-erator is for the special case where the neutrinos are inpure states, with a e,E and a x,E being the respective fla-vor amplitudes. The quantum kinetic equation (QKE)governing the evolution of the general density operator ρ E has the form [16, 54–70, 93, 94] i dρ E ( r ) dr = [ H E ( r ) , ρ E ( r )] + i C E ( r ) . (3)Here we are assuming that the neutrino density ma-trix elements and potentials carry no explicit time-dependence; that is: they may vary only as functionsof position along the neutrino trajectory—a steady-statesolution. H E ( r ) is the Hamiltonian driving coherent fla-vor evolution. The last term on the right side, C E ( r ),captures the effects of collisions. Neglecting collisions,which may be justified in some supernova regions andepochs [53, 56], results in the coherent limit in whichneutrino flavor evolution is Schr¨odinger-like: i d ρ E ( r )d r = [ H E ( r ) , ρ E ( r )] . (4)Here, the right side is trace-conserving, implying unitaryevolution. Equation (4) can also be cast in the form of astandard path-integral extremization problem [24].
1. Spin basis and polarization vectors
Equation (4) can be recast in terms of Bloch-vectors P E and H E by decomposing the density matrices ρ E andHamiltonians H E , respectively, into a basis of Pauli spinmatrices (for details see Ref. [66]): ρ E = 12 ( P E, I + P E · σ )= 12 (cid:18) P E, + P E,z P E,x − iP E,y P E,x + iP E,y P E, − P E,z (cid:19) , (5) H E = 12 ( H E, I + H E · σ )= 12 (cid:18) H E, + H E,z H E,x − iH E,y H E,x + iH E,y H E, − H E,z (cid:19) . (6)We refer to the quantities P E as “Polarization vec-tors”, whereas the vectors H E will inherit the name“Hamiltonians”. Note that the subscripts x, y, z on thevector components above do not refer to spatial coor-dinates, but rather to directions in this SU (2) “Bloch-space”. The advantage of Bloch-vector decomposition istwo-fold: (1) the dynamical variables of the system; thatis: the components of P E , are now real numbers, unlikethe complex amplitudes in the density matrices, and (2)the geometric representation makes for easier visualiza-tion of the often complex underlying dynamics. To illus-trate the second point, we write the evolution equationin terms of these Bloch vectors:d P E, d r = 0 , and d P E d r = H E ( r ) × P E ( r ) . (7)The first equation is simply a restatement of tracepreservation, and in fact, for a normalized density ma-trix, P E, = 1. The second equation, on the other hand,resembles Larmor precession of a magnetic moment, withthe Hamiltonian in this case playing the role of a mag-netic field. Note that P E,z represents the probabilityof a neutrino being detected as a ν e over ν x ; that is: P E,z = | a e,E | − | a x,E | . Or, if one were to think interms of a population of neutrinos: P E,z ( r ) = n ν e ,E ( r ) − n ν x ,E ( r ) n ν,E ( r ) , (8) where n ν,E ( r ) is the number density of neutrinos of en-ergy E at a position r , and n ν α ,E ( r ) ≡ n ν,E ( r ) | a α,E ( r ) | is the “expected” number density of neutrinos in the fla-vor α at that energy and position.
2. The Hamiltonian and equations of motion for neutrinoforward scattering
Having set up the evolution equations, let us now de-scribe the specific Hamiltonians that drive flavor evo-lution in the coherent limit—first in matrix form andthen in the Pauli spin representation. The Hamilto-nian H E consists of three contributing terms: H E ( r ) = H vac ,E + H m ( r ) + H νν ( r ). The first of these terms drivesflavor oscillations in vacuum: H vac ,E = ∆2 (cid:18) − cos 2 θ sin 2 θ sin 2 θ cos 2 θ (cid:19) , (9)where θ is the mixing angle in vacuum, describing theunitary transformation between the weak interaction (fla-vor) eignestates, and the energy (mass) eignestates. Also,∆ ≡ δm / E , where δm is the mass-squared split-ting between the two energy eigenstates. The other twocontributions arise from neutrino forward-scattering onbackground matter particles ( H m , the “matter Hamilto-nian”) and other neutrinos ( H νν , the “neutrino-neutrinoHamiltonian”), respectively. In the scenario describedabove, they assume the forms: H m = V ( r )2 (cid:18) − (cid:19) , (10) H νν = (cid:88) E µ E ( r ) ρ E ( r ) . (11)Here we have already subtracted the trace from thevacuum and matter Hamiltonian, since it has no bearingon the flavor evolution. In terms of the baryon numberdensity n B ( r ), the electron fraction Y e ( r ), and the neu-trino number density n ν,E ( r ), the potentials are V ( r ) = √ G F n B ( r ) Y e ( r ) ,µ E ( r ) = √ G F α ( r ) n ν,E ( r ) , (12)where G F is the Fermi constant and α ( r ) ≡ − cos ψ ( r ) isa factor that weights the neutrino-neutrino coupling ac-cording to the intersection angle ψ ( r ) between the twoneutrino streams. Our choice of particular functionalforms for V ( r ) and µ E ( r ) is stated and explained later inthis section (III B 4) and in Appendix C. The neutrino-neutrino Hamiltonian H νν ensures that the evolutionequations are nonlinear, since the Hamiltonians drivingthe evolution of the density matrices depend on the den-sity matrices themselves; moreover, the evolution histo-ries of the two neutrino populations are now coupled toone another.Gathering the above Hamiltonians and expressingthem in the Pauli basis, one obtains the complete setof dynamical equations for the two neutrino beams:d P ,x d r = (∆ cos 2 θ − V ( r )) P ,y + µ ( r )( P ,y P ,z − P ,z P ,y ) , d P ,y d r = − (∆ cos 2 θ − V ( r )) P ,x − ∆ sin 2 θP ,z + µ ( r )( P ,z P ,x − P ,x P z, ) , d P ,z d r = ∆ sin 2 θ P ,y + µ ( r )( P ,x P ,y − P ,y P ,x ) , (13)d P ,x d r = (∆ cos 2 θ − V ( r )) P ,y + µ ( r )( P ,y P ,z − P ,z P ,y ) , d P ,y d r = − (∆ cos 2 θ − V ( r )) P ,x − ∆ sin 2 θP ,z + µ ( r )( P ,z P ,x − P ,x P z, ) , d P ,z d r = ∆ sin 2 θ P ,y + µ ( r )( P ,x P ,y − P ,y P ,x ) , (14)where for simplicity we have assumed equal neutrinonumber densities at both energies ( n ν,E ( r ) = n ν,E ( r )),so that µ E ( r ) = µ E ( r ) = µ ( r ). For brevity, we haveused P and P in place of P E and P E .
3. Physics of the model: MSW resonance and collectiveeffects
In principle, the various Hamiltonians driving neutrinoflavor evolution can – in the adiabatic limit – be com-bined and expressed in the form of effective in-mediumoscillation parameters: H vac + H m ( r ) + H νν ( r ) ≡ ∆ m ( r )2 (cid:18) − cos 2 θ m ( r ) sin 2 θ m ( r )sin 2 θ m ( r ) cos 2 θ m ( r ) (cid:19) , (15)where ∆ m ( r ) = (cid:113) (∆ cos 2 θ − V eff ( r )) + ∆ sin θ andsin θ m ( r ) = (∆ sin θ ) / ∆ m ( r ) represent the effectivein-medium mass-squared difference and mixing angle, re-spectively. Here, V eff is taken to represent the effectivematter + collective potential experienced by a neutrino.At V eff ( r ) = ∆ cos 2 θ , the in-medium mixing angle θ m achieves its maximal value of π/ H m dominates, so that the heavier mass eigenstate essentially aligns with the ν e fla-vor state. At large radii, however, H vac takes over, andfor small mixing angles, this means the heavier masseigenstate aligning more closely with ν x . If this transi-tion is sufficiently adiabatic, a neutrino initially emittedas ν e undergoes near-complete conversion into ν x priorto its detection.In the numerical calculations discussed below wechoose neutrino energy ratios and matter potentials thatcan encompass highly adiabatic neutrino flavor trans-formation, so that neutrinos stay in instantaneous masseigenstates. Knowing what the flavor states of neutrinosare at the beginning of our calculations, that is: at theneutrino sphere, we can then determine the flavor statesat the end, without knowing the precise details of theintervening matter density profile. It is because of thisreason that adiabatic neutrino flavor evolution presentsa fundamental problem in interpreting a detected corecollapse neutrino signature: possible degeneracy of neu-trino flavor histories. That is, any number of smoothmatter density profiles, each transited by neutrinos adi-abatically, will facilitate conversion of an initial ν e intoa ν x , or vice versa . A key objective of this study is toascertain whether optimization techniques can map outdegeneracies. In Sec. VI, we suggest that introducing ad-ditional complexities in our model, including sharp fea-tures in the matter potential (such as shocks) that wouldengineer non-adiabaticity, can help break such degenera-cies [97–100].Introducing neutrino-neutrino coupling into this pic-ture gives rise to an array of nonlinear flavor-transformation phenomena. Nonlinearity can manifest asvarious modes of collective neutrino flavor oscillation—see Ref. [12] for a review. In these collective modes signif-icant fractions of the neutrinos in a range of energies andlocations may undergo simultaneous, sometimes synchro-nized coherent flavor oscillations. In essence, neutrino-neutrino forward scattering serves to “inform” a neutrinoabout the flavor states of others, and the nonlinear na-ture of the interactions guides neutrino flavor states intolock-step coherence. Determining the locations in radiusand energy of the transition in and out of such collec-tive modes, or whether they even occur at all, will be animportant objective for core collapse supernova neutrinoburst detection.In a practical sense, collectivity engendered by nonlin-ear neutrino-neutrino forward scattering potentials mayadd to the possible degeneracy in neutrino flavor his-tories, or it may tend to narrow the possibilities. Touse optimization techniques to explore this question, wewill now present simple functional forms for matter andneutrino-neutrino potentials.
4. Choice of the matter potential and the neutrino couplingterm
The above formulation of 2 × × V ( r ) is typically written in termsof the baryon density n B ( r ) and the electron fraction Y e ( r ). For simplicity we combine the two dependencesand describe V ( r ) using a single power-law V ( r ) = Cr , (16)where all constants, including the weak coupling G F , andphysical parameters such as the neutrino sphere radius,and n B and Y e at the neutrino sphere, have been ab-sorbed into the dimensionful constant C , which we treatas a parameter to be determined by the data assimilationprocedure.We also choose a simplified structure for the neutrino-neutrino coupling term µ ( r ). Here, the dependences ofthe neutrino number density n ν,E ( r ), and the effect ofthe intersection angle α ( r ) ≡ − cos ψ ( r ) are bundledtogether into a single power-law µ ( r ) = Qr , (17)where, as with the matter potential, all constant param-eters are absorbed into Q , which we will treat as a sin-gle parameter to be determined by the data assimilationprocedure. These functional forms are adopted as coarsemock-ups of the matter and neutrino densities surround-ing the neutrino sphere and are not meant to emulaterealistic profiles to an accurate degree. Nevertheless, thephysical motivation behind the choice of these functionalforms is discussed in Appendix C.The challenge posed to the data-assimilation machin-ery is to estimate the constants C and Q as well as theflavor-space trajectories P ( r ) and P ( r ) of the two neu-trino beams as they propagate outward from the neutrinosphere (radius r = 0) towards some radius R . We imag-ine that a detector sits at R . For this exploratory D.A.study, we have chosen energies that allow us to examinehow the procedure operates over different resonance lo-cations relative to the detector location. Our motivationis to probe the utility of eventually adding constraints onphysics within the envelope. The inputs to the D.A. ma-chinery are: (1) the model equations of motion and (2) In practice, we set the dependences as 1 / ( r + 0 . , to avoidinfinities at r = 0, where 0 . the measurements P z of each neutrino at r = R (givena known initial state at r = 0), which are processedthrough the action-minimization procedure detailed inthe previous section. IV. THE EXPERIMENTSA. The model-specific optimization procedure
Given: 1) the model embodied in Eqs. 13 and 14, andwith unknown parameters C (the weight of the matterpotential) and Q (the weight of the coupling potential),2) measurements of the model state variables P ,z and P ,z at r = R , and 3) an assumed initial known flavorstate ( P z ) of each neutrino at r = 0, we seek to identifythe path X = { P (0) , P (0) , . . . , P ( R ) , P ( R ) , C, Q } instate space such that the cost function of Eq. 1 attainsa minimum value. Within our model formulation, themeasurements y l are the P z components of both neutri-nos, and the components of x are the P x , P y , and P z values of both neutrinos. (Obviously, the only poten-tial “measurement” of supernova neutrinos would be ina terrestrial detector, and would correspond to an energyspectrum; see Sec. VI.)The two equality constraints (with coefficients k ) aredesigned to improve the efficiency of the search algo-rithm. The algorithm does not recognize relations amongnon-independent state variables, but rather considerseach independently. Because the model described byEquations 13 and 14 implicitly imposes P = constant,the model is overdetermined in the cartesian coordinatesystem. To minimize the computational expense, weadded these equality constraints to strictly impose uni-tarity at the start of the annealing procedure ( R m (cid:29) R f ), in which regime the model weight may not yet besufficiently strong for its implicit requirement to be wellrespected. The functions g and g are: g ( x ( n )) = P ,x + P ,y + P ,z g ( x ( n )) = P ,x + P ,y + P ,z ;the value of coefficient k in Equation 1 was taken to be1. B. The experimental designs
We designed two sets of experiments for D.A., wherethe sets are distinguished by the ratio of the neutrino en-ergies. In both cases the first neutrino (Neutrino 1) expe-riences the MSW resonance at a radius r = 1 . Q = 0). (We express both the location r and matter and neutrino potential coefficients C and Q as dimensionless quantities. We can provide dimensions,for example in cm and MeV cm , respectively, throughEqs. 18 and C2; see Appendix C). This requirement setsthe value of the matter coefficient C . In the first set ofexperiments, the ratio of neutrino energies E ν /E ν =2 .
5; in the second set of experiments, E ν / E ν = 0 . Q = 0,1, 100, and 1000. These choices were made to permit anexamination of whether the quality of results is sensitiveto neutrino coupling strength or energy ratio.For all experiments, we assumed the neutrinos to bein pure ν e flavor at r = 0, with a corresponding P z valueof +1. Here we make a note regarding our treatment ofthis bound constraint at r = 0. While the only actualmeasurement in this experiment occurs at r = R , for thepurposes of the D.A. procedure we treated the knowninitial state (at r = 0) as a “measurement” as well. Wethen added a one per cent uncertainty to that “measure-ment”. In this way, the initial known distribution wastreated by the algorithm as a bound constraint of finiterigidity. To be reasonable, we also added one-percentnoise to the measurement of P z at r = R .Finally, for all eight experiments, a search was per-formed 20 times, each beginning from a randomly-chosenlocation on the state space surface (for a description ofthe discretized space, see Appendix B, Section B 1).The simulated data were generated by forward-integrating Eqs. 13 and 14 out to r = R , via a fourth-order adaptive Runge-Kutta scheme “odeINT”, an open-access Python integrator. Ipopt uses a Simpson’s Rulemethod of finite differences to discretize the state space,a Newton’s Method to search, and a barrier method toimpose user-defined bounds that are placed upon thesearches. The integration step for the simulated data,and the discretization step, were each ∆ r = 0 . r = 0 to R ) consist of 20,000 points. The annealingprocedure took β from 0 to 30 in increments of 1, and R f, = 0 . j has two terms: r = 0 and R , and the sum on the mea-sured state variables l has two terms: P ,z and P ,z . Forthe model error, the sum on location n has 20,000 terms,and the sum on all state variables a has six terms: P ,x , P ,y , P ,z , P ,x , P ,y , and P ,z . We performed the exper-iments for various values of R m between 1 and 10,000.A link to the Python codes that we used to interfacewith Ipopt is provided in Appendix B, Sec. B 2. C. Rationale for experiments in light ofastrophysical considerations
The specific experimental designs we adopt—namely,the two energy ratios E ν /E ν and various values of cou-pling strength Q —were crafted to be analogies to inter-esting physical scenarios. In Appendix C we illustratethese analogies by providing supernova envelope exam-ples in which to “embed” the chosen neutrino energy ra-tios. The examples with a neutrino energy ratio of 2 . R , is well outside the su-pernova envelope, and well beyond any MSW resonances.With completely adiabatic flavor evolution this will corre-spond to the completely degenerate case, where the initialand final neutrino flavor states are essentially predeter-mined, and the flavor state history between these pointsis not uniquely determined. Consequently, we might ex-pect the optimization algorithm to fail to converge con-sistently to a single C -value in the case where Q = 0.Non-adiabatic flavor evolution, for example, because ofdensity ledges or shocks, would be expected to break thisdegeneracy.In the examples with a neutrino energy ratio of0 .
01, the final location R was held unchanged from the E ν /E ν = 2 . r = R . In particular, in these examples, theresonance location of Neutrino 2 is shifted beyond thefinal location r = R , in the limit where Q (cid:28) C . Phys-ically, this corresponds to a scenario in which the finallocation R is inside the supernova envelope.Our purpose behind this choice is twofold. First, weaimed to assess whether optimization techniques can cap-ture with fidelity the neutrino flavor evolution, if infor-mation is specified about the flavor content at certainlocations within the envelope. This information mightbe important for calculating nucleosynthesis or neutrinoheating, since both of these issues hang mostly on the ν e and ¯ ν e content of the local neutrino fluxes. For ex-ample, alpha-rich freeze-out in the post-accretion phasewill lead to overproduction of nuclides such as Zr, un-less the electron fraction is larger than about 0.48 [101].This process may place constraints on the electron neu-trino/antineutrino energy spectra and fluxes in the post-accretion epoch. There also exist speculative theoreticalmodels about the possible production of r -process nu-clides in core-collapse supernova environments, for ex-ample, via neutrino-spallation-induced liberation of neu-trons in the helium layer [102, 103].Our second motivation for choosing a value of R insidethe supernova envelope was to facilitate a direct com-parison between different energy ratios, in order to ex-amine the sensitivity of the D.A. procedure for differentfinal radii relative to the resonance locations. We includethese results here because our aim throughout this paperis not to capture physical realism, but rather to examinethe robustness of the D.A. machinery over various modelregimes.Finally, and before giving examples of specific super-nova conditions (Appendix C), it will prove useful tonote that the MSW resonance condition, δm cos 2 θ =2 E ν V eff ( r ), where V eff ( r ) is the flavor-diagonal potentialfrom background matter and neutrinos at location r , letsus determine the resonance location for a given neutrinoenergy. For purely matter-driven ( Q = 0) flavor evolu-tion this means that the location of the resonance is r res = (cid:18) E ν δm cos 2 θ (cid:19) / C / , (18)and consequently the ratio of the resonance locations fortwo different neutrino energies, E ν and E ν , is indepen-dent of C in this case: r res , r res , = (cid:18) E ν E ν (cid:19) / . (19) V. RESULTSA. Key findings
Before presenting details of the results, we summarizekey findings: • For six out of the eight experiments (which were de-fined in Section IV), the measurements contain suf-ficient information to qualitatively capture overallflavor evolution through the MSW resonance, giventhe estimated parameter values. That is: the flavorevolution histories are consistent with the parame-ter estimates.
These six experiments correspond tomodel regimes in which the equations of motion forthe neutrinos are strongly coupled to each other, tothe matter background, or both. • For these six experiments, the precise location ofMSW transition is estimated correctly in roughly33% of trials. • For the other two (out of eight) experiments, theD.A. result failed to capture the model evolution.
These two cases correspond to low coupling of theequations of motion of ν to both matter and to ν . • There is insufficient information in the measure-ments to break the degeneracy in allowed sets ofparameter values C and Q.
This result is expected,as broad ranges of the parameters C and Q willleave flavor evolution completely adiabatic, therebymatching the neutrino flavor values imposed at theendpoints. The corresponding picture that emergesof the state space surface is not one riddled withclearly defined local minima of varying depth, butrather a wide, relatively flat basin. • The sensitivity of model evolution to parameter val-ues may depend on: 1) energy ratio E ν /E ν , 2) thevalue of neutrino-neutrino coupling strength Q , 3)the location of the detector relative to the resonancelocation, or 4) any combination of the above . Forexample, when the detector sits outside the regionof flavor transformation, the model is insensitive tothe values of C and Q . When the measurement is made prior (in location) to complete flavor trans-formation, however, a correlation emerges betweenestimates of C and Q . This finding suggests thatthe addition of physical constraints within the su-pernova envelope (in a more complicated model)could prove useful for degeneracy breaking. • The high-frequency, low-amplitude oscillations thatmodulate the overall transformation are fit poorly,in comparison to the fit of the overall transforma-tion.
This is due to the high sensitivity of thisaspect of the model to precise estimates of C and Q . B. Predicted flavor evolution and parameterestimates
1. Predicted flavor evolution histories
Examples of flavor evolution for the eight experimentsare depicted in Figs. 1 and 2, for neutrino energy ratios E ν / E ν = 2 . .
01, respectively. For each experi-ment, polarization vector components P x ( r ), P y ( r ), and P z ( r ) are shown at left for ν and at right for ν . Therows correspond to results for Q parameter values of 0,1, 100, and 1000, respectively. Predicted flavor evolutionby the D.A. procedure is shown in red, alongside evolu-tion curves in blue corresponding to the correct evolutionobtained by forward integration, for comparison. Thebest results, over all trials, were obtained by a choice ofmeasurement weight R m of 1 and annealing parameter β between 13 and 15.For six out of the eight experiments, we found themodel dynamics to be captured well. For two out ofthe eight experiments, flavor evolution was traced poorlyover all trials. We will first discuss the six relatively suc-cessful experiments, and then separately the final two.The six relative successes were: all four experiments forthe E ν /E ν = 2 . Q ), and thetwo experiments for the E ν /E ν = 0 .
01 high coupling( Q = 100 and 1000). For all six of these experiments, thecorresponding plots on Figs. 1 and 2 represent roughlyone-third of results over all trials. For the other roughly66 per cent of trials, the overall transformation historyhad the same qualitative appearance, but the precise lo-cation r res of resonant flavor conversion was matched lessprecisely to the model evolution. By “less precisely”, wemean roughly a discrepancy between model and predic-tion captured in the top left panel in Fig.2: P x , P y , and P z for ν in the case for Q = 0.As we will describe in Sec. VI, all of these results cap-ture the expected behavior of resonant transformationmodified by neutrino self-coupling; the offsets in r res overthe trials are due to different (degenerate) sets of param-eter estimates.Next we examine the two experiments for the E ν /E ν = 0 .
01 case, for low coupling: Q = 0 and0 Radius RadiusRadius RadiusRadius RadiusRadius Radius
Q = 1000Q = 100Q = 1Q = 0
Neutrino 1 Neutrino 2
FIG. 1. Flavor evolution histories for the four experiments with neutrino energy ratio E ν /E ν = 2 .
5. For each experiment,polarization vector components P x ( r ), P y ( r ), and P z ( r ) are shown at left for ν and at right for ν . Red curves representneutrino flavor evolution histories predicted by Ipopt, given P z measurements at the endpoints. The blue curves shown forcomparison are the model solutions obtained by forward integration of the equations of motion. The rows correspond to resultsfor Q parameter values of 0, 1, 100, and 1000, respectively. The true (model) value of C was 1304.5 for all four experiments.
11. Unlike results for the other six experiments, all ofthese D.A. solutions fail to match the features of thetrue evolution of ν , including the measurements at theendpoints. This result may be interpreted in terms of theefficiency of information flow among the state variables.As noted, if one has available as measurements only thosecorresponding to a subset of the model’s total numberof state variables, then in order to obtain informationregarding the unmeasured states, the equations betweenmeasured and unmeasured states must be coupled tosome significant degree. This interpretation is borne outby observations that generally, across physical modelsin other fields, D.A. tends to perform poorly when theequations of motion are not strongly coupled [104], as isthe case here: ν is not strongly coupled to ν , and it isalso far from its resonance location.
2. Parameter estimates
The parameter estimates corresponding to the flavorevolution histories of Figs. 1 and 2 are listed in Table I.We found that for each experiment, the estimates variedacross trials, with values of C and Q that spanned thepermitted search ranges for each parameter (not shown);the search ranges are specified in the caption of Table I.Note that the search ranges were different for each valueof Q chosen; see Sec. VI. In addition, for comments onthe lack of errors on these parameter estimates and onmethods to quantify these errors, see Sec. VI. TABLE I. Parameter estimates C and Q corresponding tothe solutions whose state variable evolution is displayed inFigs. 1 and 2, for E ν /E ν = 2 . . C : 500 to 1900, for allexperiments. Permitted search range for Q = 0: -10 to 10;for Q = 1: 0 to 10; for Q = 100: 0 to 200; for Q = 1000: 500to 1900. E ν /E ν = 2 . E ν /E ν = 0 . C (model) Q (model) C (DA) Q (DA) C (DA) Q (DA)1304.5 0 1457 0.3 563 10 ∗ ∗ value is a permitted search range bound To examine possible model sensitivity to parameterestimates in different model regimes, we explored C - Q space for all trials over all eight experiments. For thefour experiments with E ν /E ν = 2 .
5, no statistical trendemerged. Note that for this case, the detector sits at a lo-cation beyond the complete flavor conversion of both neu-trinos. For E ν /E ν = 0 .
01, however, a trend emergedat Q = 100 and 1000. Figure 3 illustrates this trend. As Q is increased from 100 to 1000, it clearly emerges thata correlation between the C and Q values is required by the estimates. Note that for this case, the detector sitsat a location at which the resonant flavor conversion of ν is not complete. See Sec. VI for comments.Finally, to test whether the low noise added to themeasurements was in part responsible for the high de-generacy of these estimates, we repeated all experimentswith zero noise in the measurements. Results were es-sentially unchanged. C. Tests to ascertain success of prediction, giventhe parameter estimates
The degeneracy of parameter estimates describedabove gave rise to an important question: Could the de-viation of the state prediction from precise MSW reso-nance be attributed to the particular estimates of C and Q ? In other words, we sought to ascertain whether - overall trials for the six out of eight rather successful experi-ments - the evolution history was traced correctly giventhe respective parameter estimates. To this end, we examined the value of A over all val-ues of the annealing parameter β . The top panel of Fig. 4shows the value of the action A over the annealing proce-dure corresponding to the first solution depicted in Fig. 1,where E ν /E ν = 2 . Q = 0. Thisis a purely matter-driven neutrino evolution case. The y -axis is the base-ten logarithm of A , and the x -axis isthe parameter β defined by: R f = R f, β ; R f, = 0 . R m = 1.This top panel of Fig. 4 is representative of the A -versus- β plots for all of the most-precisely-fit flavor evo-lution solutions presented in Figs. 1 and 2; that is, forroughly 33 per cent of trials in six out of the eight ex-periments. Specifically, the logarithm of A holds at aconstant value of -4.0 (up to machine precision) for allvalues of β .The bottom panel of Fig. 4 shows an overlaid plot of A -versus- β for all 20 trials corresponding to the first ex-periment (again, for E ν /E ν = 2 . Q = 0.) Anexamination of each of the 20 results individually re-vealed that the value log A = − . R m (cid:29) R f and R f (cid:29) R m ) on each plotwas correlated with the degree to which the correspond-ing flavor history solution matched the precise locationof transformation.As there occurs no evolution of action with increasing R f , we conclude: the dominant contribution to the ac-tion is the measurement error, and the model dynamicsare obeyed well. Further, given that A has the same min-imum value for all parameter sets, we infer: all parametersets - and corresponding evolution histories - equally-wellcapture the given measurements. Finally, we note two systematic features that we haveidentified for obtaining best-fit solutions to the flavor his-tories. Namely: all best fits, as exemplified in Figs. 1 and2
Neutrino 1 Neutrino 2
Q = 0Q = 1Q = 100Q = 1000
Radius RadiusRadius RadiusRadius RadiusRadius Radius
FIG. 2. Flavor evolution histories for the four experiments with neutrino energy ratio E ν /E ν = 0 .
01. For each experiment,polarization vector components P x ( r ), P y ( r ), and P z ( r ) are shown at left for ν and at right for ν . Red curves representneutrino flavor evolution histories predicted by Ipopt, given P z measurements at the endpoints. The blue curves shown forcomparison are the model solutions obtained by forward integration of the equations of motion. The rows correspond to resultsfor Q parameter values of 0, 1, 100, and 1000, respectively. The true (model) value of C was 1304.5 for all four experiments.Note the poor match to the flavor evolution of ν for the low-coupling cases ( Q = 0 and 1). See text for comments. (model)(model) FIG. 3. Relation between Ipopt-predicted values of param-eters C (matter coupling strength) and Q ( ν - ν couplingstrength), for E ν /E ν = 0 .
01, for: model Q = 100 ( top )and model Q = 1000 ( bottom ). Symbol representations areas follows. Blue x’s: estimates over all trials for all valuesof annealing parameter β ; solid green squares: estimates overall trials for the values of β that consistently yielded the bestfits to flavor evolution: β = 13 through 15; open black circle:the estimate corresponding to the plot of flavor evolution onFig. 2; solid red star: the true model value. A weak positivecorrelation between C and Q appears for Q = 100, and getsstronger at Q = 1000. The same plot (of C - Q space) for theexperiments in which E ν /E ν = 2 . Q (not shown). See Sec. VI forcomments.
2, correspond to values of β between 13 and 15 (out of thefull span of zero to 30), and R m = 1. The optimal valuesof these user-defined parameters vary across dynamicalmodels, and for any given model it is important to iden-tify a systematically reliable choice. Currently, this canonly be done via trial-and-error; see Appendix B, Sec-tion B 2 b. VI. DISCUSSION
There are many issues in nonlinear compact objectneutrino flavor evolution that are difficult to treat withstandard initial-value-problem integrations of the neu-
FIG. 4. The action A over the annealing procedure. The y -axis is log A , and the x -axis is the parameter β defined by: R f = R f, β ; R f, = 0 .
01 and R m = 1. Top : A -versus- β forthe first solution depicted in Fig. 1, where E ν /E ν = 2 . Q = 0. Specifically, log A holds constantat -4.0 to machine precision. This panel is representative ofthe roughly 33 per cent of solutions that fit flavor evolutionwell. Bottom : An overlaid plot of A -versus- β for all 20 solu-tions corresponding to the first experiment (for E ν /E ν = 2 . Q = 0.) The degree of scatter at the extremes ( β ∼ trino flavor quantum kinetic equations as described inSec. III. This concern constitutes our reason for consid-ering numerical path integral-inspired approaches to thisproblem. To this end, here we have investigated the effi-cacy with which a particular D.A. algorithm (1) capturesneutrino flavor evolution, and (2) identifies the sensitiv-ity of that evolution to parameter values. It is significantthat using one particular algorithm (Ipopt), we were ablecapture key features, like MSW resonance locations, ofa simple supernova envelope neutrino flavor oscillationmodel. Furthermore, our results do indeed capture thedegeneracy inherent in the highly adiabatic versions ofour model, and they reveal that certain choices, for ex-ample: neutrino-neutrino coupling strength Q , neutrinoenergy ratio, and detector location, affect model sensitiv-ity to the unknown parameters to be estimated. Thesefindings are encouraging. In addition, and as discussed inthe previous section, there existed particular realizations4of parameter choices in our simple model that yieldedan MSW resonance prediction that was offset from thetrue location. It will be important to explore methodsof breaking the parameter degeneracies that caused theoffsets.In this section we will: 1) examine the significanceof the results, 2) consider immediate enhancements tothe model and to the D.A. experiment formulation, and3) describe a tentative plan for ultimately investigat-ing D.A. as an alternative method to solving the back-scattering (halo) problem. A. The physics captured by the results
Even in this simplistic two-neutrino model an array offlavor effects is evident. Perhaps the most prominent ofthese is the MSW conversion that occurs for one or bothof the neutrinos in the scenarios depicted in Figs. 1 and2. As the value of Q increases, the nonlinear coupling be-tween ν and ν grows stronger and leads to modificationsboth in the locations of the resonances and in the flavor-space trajectories traversed by the neutrinos as they passthrough resonance. The influence of nonlinearity is espe-cially conspicuous in the scenario with Q = 1000 and anenergy ratio of E ν /E ν = 0 .
01, wherein the resonance of ν , which would have occurred at r ∼ ν , which would have otherwise occurred at r > Q values that we have considered.In the most successful cases the solutions also exhibitmore subtle features of flavor transformation. The smallexcursion in P y that is common to many of the resultsreflects the slight deviation from perfect adiabaticity: thepolarization vector swings away from the xz -plane as itstruggles to keep up with the Hamiltonian vector. Super-imposed on top of these excursions are small-amplitudehigh-frequency oscillations, which correspond to fast pre-cession of P about H , and although there are quanti-tative discrepancies between the model and predictioncurves, the correct qualitative behavior is visible. Takenas a whole, Figs. 1 and 2 suggest that even with just twomeasurements per neutrino the data assimilation proce-dure is nevertheless sensitive to the fine features of flavortransformation. B. Model sensitivity to parameter values
We now comment on our finding that the model’s sen-sitivity to parameter values depends on detector location.As noted, we found that for experiments where the de-tector sat outside the range of resonant flavor transfor-mation of both neutrinos, there existed high degeneracyin permitted values of parameters C and Q . This de-generacy is a result of the flavor conversion of neutrinos having completed prior to detection, thus rendering ameasurement of P z at that location relatively insensitiveto both C and Q independently. In the cases where themeasurement occurred prior (in location) to the completetransformation of ν , however, there emerged a correla-tion between the estimated values of C and Q . This resultoccurred for the energy ratio E ν /E ν = 0 .
01. The corre-lation between C and Q is strongly evident for Q = 1000and more weakly for Q = 100.Our interpretation of these findings is as follows. Forthe Q = 1000 case (Fig. 2: fourth row, right panel), ν isstill evolving through resonance at the location of mea-surement . In this regime, the value of P ,z is highlysensitive to the value of Q at a particular location. Be-cause both V ( r ) and µ ( r ) are taken to scale as 1/ r (seeEqs. 16 and 17), the estimated value of Q essentially setsa corresponding value of C . For the case of Q = 100,the observed correlation between C and Q is weaker be-cause: 1) while some flavor oscillations of ν are apparent, ν is not yet evolving through resonance (Fig. 2: thirdrow, right panel); and 2) Since Q (cid:28) C , the dependenceof resonance location on Q is minimal. This discoveryof detector-location-dependent model sensitivity to pa-rameter values (Sec. IV C) suggests that it will be use-ful to add to a more realistic model additional knownphysics within the supernova envelope—as constraints,rather than as strict measurements (see immediately be-low). C. Next steps
Ideally, we would next like to explore a model of simi-lar simplicity, but with the addition of a back-scatteringterm. Given, however, the degeneracy of parameter esti-mates in this simple model - and hence an unreliable pre-dicted precision of the MSW transition, it is clear that wefirst must improve the current D.A. procedure. Improve-ments may involve amending the model, the D.A. formu-lation, or both. Here we outline a tentative stepwise planfor the ultimate employment of D.A. upon a realistically-sized flavor evolution model that includes terms that areproblematic for traditional numerical integration.1.
Tailor the D.A. procedure to the point at which weare consistently obtaining precise estimates of pa-rameters and MSW transition location for a statis-tically significant fraction of trials - for example,over 99 per cent.
This step may require any of thefollowing amendments, to be discussed in detail be-low: 1) addition to the model of more complicated Note that in this regime the neutrino-neutrino potential is nega-tive , since ν has already gone through resonance and thereforehas P ,z = −
1. This leads to the resonance of ν being pulledfarther in as compared to the case with no ν – ν coupling, an ef-fect that was pointed out in Refs. [44, 45] and is akin to the“matter-neutrino resonance” (MNR) [25, 105–109] add to the model a term that presentsformidable difficulties to traditional numerical inte-gration; namely: back-scattering. No longer able tosolve the model via standard forward integration,we obtain “blind” D.A. solutions. We are now per-mitting the D.A. to lead us in an understanding ofthe new physics that emerges.
3. Ultimately, having identified a D.A. protocol thatworks reliably for a toy model, we scale the model inneutrino number, consider multi-angle calculations,and take the measurements to correspond to energyspectra.
D. Immediate next step: Seek methods fordegeneracy-breaking
1. Amending the model
Degeneracies in outcomes for different neutrino fla-vor evolution histories, for example stemming from adi-abatic evolution through MSW resonances, are amongthe sought-after targets of a numerical calculation. Inmany ways, however, they offer a computational chal-lenge. If we think physically about what could cause non-adiabatic neutrino flavor evolution in a supernova enve-lope we immediately think of physics that gives abruptjumps or ledges in the neutrino forward scattering poten-tials. Shocks or entropy jumps associated with supernovaprogenitor star fossil nuclear burning shells can producethese features in the matter potential. A goal might beto reveal such features by taking a detected Galactic su-pernova core collapse neutrino burst signature, and thenemploying simulations to reverse-engineer the neutrinoflavor histories in the supernova envelope. Such jumpsare tantamount to alterations in the neutrino forward-scattering potentials. Regarding our ongoing study ofD.A. applied to this problem, then, one method to break parameter degeneracies may be to add to the model morecomplicated potentials that capture the physics just de-scribed. Another goal might be to ascertain whether evi-dence of nonlinear neutrino flavor evolution, for example:collective oscillations, appears in such a detected signa-ture.Another method of degeneracy-breaking may be theaddition of constraints within the supernova envelope -that is: between locations r = 0 and R . We have al-ready stipulated that we know the neutrino fluxes, en-ergy spectra, and flavor states at the neutrino sphere.What if, in addition, we stipulate that at a certain epoch(time-post-bounce) in the supernova, and in a range of lo-cation (radius), the polarization vector components andenergy spectra lie in given ranges? That is: we mightwant to examine how well we can fit detected neutrinodata with assumed neutrino sphere flavor and energy dis-tributions, while also demanding, for example, neutron-rich conditions in a certain region in the envelope. Thisneutron excess might arise from the ν e + n → p + e − ,¯ ν e + p → n + e + competition [42], or neutral-and-charged-current neutrino spallation of neutrons from nuclei. Wecan incorporate such considerations in the D.A. proce-dure as constraints on state variable values.Finally, we may add neutrinos, as well as anti-neutrinos, to the model. As neutrinos in any givenmodel are coupled all-to-all, the addition of neutrinoswould tighten the coupling among the model PDEs.In short, prior to adding back-scattering, we considerit important to examine embellishments to currentmodel for degeneracy-breaking. First, we will add morecomplicated potentials (to represent shocks and changesin electron fraction and density at particular locationsr). A more complicated model will not only capturemore faithfully the physics of interest, but it will furnishadditional parameters to be estimated, thereby possiblyproviding more power to the measurements in distin-guishing among solutions. Second, we will add moremeasurements, in terms of neutrino number, to examinethe measurement-dependence of degeneracy-breaking.Third, we will examine the model dynamics at variousepochs post-core bounce. Fourth, and separate from“measurements”, we will add constraints on the values ofvariables within particular ranges of location at variousepochs.
2. Amending the D.A. procedure
In addition to model embellishments, various modifi-cations to the D.A. procedure may improve parameterestimation.First we note our choices for various user-definedsearch parameterizations: the relative weighting terms R m and R f , the coefficient k of the equality constraints,and the choice of discretized step size in radius r . Thereexists no universal optimal choice for any of these values.6The optimal choices are model-dependent and require ex-tensive trial and error to identify. In this paper we ex-perimented with various choices and took the set thatyielded the most reliable results, but our experimenta-tion was not exhaustive. As one example: the relativeweighting terms R f and R m were taken to be equal atall points on any given path. We might examine whetherloosening this constraint changes the results. (For com-ments on choosing an optimal ratio of R m / R f , see Ap-pendix B, Section B 2 b.) Finally, a value of R m = 1.0consistently produced the best results, and this was thevalue used for all results presented in this paper.Second, the user-defined search ranges for parametervalues C and Q differed over all eight experiments, de-pending on the true model values of these quantities. Ob-viously, such bias has no place in a true D.A. experiment- that is: in an experiment where the measurements havebeen generated by a real physical process whose associ-ated model parameter values are, of course, not known.A prerequisite to taking on real (experimental) measure-ments will be the elimination of this bias in simulatedexperiments.More broadly, there exist many formulations for in-corporating measurements and model evolution into anoptimization framework. To discuss the myriad of exist-ing techniques is beyond the scope of this paper. Here wenote just two as examples. The first is variational syn-chronization (or “nudging”), where the model is imposedas a strong constraint upon the cost function rather thanas a term within it. This formulation has demonstratedbetter performance for some problems in neuroscience([110]).As a second example, we cite Monte Carlo (MC) algo-rithms. This is a class entirely separate from variationalapproaches. MC methods have three notable advantagesover the variational method. First, they are more effi-cient at searching a wider area, thereby offering a moreglobal view of the action surface. Second, the map yieldsnot only the depth of a minimum but also its width -thereby offering a quantification of errors on parameterestimates . Third, MC methods are more readily paral-lelizable. This feature may be valuable when consideringscalability. E. Scalability
Ultimately, we will aim to build a D.A. formulation fora large-scale model. Because of the enormity of realisticastrophysical scenarios in terms of model degrees of free-dom, the practicality of scalability is an important con-sideration. In numerical weather prediction, large modelsare being handled by GPUs. The model currently used at On a related note, work is ongoing on a MC method to quan-tify errors on estimates obtained via the variational approachdescribed in this paper [111]. the European Centre for Medium Range Weather Fore-casts, for example, contains 10 degrees of freedom, 10 of which are measured [112]. This size is comparable tothe size of a neutrino flavor evolution model we would ul-timately aim to build. See Ref. [113] for information ontheir computing facilities. Current state-of-the-art nu-merical radiation hydrodynamics simulations of core col-lapse supernovae that incorporate detailed equation ofstate physics and include Boltzmann neutrino transportmight utilize 10 degrees of freedom. Simulations of thissize may be augmented to treat the halo problem, but fullquantum kinetic treatments of neutrino flavor and spinevolution may remain elusive, even at the exascale [114].As noted, MC frameworks for D.A. may be better suitedfor large-scale simulations, as they are readily paralleliz-able. See Ref. [115] for a formulation of the MC algorithmusing a cost function similar to that used in this paper;see Ref. [116] for an exploration of GPU processing ca-pabilities for such MC structures. Finally, we note thatother kinds of statistical approaches to gas dynamics andmagnetohydrodynamics are being pursued, for example,Gaussian Process Modeling [117]. Whether any of thesemight be employed to tackle neutrino flavor evolution isan open question. VII. SUMMARY
Our exploration of optimization-based data assimila-tion techniques for treating the neutrino flavor transfor-mation problem in supernovae has yielded insights intothe utility of these approaches. Advantages include: 1)the ability of D.A. to search efficiently over ranges ofinput model parameters, and 2) the integration-blindfeature offered by an optimization framework. Obvi-ously, we have only scratched the surface of this prob-lem. We plan immediate modifications to both themodel and D.A. procedure, to obtain more precise andreliable results in advance of considering the direction-changing scattering problem in a toy model. Followingthat achievement, we envision scaling up the sophistica-tion of the model, by adding more neutrinos, with manymore neutrino energy bins, and a more sophisticated ge-ometry and realistic non-smooth matter density profilesand direction-changing neutrino scattering. We also envi-sion treating other supernova epochs, for example, shockbreak-out and the associated neutronization burst.
ACKNOWLEDGMENTS
We thank Nirag Kadakia, Paul Rozdeba, Sasha Shir-man, and Daniel Breen for discussions on data assimila-tion and experimental design. We also appreciate Jung-Tsung Li, James Tian, Sebastien Tawa, John Cherry,Brad Keister, Baha Balentekin, and Shashank Shalgarfor discussions on the relevance of D.A. to neutrino astro-physics. Thanks also to Evan Grohs and an anonymous7reviewer, for valuable feedback that has significantly im-proved the paper. This work was supported in part byNSF grants PHY-1307372 and PHY-1614864. Finally,thanks to the village of Doylestown, Ohio—for every-thing.
Appendix A: A path-integral based formulation ofstatistical data assimilation1. Summary of purpose and strategy
Here we lay out a derivation of the cost function A used in this paper. We begin by seeking the probabilityof obtaining a path X in the model’s state space givenobservations Y , or: P ( X | Y ). If we write: P ( X | Y ) = e − A ( X , Y ) , the equation above will then mean: the path X for whichthe probability (given Y ) is greatest is the path that min-imizes A . Now, if A is sufficiently large (where “suf-ficiently” must be defined by the results of a particularD.A. experiment using a particular model), we can useLaplace’s method to estimate the minimizing path on thesurface of A .A formulation for A will permit us to obtain the ex-pectation value of any function G ( X ) on a path X ; ex-pectation values are the quantities of interest when theproblem is statistical in nature. We can write the expec-tation value of G ( X ) as: (cid:104) G ( X ) (cid:105) = (cid:82) d X G ( X ) e − A ( X , Y ) (cid:82) d X e − A ( X , Y ) . That is: the expectation value can be expressed as aweighted sum over all possible paths, where the weightsare exponentially sensitive to A . The RMS variation,and higher moments of G ( X ), can be calculated by tak-ing the x a to the appropriate higher exponents. If thequantity of interest is the path X itself, then we choose G ( X ) = X .It remains, then, to write a functional form for A .This will take place in two steps. First we shall considerhow measurements and model dynamics enter into theprocess state and parameter estimation. This we will dovia an examination of Bayesian probability theory andMarkov chain transition probabilities, for the effect ofmeasurements and model dynamics, respectively. Sec-ond, we shall make four simplifying assumptions: 1) themeasurements taken at different times are independent;2) both measurement and model errors have Gaussiandistributions; 3) each measurement is taken to corre-spond directly to one model state variable; 4) the min-imizing path is independent of the guess - in state andparameter space - of the initial path.In what follows, we shall describe this strategy in somedetail (for a full treatment, see Ref.[1]). To remind the Laplaces method was developed to approximate integrals of theform: (cid:82) e Mf ( x ) dx . For sufficiently high values of the coefficient M , significant contributions to the integral will come only frompoints in a neighborhood around the minimum, which can thenbe estimated. D PDEs,each of which represents the evolution of one of themodel’s D state variables. From the corresponding phys-ical system, we are able to measure L quantities, eachof which corresponds to one of the model’s D state vari-ables. Typically the measurements are sparse ( L (cid:28) D ),and the sampling may be infrequent or irregular.
2. Considering model dynamics only (nomeasurements yet)
We shall first examine this formulation by consider-ing the model’s time evolution in the absence of mea-surements. We represent the model’s path through statespace as the set X = { x ( t ) , x ( t ) , . . . , x ( t N ) , p } , where t N is the final “time point” and the vector x ( t ) containsthe values of the D total state variables, and p are theunknown parameters (here, the phrasing “time” can alsobe taken to represent other grid parameterizations; forinstance: location). a. Assuming that a Markov process underlies the dynamics If we assume that the dynamics are memory-less, orMarkov, then x ( t ) is completely determined by x ( t − ∆ t ),where t − ∆ t means: “the time immediately preceding t ” and an appropriate discretization of time ∆ t for ourparticular model has been chosen. A Markov processcan be described in the continuous case by a differentialequation, or as a set of differential equations:d x a ( t )d t = F a ( x ( t ) , p ); a = 1 , , . . . , D, and we note that the model is an explicit function of thestate variables x ( t ) and the unknown parameters p . It isin this way that the unknown parameters are consideredto be on equal footing with the variables; namely: theyare variables with trivial dynamics.In discrete time, that relation can be written in variousforms. For our purposes, we use the trapezoidal rule: x a ( n + 1) = x a ( n ) + ∆ t F a ( x ( n + 1)) + F a ( x ( n ))] , where for simplicity we have taken n and n + 1 to repre-sent the values of t n and t n +1 . b. Permitting stochasticity in the model and recasting itsevolution in terms of probabilities We are interested in ascertaining the model evolutionfrom time step to time step, where now we allow for somestochasticity in the model dynamics. In this scenario,the evolution can be formulated in terms of “transitionprobabilties”, e.g., P ( x ( n + 1) | x ( n ))—the probability of the system reaching a particular state at time n + 1 givenits state at time n . If the process were deterministic, thenin our case P ( x ( n + 1) | x ( n )) would simply reduce to: δ D ( x ( n + 1) − x ( n ) − ∆ t [ F ( x ( n + 1)) + F ( x ( n ))]). Wewill revisit to this expression later in this section, under Approximating the Action .For a Markov process, the transition probability fromstate x ( n ) to state x ( n + 1) represents the probabilityof reaching state x ( n + 1) given x ( n ) and x at all priortimesteps. Or: P ( x ( n + 1) | x ( n )) = P ( x ( n + 1) | x ( n ) , x ( n − , . . . , x (0))so that P ( X ) ≡ P ( x (0) , x (1) , . . . , x ( N ))= N − (cid:89) n =0 P ( x ( n + 1) | x ( n )) P ( x (0)) . We now write P ( X ) ≡ e − A ( X ) , where A is the action defined on the model’s path X instate space (or: the path that minimizes the Action is thepath most likely to occur ) . Then the model term of theAction, A , model , can be written: A , model = − (cid:88) log[ P ( x ( n + 1) | x ( n ))] − log[ P ( x (0))] , where the second term represents uncertainty in initialconditions.
3. Now with measurements
We now consider the effect of measurements. Let usdefine a complete set of measurements Y to be the setof all vectors y ( n ) at all times n —the analog of X forthe complete set of state variable values. We shall ex-amine the effect of these measurements upon a model’sdynamics by invoking the framework of “conditional mu-tual information” (CMI); for a useful definition of CMI,see Ref. [118] . The reader might find it of interest to note the quantum-mechanical analog of the transition probability, which involvesthe trivial addition of the term i (cid:126) in the exponent: P ( x ( n +1) | x ( n )) = e i (cid:126) A ( t n +1 ,t n ) , where A here is the classical action. The reader may find an intuitive understanding of our use ofthe CMI by the following consideration. The overall informa-tion, in bits, in a set A is defined as the Shannon entropy H ( A ) = − (cid:80) A P ( A ) log[ P ( A )]. The CMI is a means to quan-tify the amount of information, in bits, that is transferred alonga model trajectory within a particular temporal window. Thatinformation is equivalent to: − N (cid:88) n =0 log[ P ( x ( n ) | y ( n ) , Y ( n − x ( n ) , y ( n ) | Y ( n − x ( n ) upon observing event y ( n ), conditioned on having previously observed event(s) Y ( n − x ( n ) , y ( n ) | Y ( n − (cid:20) P ( x ( n ) , y ( n ) | Y ( n − P ( x ( n ) | Y ( n − P ( y ( n ) | Y ( n − (cid:21) .
4. The complete Action
With measurement considerations included, the actionnow becomes: A ( X , Y ) = − (cid:88) log[ P ( x ( n + 1) | x ( n ))] − log[ P ( x (0))] − (cid:88) CMI( x ( n ) , y ( n ) | Y ( n − , where the first and second terms represent the model dy-namics including initial conditions, and the third termrepresents the transfer of information from measure-ments. The summations are over time. As noted, this for-mulation positions us to calculate the expectation valueof any function G ( X ) on the path X .We now offer an interpretation of the measurementterm. The measurement term can be considered to be anudging (or synchronization) term. While nudging termsare often introduced rather artificially in the interest ofmodel control, however, we have shown that the mea-surement term arises naturally through considering theeffects of the information those measurements contain.For this reason, we prefer to regard the measurementterm as a guiding potential. In the absence of the poten-tial, we live in a state space restricted only by our model’sdegrees of freedom. The introduction of the measure-ments guides us to a solution within a sub space in whichthose particular measurements are possible.
5. Approximating the Action
We now seek to simplify the Action formulation for thepurposes of calculation. a. The measurement term
Regarding the measurement term, we make four as-sumptions: • The measurements taken at different times are in-dependent of each other. This permits us to writethe CMI simply as: P ( x ( n ) | y ( n )). Or: A ( X , Y ) = − log[ P ( X | Y )] . • There may be an additional relation between themeasurements and the state variables to whichthose measurements correspond, which can be ex-pressed with the use of some transfer function h l : h l ( x ( n )) = y l ( n ). • For each of the L measured state variables, weallow for a noise term θ l at each timepoint, foreach measurement y l that corresponds to a statevariable x l : y l ( n ) = h l ( x ( n )) + θ l ( n ). In thiscase, then, P ( x ( n ) | y ( n )) is simply some functionof h ( x ( n )) − y ( n ) at each timepoint. • The measurement noise has a Gaussian distribu-tion.Taking these assumptions, we arrive at:CMI( x ( n ) , y ( n ) | Y ( n − − L (cid:88) l,k =1 ( h l ( x ( n )) − y l ( n )) [ R m ( n )] lk h k ( x ( n )) − y k ( n )) , where R m is the inverse covariance matrix of the mea-surements y l . b. The model term We simplify the model term by assuming that themodel may have errors, which will broaden the delta func-tion in the expression noted earlier for the deterministiccase. If we assume that the distribution of errors is Gaus-sian, then δ D ( z ) becomes: (cid:114) det R f (2 π ) D e [ − z R f z ], where R f is the inverse covariance matrix for the model’s state vari-ables.Taking both approximations together, assuming thatthe transfer function h l is simply unity, and assumingthat the minimizing path is independent of considera-tions of initial conditions, we obtain: A = N − (cid:88) n D (cid:88) a R fa (cid:18) x a ( n + 1) − x a ( n ) t n +1 − t n − f a ( x ( n )) (cid:19) + (cid:88) j L (cid:88) l R ml y l ( j ) − x l ( j )) , where f a ( x ( n )) ≡ [ F a ( x ( n )) + F a ( x ( n + 1))]. The first(model) term involves a summation over all D state vari-ables, and the second (measurement) term involves asummation over the L measured quantities. Note thathere we write the model error term in a simpler, moregeneral manner than the specific formulation used in thispaper (Eq. 1).Finally, we allow in the cost function the addition ofequality constraints, of the general form k g ( x ( n )), where0the coefficient k set the strength of the constraint func-tion g . The specific equality constraints chosen in thispaper are described in Sec. II. Appendix B: Details of the D.A. procedure1. The discretized search space
The optimization procedure searches a ( D ( N +1)+ p )-dimensional state space, where D is the number of statevariables of a model, N is the number of discretized steps,and p is the number of unknown parameters. Note thateach location point is considered a separate dimension.Thus the action, instead of being a functional of D func-tions, is a function of ( D ( N + 1) + p ) variables.Ipopt, the specific algorithm used in this paper, em-ploys a Newton’s, or descent-only, search. The spatialresolution is set by a user-defined step size. The userprovides the objective function, model, the Jacobean andHessian matrices of the model, permitted search rangesof variables and unknown parameters, and discretizedstep size. The algorithm iteratively searches for a pathin the state space that minimizes the action subject tothe requirements that the first derivative of the objectivefunction at the minimizing path along any direction bezero and that its second derivative along any direction bepositive definite. The resulting “path” is a set of statevectors, one at each discretized step, and specific valuesof the unknown parameters. Each path corresponds to asingle point in the ( D ( N + 1) + p )-dimensional space. Inthis way, the model parameters are considered on equalfooting with the state variables; namely: the unknownparameters are state variables with trivial dynamics. Fi-nally, to impose user-defined bounds placed upon thesearches, Ipopt uses a barrier method. For details, see[90].
2. Specific choices governing D.A. experiments inthis paper a. Interface with Ipopt
Ipopt requires a user interface to discretize statespace and calculate the model equations of motion,Jacobean, and Hessian matrices that are used in theminimization procedure. We used a suite of Pythoncodes to generate this interface; it is available here:https://github.com/yejingxin/minAone. b. Choosing R f / R m for best results As noted in Section VI, there exists no universal rulefor choosing an optimal ratio of model and measurementweights. An optimal value is model-dependent and must be identified via trial-and-error. Generally, for many bio-physical models of neurons, small neuronal networks, at-mospheres, and chaotic Lorenz-63 and Lorenz-96 models,a value of β between 10-20 is found to be ideal (privatecommunications 2017). The reader may compare thisrange to our identification of β ∈ [13 , R m (cid:29) R f and R f (cid:29) R m ) are expected for any model, for the following rea-sons. For low R f , the model constraints are not yet suf-ficiently strict to require a converging solution. For high R f , the failure of solutions has at least two potentialcauses. First, one encounters numerical problems withconsidering “infinite” model weight. The problem is ill-conditioned when it involves a matrix whose elements areso large that the matrix is not invertible. The optimizingsolution may thus become overly sensitive to changes inthe state vector. Rounding error may render these so-lutions invalid. A second possible cause is discretizationerror at high R f . In taking a discretized derivative, oneretains only the first term in a Taylor series. As the mul-tiplicative factor grows, the higher-order terms - whichare ignored - will become important. Appendix C: Embedding the model into a simplifiedastrophysical system1. Forms for matter and coupling potentials
The cubic radial dependence of the matter potential isactually close to the expected density run in the super-nova envelope in some cases. For example, some secondsafter a supernova explosion, perhaps 3 to 10 s post-corebounce, we can be left with a tenuous, near-hydrostaticenvelope sitting in a gravitational potential well domi-nated by the hot, proto-neutron star. This envelope isbeing heated to high entropy by the intense neutrino ra-diation from the neutrino sphere, and driven off. This isthe “neutrino-driven wind” epoch. It is a candidate sitefor r -process nucleosynthesis, but one fraught with chal-lenges stemming from uncertain neutrino flavor transfor-mation physics and the “alpha effect”: the interactionbetween charged current ν e and ¯ ν e captures and aggres-sive alpha particle formation in the high entropy wind. Inturn, the entropy of this wind is a complicated functionof neutrino heating and flavor histories.We can approximate the wind regime envelope as:(1) being in hydrostatic equilibrium, with enthalpy perbaryon equal to the local gravitational binding energyper baryon; and (2) with the entropy of the material be-ing carried entirely by relativistic particles, namely pho-tons and electron-positron pairs. The latter assumptionis tantamount to the entropy being high. We can com-bine (1) and (2) and find for a constant entropy envelopethe baryon density dependence on radius r :1 n B ( r ) = ρ ( r ) N A = (cid:18) π (cid:19) g (cid:34) M NS m p m (cid:35) s r , (C1)where the baryon mass (energy) density is ρ , Avogadro’snumber is N A , and M NS , m p , and m pl are, respectively,the neutron star mass, proton rest mass, and the Planckmass. Here s is the entropy-per-baryon in units of Boltz-mann’s constant k b , and g is the statistical weight inrelativistic particles. In terms of this parametrization ofthe density run, our constant C in the expression for V ( r )is [119] C = √ G F Y e (cid:18) π (cid:19) (cid:18) g / (cid:19) (cid:34) M NS m p m (cid:35) s ≈ . × MeV cm (cid:18) g / (cid:19) (cid:20) M NS . (cid:12) (cid:21) Y e s , (C2)where Y e is the electron fraction, and the entropy-per-baryon in units of 100 k b is s .In a spherical geometry, with neutrinos emitted from a sharp neutrino sphere, the radial dependence of the ν – ν potential is µ ( r ) ∼ /r , as both the neutrino numberflux n ν ( r ) and the angle factor α ( r ) each dilute as 1 /r in the far-field limit. In more complicated models, includ-ing those that incorporate back-scattering, we expect adifferent radial dependence than 1 /r . Specifically, weexpect the neutrino potential to drop less quickly withradius than in conventional bulb models. Here then,for simplicity and to enable a direct comparison to thematter potential, we choose µ ( r ) ∼ /r . That is: ourmotivation here was to introduce nonlinearity in a sim-ple manner, while avoiding the use of different functionalforms for V ( r ) and µ ( r ).
2. Neutrino energy ratios set within anastrophysical context
For the energy ratio E ν /E ν = 2 . s = 1, g = 11 /
2, and Y e = 0 . > E ν = 25 MeV and E ν = 10 MeV are289 km and 213 km, respectively, and the ratio in Eq. 19is ≈ .
4. Note that the corresponding resonance widths, δr = | V / ( dV /dr ) | res tan 2 θ ∼ ( r res /
3) sin 2 θ ∼
10 km for θ = 0 .
1, are small enough that the resonances are wellseparated for these neutrino energies. We can also con-sider the same neutrino energies, but now with a smallerentropy, s = 0 .
1, a slightly smaller electron fraction, Y e = 0 .
35, and g -factor, g = 2. These choices willvery crudely mock up an earlier accretion phase super-nova envelope. In this case the resonance locations are 4254 km and 3135 km, respectively. If we consider thesame envelope parameters but now take neutrino ener-gies E ν = 2 . E ν = 1 MeV, we obtain reso-nance locations at 1975 km and 1455 km, respectively. Inall of these cases neutrino flavor evolution through theseresonances will be adiabatic.If we take the neutrino energy ratio of 0 .
01, with E ν = 0 . E ν = 50 MeV, and the wind-likehigher entropy conditions described above, we obtain res-onance locations at 79 km and 364 km, respectively, forthe Q = 0 case. In this case, our experimental setupwould put the final location R between these resonances,inside the supernova envelope. We study this scenario,for multiple values of coupling strength Q , in order toexamine collective effects and explore the sensitivity ofthe D.A. procedure to flavor information deep in the su-pernova envelope. Appendix D: Evolution of massive stars, weakinteractions, and neutrino flavor physics
The following is a pedagogical overview of neutrinophysics in core collapse supernovae [120–122].
1. Evolution of massive stars and the weakinteraction
The weak interaction, the nuclear force responsible forchanging neutrons to protons and vice versa, is the keyto why stars shine, and why big stars collapse, explode,and synthesize the elements. The sun and stars like itburn hydrogen into helium, combining four protons intoa helium nucleus, and thereby turning two of those pro-tons into neutrons along the way. The fundamental weakreaction in the sun turns two protons into a deuteriumnucleus with the emission of a positron and an accompa-nying electron-flavor neutrino, p + p → D + e + + ν e .Neutrinos experience only gravitation and the weakforce, making them very “slippery,”; that is: able to es-cape from deep inside a dense star, and carry away en-ergy. The weak interaction is aptly named, being sometwenty orders of magnitude weaker than electromagneticforces, at the relevant energy scales. Indeed, hydrogenburning in the sun is desperately slow. It will take 10 years for the sun to burn through all of its hydrogen.In more massive stars, however, weak interactions, alongwith attendant neutrino emission, combined with grav-itation, can nevertheless engineer their violent destruc-tion.Stars some ten or more times the mass of the sun( M ≥
10 M (cid:12) ) evolve in millions of years through a se-ries of nuclear burning epochs: hydrogen to helium, tocarbon and oxygen, to magnesium, to silicon. Finally, sil-icon burns to “iron,” forming a core with mass ∼ . (cid:12) composed of relatively neutron-rich iron-peak nuclei (forexample, Fe, Ca, etc.). From core carbon burning2onward in these objects, the energy carried away byneutrinos exceeds that radiated by photons! Neutrinoscarry away the heat generated by nuclear reactions, forc-ing the star to contract and release more gravitationalbinding energy, accelerating nuclear burning, and so on.This jams the electrons in the star into a smaller andsmaller volume, and the Pauli principle implies that theyare consequently forced into higher and higher energystates—the electrons become relativistically degenerate.In turn, the energy dependent weak interactions, for ex-ample, electron capture on protons to make neutrons( e − + p → n + ν e ) proceed faster. Though this iron corehas a density more than ten orders of magnitude higherthan that of water, it is essentially transparent to theseneutrinos.The end result is that neutrinos leave and refrigeratethe core. Though the core has a temperature of nearly10 K ( ∼ ∼ ∼ g cm − , roughly one percent ofnuclear matter density, it becomes opaque to neutrinos.The neutrinos are trapped and quickly come into thermaland chemical equilibrium with the matter. As the col-lapse proceeds, the outer portions of the core are fallingin supersonically. When the inner part of the core reachesnuclear density, the nucleons touch, and this region stopsabruptly. The outer, supersonic part of the core slamsinto this “brick wall,” generating a shock wave that prop-agates outward through the outer core.In broad brush terms, the ∼ . (cid:12) core collapsesfrom a configuration with a radius like that of the earth( ∼ cm) to one with a radius of roughly 45 km inabout one second. Within another second or two it quasi-statically shrinks down to a radius of 10 km. The up-shot is a prodigious gravitational binding energy change,amounting to about ten percent of the entire rest mass ofthe core. One percent of this energy largess resides in thebulk in-fall kinetic energy of the core (and consequentlythe initial energy in the shock wave), and the other 99percent, some 10 erg, is in the trapped seas of neutrinosof all kinds.At the edge of the proto-neutron star, deemed the“neutrino sphere,” the matter density, and opacity to neutrinos, drops off dramatically and neutrinos canmore or less freely stream away, mostly unhindered bydirection-changing or inelastic collisions with particlesthat carry weak charge, for example, neutrons and pro-tons, electrons, and other neutrinos. The average en-ergies of the neutrinos streaming out are of order ∼
10 MeV. With a gravitational binding energy of 10 erg( ∼ MeV), this amounts to some 10 neutrinos car-rying this energy away in a matter of a few seconds.These are titanic neutrino fluxes.The shock wave that propagates through the super-nova envelope is associated with an entropy jump acrossthe shock front—the material that the shock plows intohas an entropy-per-baryon ∼ k B , whereas the materialbehind the shock has an entropy-per-baryon of ∼ k B .As a result, the passing of the shock wave through theenvelope results in the dissociation of nuclei into mostlyfree nucleons, a process that costs ∼ ν e + n → p + e − and ¯ ν e + p → n + e + ) maydeposit enough energy in the matter behind the shock tore-energize it and get it moving again with an energy of10 erg, resulting in a supernova explosion. This processcan be aided by hydrodynamic motion of the neutrino-heated material. In the end, about one percent of thetotal neutrino energy needs to be deposited in this ma-terial to get an explosion.
2. Collective neutrino flavor transformations insupernovae
It is known that neutrinos come in three “flavors”, ν e , ν µ , and ν τ , corresponding to each of the three chargedleptons. These flavors denote weak-interaction eigen-states, essentially determining how these particles inter-act in matter. Each neutrino has an antiparticle, imply-ing that there are six kinds of neutrinos: ν e , ¯ ν e , ν µ , ¯ ν µ , ν τ , ¯ ν τ . These particles are spin-1 /
2, electrically neutral,and have very small rest masses. We do not know whatthe masses are, but the differences of the squares of thesemasses are measured: the so-called solar mass-squaredsplitting δm (cid:12) = m − m ≈ . × − eV , and the at-mospheric mass-squared splitting δm = m − m ≈ . × − eV , where m , m , and m are the neutrinomass eigenvalues corresponding to the energy eigenstates(sometimes called “mass” states) of the neutrinos. Ex-periment shows that these neutrino mass states are notcoincident with the flavor states and this can have conse-quences for the core-collapse supernova mechanism andfor neutrino detection.The fact that neutrino mass states are not coincidentwith flavor states means that neutrinos emitted initiallyin one flavor state can transform into another as theypropagate, with consequences for the way these parti-3cles effect heating, nucleosynthesis, etc. Flavor transfor-mations are modified in the presence of potentials aris-ing from neutrino forward scattering on particles thatcarry weak charge, such as leptons, nucleons, and otherneutrinos. As the neutrinos stream away through thelower density material above the neutrino sphere, theyacquire through forward scattering an “index of refrac-tion”, equivalent to an effective mass in medium. This isanalogous to the way photons acquire an index of refrac-tion and effective mass propagating through a transpar-ent medium like glass. Unlike this optical case, however,the “medium” through which the supernova neutrinospass consists, in part, of other neutrinos! This makes theneutrino flavor transformation problem fiercely nonlin-ear: the potentials that determine how neutrinos changetheir flavors depend on the flavor states of the neutrinos.These nonlinear effects become important in environ-ments where the neutrino fluxes are substantial, such ascore-collapse supernovae, compact object mergers, andalso the early universe. A complete treatment of flavor-transformation physics in these environments is impor-tant, because the charged-current weak interactions areflavor-dependent at typical temperatures ( ∼ MeV)—the ν e ’s participate, but ν µ ’s and ν τ ’s do not as there areno µ or τ leptons around to scatter on. As a result, theeffective scattering cross-sections for ν e are larger thanthose for ν µ and ν τ , resulting in different energy deposi-tion rates—relevant for the supernova explosion mech-anism. Moreover, the charged-current weak processes ν e + n → p + e − and ¯ ν e + p → n + e + determine the n/p ratio, and therefore knowing the flavor content isessential for evaluating the nucleosynthesis prospects inthese environments.Thermal processes during the core collapse manufac-ture neutrino-antineutrino pairs of all flavors and thesethermalize with the electron capture-created ν e ’s. Thenet result is a rough equipartition of energy among allsix types of neutrinos. Neutrinos of different flavors,however, have correspondingly different interactions inthe matter near the neutrino sphere. The result is thatelectron-flavor neutrinos, with the largest interactions,decouple furthest out, where it is coolest, and have loweraverage energies as a consequence. µ and τ flavor neutri-nos and their antiparticles have no charged-current weakinteractions, and so these neutrinos decouple deeper in,where it is hotter. Consequently, these are on averagemore energetic. Electron antineutrinos have energies inbetween those of the electron neutrinos and the µ or τ flavor neutrinos.Neutrinos diffuse out of the hot proto-neutron star corewith a typical random walk time of seconds. This ratherlong diffusion time also sets the timescale over which neutrino spectral parameters and fluxes change. Thetimescale for these changes can then be long compared toneutrino transit times across regions of interest. Numer-ical studies of supernova neutrino flavor evolution havetraditionally sought to take advantage of this situationby seeking stationary, time-independent solutions to theevolution equations, wherein the neutrino fluxes/spectradepend only on position. These numerical studies, inwhich some ∼ nonlinearly-coupled Schr¨odinger-likeequations are solved on a supercomputer, have yieldedunexpected and surprising results [2–47]. Nonlinearity inthe neutrino flavor potentials can give rise to collectiveneutrino flavor oscillations, where significant populationsof neutrinos in the supernova envelope can execute syn-chronized or other organized and simultaneous changesin flavor, across a range of neutrino energies and in alarge region of space or time.One of the limitations of current simulations of neu-trino flavor evolution in supernovae is the failure toaccount for potentials arising from neutrino direction-changing scattering. This is the “neutrino halo” effect.Even though a relatively small fraction of neutrinos un-dergo direction-changing scattering, they could neverthe-less contribute significantly to the forward-scattering po-tential felt by the outward-streaming neutrinos. This isa consequence of the peculiar intersection-angle depen-dence of the weak-interaction potential. In certain re-gions of the envelope, and for certain epochs, it has beenshown that the potential term arising from the halo neu-trinos could in fact be the dominant term [123, 124]. Acomplete treatment of neutrino flavor evolution that in-cludes the effects of both forward and direction-changingscattering, necessitates the use of the so-called “Quan-tum Kinetic Equations” (QKE) [16, 54–70, 93, 94]. Inhigh-density regions, where the scattering rates are largeso that quantum mechanical phases do not have anytime to build up, the QKEs reduce to a Boltzmann-likeform. In the other limit, where the neutrinos essentiallyfree-stream and only experience coherent forward scat-tering, the QKEs reduce to a Liouville-von Neumann(Schr¨odinger-like) equation.If in the future we are lucky enough to detect the neu-trino burst from a Galactic core collapse event, we willwant to know whether the detected signal indicates thatthe simple forward-scattering-based optical analogy issufficient to explain the neutrino flavor data, or whetherthe halo must be invoked. A key objective will be to usethis signal to potentially extract information regardingthe conditions in the envelope and to ascertain whethercollective oscillations, and their signatures like spectralswaps/splits occurred. These issues prompt the explo-ration of alternative calculation techniques. [1] H. Abarbanel, Predicting the Future: Completing Mod-els of Observed Complex Systems (Springer Science &Business Media, 2013). [2] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z.Qian, Phys. Rev. D , 105014 (2006), arXiv:astro-ph/0606616. [3] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian,Physical Review Letters , 241101 (2006), arXiv:astro-ph/0608050.[4] H. Duan, G. M. Fuller, and Y.-Z. Qian, Phys. Rev. D , 123004 (2006), arXiv:astro-ph/0511275.[5] H. Duan, G. M. Fuller, and Y.-Z. Qian, Phys. Rev. D , 085013 (2007), arXiv:0706.4293.[6] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z.Qian, Phys. Rev. D , 125005 (2007), arXiv:astro-ph/0703776.[7] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z.Qian, Physical Review Letters , 241802 (2007),arXiv:0707.0290.[8] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z.Qian, Physical Review Letters , 021101 (2008),arXiv:0710.1271.[9] H. Duan, G. M. Fuller, and Y.-Z. Qian, Phys. Rev. D , 085016 (2008), arXiv:0801.1363.[10] H. Duan, G. M. Fuller, and J. Carlson, Comput. Sci.Dis. , 015007 (2008), arXiv:0803.3650 [astro-ph].[11] J. F. Cherry, G. M. Fuller, J. Carlson, H. Duan,and Y.-Z. Qian, Phys. Rev. D , 085025 (2010),arXiv:1006.2175 [astro-ph.HE].[12] H. Duan, G. M. Fuller, and Y.-Z. Qian, Annual Re-view of Nuclear and Particle Science , 569 (2010),arXiv:1001.2799 [hep-ph].[13] H. Duan and A. Friedland, Physical Review Letters ,091101 (2011), arXiv:1006.2359 [hep-ph].[14] J. F. Cherry, M.-R. Wu, J. Carlson, H. Duan, G. M.Fuller, and Y.-Z. Qian, ArXiv e-prints (2011),arXiv:1108.4064 [astro-ph.HE].[15] J. F. Cherry, M.-R. Wu, J. Carlson, H. Duan, G. M.Fuller, and Y.-Z. Qian, Phys. Rev. D , 125010 (2012),arXiv:1109.5195 [astro-ph.HE].[16] Y. Zhang and A. Burrows, Phys. Rev. D , 105009(2013), arXiv:1310.2164 [hep-ph].[17] A. B. Balantekin, in American Institute of Physics Con-ference Series , American Institute of Physics Confer-ence Series, Vol. 1594, edited by S. Jeong, N. Imai,H. Miyatake, and T. Kajino (2014) pp. 313–318,arXiv:1401.5818 [astro-ph.SR].[18] C. Lunardini, in
American Institute of Physics Confer-ence Series , American Institute of Physics ConferenceSeries, Vol. 1666 (2015) p. 070001.[19] A. B. Balantekin, in
American Institute of PhysicsConference Series , American Institute of PhysicsConference Series, Vol. 1743 (2016) p. 040001,arXiv:1605.00578 [nucl-th].[20] E. Akhmedov and A. Mirizzi, Nuclear Physics B ,382 (2016), arXiv:1601.07842 [hep-ph].[21] S. Chakraborty, R. S. Hansen, I. Izaguirre, and G. G.Raffelt, Journal of Cosmology and Astroparticle Physics , 028 (2016), arXiv:1507.07569 [hep-ph].[22] R. Barbieri and A. Dolgov, Nuclear Physics B , 743(1991).[23] C. Volpe, ArXiv e-prints (2016), arXiv:1609.06747[astro-ph.HE].[24] A. B. Balantekin and Y. Pehlivan, Journal of Physics GNuclear Physics , 47 (2007), arXiv:astro-ph/0607527.[25] L. Johns, M. Mina, V. Cirigliano, M. W. Paris,and G. M. Fuller, Phys. Rev. D , 083505 (2016),arXiv:1608.01336 [hep-ph].[26] G. Raffelt, S. Sarikas, and D. de Sousa Seixas, ArXive-prints (2013), arXiv:1305.7140 [hep-ph]. [27] S. Sarikas, I. Tamborra, G. Raffelt, L. H¨udepohl,and H.-T. Janka, Phys. Rev. D , 113007 (2012),arXiv:1204.0971 [hep-ph].[28] G. G. Raffelt and A. Y. Smirnov, Phys. Rev. D ,081301 (2007), arXiv:0705.1830.[29] S. Hannestad, G. G. Raffelt, G. Sigl, and Y. Y. Y.Wong, Phys. Rev. D , 105010 (2006), arXiv:astro-ph/0608695.[30] B. Dasgupta, A. Dighe, G. G. Raffelt, and A. Y.Smirnov, Physical Review Letters , 051105 (2009),arXiv:0904.3542.[31] D. N¨otzold and G. Raffelt, Nuclear Physics B , 924(1988).[32] S. Pastor, G. Raffelt, and D. V. Semikoz, Phys. Rev. D , 053011 (2002), arXiv:hep-ph/0109035.[33] B. Dasgupta, A. Dighe, A. Mirizzi, and G. Raffelt,Phys. Rev. D , 033014 (2008), arXiv:0805.3300.[34] J. F. Cherry, J. Carlson, A. Friedland, G. M. Fuller,and A. Vlasenko, Phys. Rev. D , 085037 (2013),arXiv:1302.1159 [astro-ph.HE].[35] J. F. Cherry, A. Vlasenko, G. M. Fuller, and J. Carlson,preprint (2012).[36] B. Dasgupta and A. Mirizzi, Phys. Rev. D , 125030(2015), arXiv:1509.03171 [hep-ph].[37] S. Abbar and H. Duan, Physics Letters B , 43(2015), arXiv:1509.01538 [astro-ph.HE].[38] H. Duan and S. Shalgar, Physics Letters B , 139(2015), arXiv:1412.7097 [hep-ph].[39] H. Duan, International Journal of Modern Physics E ,1541008 (2015), arXiv:1506.08629 [hep-ph].[40] S. Abbar, H. Duan, and S. Shalgar, Phys. Rev. D ,065019 (2015), arXiv:1507.08992 [hep-ph].[41] H. Duan and S. Shalgar, Journal of Cosmology andAstroparticle Physics , 084 (2014), arXiv:1407.7861[hep-ph].[42] Y.-Z. Qian, G. M. Fuller, G. J. Mathews, R. W. Mayle,J. R. Wilson, and S. E. Woosley, Physical Review Let-ters , 1965 (1993).[43] Y.-Z. Qian and G. M. Fuller, in American Astronomi-cal Society Meeting Abstracts , Bulletin of the AmericanAstronomical Society, Vol. 24 (1992) p. 1264.[44] Y.-Z. Qian and G. M. Fuller, Phys. Rev. D , 1479(1995), arXiv:astro-ph/9406073.[45] Y.-Z. Qian and G. M. Fuller, Phys. Rev. D , 656(1995), astro-ph/9502080.[46] A. Banerjee, A. Dighe, and G. Raffelt, Phys. Rev. D , 053013 (2011), arXiv:1107.2308 [hep-ph].[47] A. Friedland, Physical Review Letters , 191102(2010), arXiv:1001.0996 [hep-ph].[48] R. F. Sawyer, Phys. Rev. D72 , 045003 (2005),arXiv:hep-ph/0503013 [hep-ph].[49] R. F. Sawyer, Phys. Rev. Lett. , 081101 (2016),arXiv:1509.03323 [astro-ph.HE].[50] I. Izaguirre, G. Raffelt, and I. Tamborra, Phys. Rev.Lett. , 021101 (2017), arXiv:1610.01612 [hep-ph].[51] B. Dasgupta, A. Mirizzi, and M. Sen, JCAP , 019(2017), arXiv:1609.00528 [hep-ph].[52] R. F. Sawyer, Phys. Rev.
D79 , 105003 (2009),arXiv:0803.4319 [astro-ph].[53] J. F. Cherry, J. Carlson, A. Friedland, G. M. Fuller,and A. Vlasenko, Physical Review Letters , 261104(2012), arXiv:1203.1607 [hep-ph].[54] A. Vlasenko, G. M. Fuller, and V. Cirigliano, ArXive-prints (2014), arXiv:1406.6724 [astro-ph.HE]. [55] C. Volpe, D. V¨a¨an¨anen, and C. Espinoza, Phys. Rev.D , 113010 (2013), arXiv:1302.2374 [hep-ph].[56] A. Vlasenko, G. M. Fuller, and V. Cirigliano, Phys.Rev. D , 105004 (2014), arXiv:1309.2628 [hep-ph].[57] J. Serreau and C. Volpe, Phys. Rev. D , 125040(2014), arXiv:1409.3591 [hep-ph].[58] P. A. Andreev, Physica A Statistical Mechanics and itsApplications , 108 (2015).[59] C. Volpe, International Journal of Modern Physics E , 1541009 (2015), arXiv:1506.06222 [astro-ph.SR].[60] V. Cirigliano, G. M. Fuller, and A. Vlasenko, PhysicsLetters B , 27 (2015), arXiv:1406.5558 [hep-ph].[61] A. Kartavtsev, G. Raffelt, and H. Vogel, Phys. Rev. D , 125020 (2015), arXiv:1504.03230 [hep-ph].[62] A. Dobrynina, A. Kartavtsev, and G. Raffelt, Phys.Rev. D , 125030 (2016), arXiv:1605.04512 [hep-ph].[63] A. Chatelain and C. Volpe, ArXiv e-prints (2016),arXiv:1611.01862 [hep-ph].[64] C. Volpe, Journal of Physics Conference Series ,062068 (2016), arXiv:1601.05018 [astro-ph.HE].[65] D. N. Blaschke and V. Cirigliano, Phys. Rev. D ,033009 (2016), arXiv:1605.09383 [hep-ph].[66] G. Sigl and G. Raffelt, Nuclear Physics B , 423(1993).[67] G. Raffelt, G. Sigl, and L. Stodolsky, Physical ReviewLetters , 2363 (1993), arXiv:hep-ph/9209276.[68] G. Raffelt and G. Sigl, Astroparticle Physics , 165(1993), arXiv:astro-ph/9209005.[69] P. Strack and A. Burrows, Phys. Rev. D , 093004(2005), arXiv:hep-ph/0504035.[70] C. Y. Cardall, Phys. Rev. D , 085017 (2008),arXiv:0712.1188.[71] K. Abe, T. Abe, H. Aihara, Y. Fukuda, Y. Hay-ato, K. Huang, A. K. Ichikawa, M. Ikeda, K. Inoue,H. Ishino, Y. Itow, T. Kajita, J. Kameda, Y. Kishi-moto, M. Koga, Y. Koshio, K. P. Lee, A. Minamino,M. Miura, S. Moriyama, M. Nakahata, K. Nakamura,T. Nakaya, S. Nakayama, K. Nishijima, Y. Nishimura,Y. Obayashi, K. Okumura, M. Sakuda, H. Sekiya,M. Shiozawa, A. T. Suzuki, Y. Suzuki, A. Takeda,Y. Takeuchi, H. K. M. Tanaka, S. Tasaka, T. Tomura,M. R. Vagins, J. Wang, and M. Yokoyama, ArXiv e-prints (2011), arXiv:1109.3262 [hep-ex].[72] Hyper-Kamiokande proto-collaboration, :, K. Abe,K. Abe, H. Aihara, A. Aimi, R. Akutsu, C. Andreopou-los, I. Anghel, L. H. V. Anthony, and et al., ArXive-prints (2016), arXiv:1611.06118 [hep-ex].[73] J. Ahrens, X. Bai, G. Barouch, S. W. Barwick, R. C.Bay, T. Becka, K.-H. Becker, D. Bertrand, A. Biron,J. Booth, O. Botner, A. Bouchta, M. M. Boyce,S. Carius, A. Chen, D. Chirkin, J. Conrad, J. Cooley,C. G. S. Costa, D. F. Cowen, E. Dalberg, T. DeYoung,P. Desiati, J.-P. Dewulf, P. Doksus, J. Edsj¨o, P. Ek-str¨om, T. Feser, M. Gaug, A. Goldschmidt, A. Hall-gren, F. Halzen, K. Hanson, R. Hardtke, M. Hellwig,H. Heukenkamp, G. C. Hill, P. O. Hulth, S. Hundert-mark, J. Jacobsen, A. Karle, J. Kim, B. Koci, L. K¨opke,M. Kowalski, J. I. Lamoureux, H. Leich, M. Leuthold,P. Lindahl, I. Liubarsky, P. Loaiza, D. M. Lowder,J. Madsen, P. Marciniewski, H. S. Matis, T. C. Miller,Y. Minaeva, P. Mioˇcinovi´c, P. C. Mock, R. Morse,T. Neunh¨offer, P. Niessen, D. R. Nygren, H. Ogelman,C. P´erez de los Heros, R. Porrata, P. B. Price, K. Rawl-ins, C. Reed, W. Rhode, S. Richter, J. Rodr´ıguez Mar- tino, P. Romenesko, D. Ross, H.-G. Sander, T. Schmidt,D. Schneider, R. Schwarz, A. Silvestri, M. Solarz, G. M.Spiczak, C. Spiering, N. Starinsky, D. Steele, P. Steffen,R. G. Stokstad, O. Streicher, P. Sudhoff, I. Taboada,L. Thollander, T. Thon, S. Tilav, M. Vander Donckt,C. Walck, C. Weinheimer, C. H. Wiebusch, R. Wis-chnewski, H. Wissing, K. Woschnagg, W. Wu, G. Yodh,and S. Young, Astroparticle Physics , 345 (2002),astro-ph/0105460.[74] IceCube Collaboration, R. Abbasi, Y. Abdou, T. Abu-Zayyad, M. Ackermann, J. Adams, J. A. Aguilar,M. Ahlers, M. M. Allen, D. Altmann, and et al., ArXive-prints (2011), arXiv:1108.0171 [astro-ph.HE].[75] DUNE Collaboration, R. Acciarri, M. A. Acero,M. Adamowski, C. Adams, P. Adamson, S. Adhikari,Z. Ahmad, C. H. Albright, T. Alion, and et al., ArXive-prints (2015), arXiv:1512.06148 [physics.ins-det].[76] R. Laha, J. F. Beacom, and S. K. Agarwalla, ArXive-prints (2014), arXiv:1412.8425 [hep-ph].[77] F. An, G. An, Q. An, V. Antonelli, E. Baussan, J. Bea-com, L. Bezrukov, S. Blyth, R. Brugnera, M. BuizzaAvanzini, J. Busto, A. Cabrera, H. Cai, X. Cai,A. Cammi, G. Cao, J. Cao, Y. Chang, S. Chen, S. Chen,Y. Chen, D. Chiesa, M. Clemenza, B. Clerbaux, J. Con-rad, D. D’Angelo, H. De Kerret, Z. Deng, Z. Deng,Y. Ding, Z. Djurcic, D. Dornic, M. Dracos, O. Drapier,S. Dusini, S. Dye, T. Enqvist, D. Fan, J. Fang, L. Favart,R. Ford, M. G¨oger-Neff, H. Gan, A. Garfagnini, M. Gi-ammarchi, M. Gonchar, G. Gong, H. Gong, M. Go-nin, M. Grassi, C. Grewing, M. Guan, V. Guarino,G. Guo, W. Guo, X.-H. Guo, C. Hagner, R. Han,M. He, Y. Heng, Y. Hsiung, J. Hu, S. Hu, T. Hu,H. Huang, X. Huang, L. Huo, A. Ioannisian, M. Jeitler,X. Ji, X. Jiang, C. Jollet, L. Kang, M. Karagounis,N. Kazarian, Z. Krumshteyn, A. Kruth, P. Kuusiniemi,T. Lachenmaier, R. Leitner, C. Li, J. Li, W. Li, W. Li,X. Li, X. Li, Y. Li, Y. Li, Z.-B. Li, H. Liang, G.-L.Lin, T. Lin, Y.-H. Lin, J. Ling, I. Lippi, D. Liu, H. Liu,H. Liu, J. Liu, J. Liu, J. Liu, Q. Liu, S. Liu, S. Liu,P. Lombardi, Y. Long, H. Lu, J. Lu, J. Lu, J. Lu,B. Lubsandorzhiev, L. Ludhova, S. Luo, V. Lyashuk,R. M¨ollenberg, X. Ma, F. Mantovani, Y. Mao, S. M.Mari, W. F. McDonough, G. Meng, A. Meregaglia,E. Meroni, M. Mezzetto, L. Miramonti, T. Mueller,D. Naumov, L. Oberauer, J. P. Ochoa-Ricoux, A. Ol-shevskiy, F. Ortica, A. Paoloni, H. Peng, J.-C. Peng,E. Previtali, M. Qi, S. Qian, X. Qian, Y. Qian, Z. Qin,G. Raffelt, G. Ranucci, B. Ricci, M. Robens, A. Ro-mani, X. Ruan, X. Ruan, G. Salamanna, M. Shaevitz,V. Sinev, C. Sirignano, M. Sisti, O. Smirnov, M. Soiron,A. Stahl, L. Stanco, J. Steinmann, X. Sun, Y. Sun,D. Taichenachev, J. Tang, I. Tkachev, W. Trzaska,S. van Waasen, C. Volpe, V. Vorobel, L. Votano, C.-H. Wang, G. Wang, H. Wang, M. Wang, R. Wang,S. Wang, W. Wang, Y. Wang, Y. Wang, Y. Wang,Z. Wang, Z. Wang, Z. Wang, Z. Wang, W. Wei, L. Wen,C. Wiebusch, B. Wonsak, Q. Wu, C.-E. Wulz, M. Wurm,Y. Xi, D. Xia, Y. Xie, Z.-z. Xing, J. Xu, B. Yan,C. Yang, C. Yang, G. Yang, L. Yang, Y. Yang, Y. Yao,U. Yegin, F. Yermia, Z. You, B. Yu, C. Yu, Z. Yu, S. Za-vatarelli, L. Zhan, C. Zhang, H.-H. Zhang, J. Zhang,J. Zhang, Q. Zhang, Y.-M. Zhang, Z. Zhang, Z. Zhao,Y. Zheng, W. Zhong, G. Zhou, J. Zhou, L. Zhou,R. Zhou, S. Zhou, W. Zhou, X. Zhou, Y. Zhou, Y. Zhou, and J. Zou, Journal of Physics G Nuclear Physics ,030401 (2016), arXiv:1507.05613 [physics.ins-det].[78] J.-S. Lu, Y.-F. Li, and S. Zhou, Phys. Rev. D ,023006 (2016), arXiv:1605.07803 [hep-ph].[79] J. T. Betts, Practical methods for optimal control andestimation using nonlinear programming , Vol. 19 (Siam,2010).[80] R. Kimura, Journal of Wind Engineering and IndustrialAerodynamics , 1403 (2002).[81] E. Kalnay, Atmospheric modeling, data assimilation andpredictability (Cambridge university press, 2003).[82] G. Evensen,
Data assimilation: the ensemble Kalmanfilter (Springer Science & Business Media, 2009).[83] W. G. Whartenby, J. C. Quinn, and H. D. Abarbanel,Monthly Weather Review , 2502 (2013).[84] D. Rey, J. An, J. Ye, and H. D. Abarbanel, NonlinearProcesses in Geophysics (Accepted 2016).[85] C. D. Meliza, M. Kostuk, H. Huang, A. Nogaret,D. Margoliash, and H. D. Abarbanel, Biological cy-bernetics , 495 (2014).[86] N. Kadakia, E. Armstrong, D. Breen, U. Morone,A. Daou, D. Margoliash, and H. DI Abarbanel, arXivpreprint arXiv:1608.04005 (2016).[87] D. Breen, S. Shirman, E. Armstrong, N. Kadakia, andH. Abarbanel, arXiv preprint arXiv:1608.04433 (2016).[88] B. A. Toth, M. Kostuk, C. D. Meliza, D. Margoliash,and H. D. Abarbanel, Biological cybernetics , 217(2011).[89] A. Tarantola,
Inverse problem theory and methods formodel parameter estimation (siam, 2005).[90] A. W¨achter and L. T. Biegler, Mathematical Program-ming , 25 (2006).[91] A. B. Balantekin and G. M. Fuller, Physics Letters B , 195 (2000), arXiv:hep-ph/9908465.[92] D. O. Caldwell, G. M. Fuller, and Y.-Z. Qian, Phys.Rev. D , 123005 (2000), arXiv:astro-ph/9910175.[93] M. A. Rudzsky, Astrophysics and Space Science , 65(1990).[94] B. H. J. McKellar and M. J. Thomson, Phys. Rev. D , 2710 (1994).[95] L. Wolfenstein, Phys. Rev. D , 2369 (1978).[96] S. P. Mikheyev and A. Y. Smirnov, Yad. Fiz. (1985).[97] R. C. Schirato and G. M. Fuller, ArXiv Astrophysicse-prints (2002), astro-ph/0205390.[98] R. Tom`as, M. Kachelrieß, G. Raffelt, A. Dighe, H.-T.Janka, and L. Scheck, Journal of Cosmology and As-troparticle Physics , 015 (2004), astro-ph/0407132.[99] J. Xu, L.-J. Hu, R.-C. Li, X.-H. Guo, and B.-L. Young,ArXiv e-prints (2014), arXiv:1412.7240 [hep-ph].[100] G. Fogli, E. Lisi, A. Mirizzi, and D. Montanino, Jour-nal of Cosmology and Astroparticle Physics , 012 (2006).[101] R. D. Hoffman, S. E. Woosley, G. M. Fuller, and B. S.Meyer, Astrophys. J. , 478 (1996).[102] R. I. Epstein, S. A. Colgate, and W. C. Haxton, Phys.Rev. Lett. , 2038 (1988).[103] P. Banerjee, Y.-Z. Qian, A. Heger, and W. Haxton, Proceedings, 13th International Symposium on Originof Matter and Evolution of the Galaxies (OMEG2015):Beijing, China, June 24-27, 2015 , EPJ Web Conf. ,06001 (2016).[104] Private communication.[105] A. Malkus, A. Friedland, and G. C. McLaughlin, ArXive-prints (2014), arXiv:1403.5797 [hep-ph].[106] M.-R. Wu, H. Duan, and Y.-Z. Qian, Physics LettersB , 89 (2016), arXiv:1509.08975 [hep-ph].[107] D. V¨a¨an¨anen and G. C. McLaughlin, Phys. Rev. D ,105044 (2016), arXiv:1510.00751 [hep-ph].[108] A. Malkus, G. C. McLaughlin, and R. Surman, Phys.Rev. D editedby J. Carlson and M. J. Savage.[115] J. C. Quinn and H. D. Abarbanel, Quarterly Journal ofthe Royal Meteorological Society , 1855 (2010).[116] J. C. Quinn and H. D. Abarbanel, Journal of Compu-tational Physics , 8168 (2011).[117] A. Reyes, D. Lee, C. Graziani, and P. Tzeferacos, ArXive-prints (2016), arXiv:1611.08084 [physics.comp-ph].[118] A. D. Wyner, Information and Control , 51 (1978).[119] G. M. Fuller and Y.-Z. Qian, Phys. Rev. D , 023004(2006), arXiv:astro-ph/0505240.[120] H. A. Bethe, Reviews of Modern Physics , 801 (1990).[121] A. Burrows, J. Hayes, and B. A. Fryxell, Astrophys. J. , 830 (1995), astro-ph/9506061.[122] H.-T. Janka, K. Langanke, A. Marek, G. Mart´ınez-Pinedo, and B. M¨uller, Physics Reports , 38 (2007),astro-ph/0612072.[123] J. F. Cherry, J. Carlson, A. Friedland, G. M. Fuller,and A. Vlasenko, Physical Review Letters , 261104(2012), arXiv:1203.1607 [hep-ph].[124] J. F. Cherry, J. Carlson, A. Friedland, G. M. Fuller,and A. Vlasenko, Phys. Rev. D87