{'fraction_non_alphanumeric': 0.03968238963985254, 'fraction_numerical': 0.01105964646942055, 'mean_word_length': 4.778457504915884, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Extended evolution equations for neutrino propagation in astrophysical and cosmological environments

We derive the evolution equations for a system of neutrinos interacting among themselves and with a matter background, based upon the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy. This theoretical framework gives an (unclosed) set of first-order coupled integro-differential equations governing the evolution of the reduced density matrices. By employing the hierarchy, we first rederive the mean-field evolution equations for the neutrino one-body density matrix associated with a system of neutrinos and antineutrinos interacting with matter and with an anisotropic neutrino background. Then, we derive extended evolution equations to determine neutrino flavor conversion beyond the commonly used mean-field approximation. To this aim we include neutrino-antineutrino pairing correlations to the two-body density matrix. The inclusion of these new contributions leads to an extended evolution equation for the normal neutrino density and to an equation for the abnormal one involving the pairing mean-field. We discuss the possible impact of neutrino-antineutrino correlations on neutrino flavor conversion in the astrophysical and cosmological environments, and possibly upon the supernova dynamics. Our results can be easily generalized to an arbitrary number of neutrino families.

We derive the evolution equations for a system of neutrinos interacting among themselves and with a matter background, based upon the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy. This theoretical framework gives an (unclosed) set of first-order coupled integro-differential equations governing the evolution of the reduced density matrices. By employing the hierarchy, we first rederive the mean-field evolution equations for the neutrino one-body density matrix associated with a system of neutrinos and antineutrinos interacting with matter and with an anisotropic neutrino background. Then, we derive extended evolution equations to determine neutrino flavor conversion beyond the commonly used mean-field approximation. To this aim we include neutrino-antineutrino pairing correlations to the two-body density matrix. The inclusion of these new contributions leads to an extended evolution equation for the normal neutrino density and to an equation for the abnormal one involving the pairing mean-field. We discuss the possible impact of neutrino-antineutrino correlations on neutrino flavor conversion in the astrophysical and cosmological environments, and possibly upon the supernova dynamics. Our results can be easily generalized to an arbitrary number of neutrino families.

I. INTRODUCTION
Oscillations between quantum states is a widespread phenomenon, appearing in different physical contexts, like the Rabi oscillations in optics, or the K 0 −K 0 and neutrino oscillations in particle physics. The propagation in a medium can produce resonant conversion between quantum states, such as when neutrinos change their flavor while traveling in a star, or when photons modify their polarization in a birefringent medium. The basic oscillation phenomenon can be modified sometimes in surprising ways, if the interaction with a medium introduces complexity. The study of how neutrinos change their flavor in stellar environments and in the early universe has uncovered numerous such examples. Complexity arises from the non-linear and the many-body character of the problem.
Neutrinos are elementary particles having non-zero mixings, as first conjectured by Pontecorvo [1] and discovered in 1998 by Super-Kamiokande [2]. Numerous experiments have contributed to the measurement of the neutrino mixing angles and squared-mass differences that shape the way they change their flavor while traveling [3]. In particular, precisely known are the neutrino mixing angles of the Maki-Nakagawa-Sakata-Pontecorvo (MNSP) matrix that relates the interaction (flavor) to the propagation (mass) eigenstate basis [4], the squaredmass difference values and one sign, while the other sign remains unknown (the hierarchy problem). Addressing the question of the value of the (Dirac or Majorana) CP violating phases is one of the major future goals, jointly with the determination of the neutrino absolute mass and of the neutrino (Dirac or Majorana) nature [5].
Astrophysical and cosmological environments produce copious amounts of neutrinos. Therefore, a precise knowledge of neutrino flavor conversion in media is required, to assess e.g. the neutrino impact on the supernova dynamics and on (stellar or primordial) nucleosynthesis processes, to interpret the signal associated with solar neutrinos, to predict the one produced by core-collapse supernovae, or to understand how neutri-nos change their flavor while traversing the Earth. It is experimentally established that the origin of the solar neutrino deficit is a resonant flavor conversion induced by neutrinos interacting with matter. This is the well known Mikheev-Smirnov-Wolfenstein (MSW) effect [6,7]. More precisely, the MSW effect produces a deficit of the high energy ( 8 B) neutrinos, while averaged vacuum oscillations account for the one of low energy ( 7 Be, pp and pep) ν. While the MSW is the reference phenomenon to understand how neutrinos change their flavor in media, various other phenomena impacting flavor conversion can occur, depending on the specific environment under consideration. For example, Pantaleone has first pointed out the presence of a non-linear refractive index due to the neutrino interaction with other neutrinos whenever the neutrino number density is large [8]. This is relevant for neutrino evolution in the early universe, in core-collapse supernovae, and for low energy neutrinos produced in accretion-disk black-hole scenarios. Indeed, simulations implementing the neutrino-neutrino interaction show the emergence of new phenomena as pointed out by Samuel [9], that can be interpreted as the synchronization of effective spins [10], a flavor [11] or gyroscopic pendulum [12], or a magnetic resonance phenomenon [13] (see e.g. [14] for a review). Features associated with the explosion dynamics of a core-collapse supernova, such as the presence of shock-waves and of turbulence, produce new interference phenomena, like multiple MSW resonances [15] and eventually depolarization [16,17]. The spectral and location changes of the supernova neutrino fluxes impact supernova observations and are currently being investigated. In the cosmological context, neutrino flavor conversion is also important at the epoch of Big-bang nucleosynthesis (see e.g. [18,19] for a review). Numerous works have investigated for example the effects of mixings among active flavors [20][21][22], with a possible non-zero leptonic CP violating phase [23,24] or between active and sterile neutrinos [24][25][26][27] on the primordial element abundance(s).
Numerous works have aimed at formulating theoretically the equations of motion that describe the evolution of particles with mixings in a medium. The MSW effect is usually accounted for, by using an effective Hamiltonian that is linear in the weak coupling constant, and depends upon the matter number density [5]. In [8] a generalization of the neutrino evolution equations is made in a similar way, to implement interaction of neutrinos with themselves. Pantaleone already emphasizes the complexity inherent with the non-linearity and the many-body character of the problem. Ref. [9] has given a mean-field equation including such interaction terms. Reference [28] has first derived evolution equations for neutrino density matrices including such interaction terms. Refs. [29] and [30] have derived neutrino evolution equations beyond the mean-field approximation and including two-body collision terms with the "molecular chaos" assumption that neglects the building up of correlations in the collision term. Such Boltzmann equations are formulated in terms of matrices of neutrino densities. Ref. [32] has also derived the neutrino Boltzmann evolution equation using first quantization, having in mind the case of the early universe. Ref. [33] has generalized the equations of [30] to the three flavor case (without the collision terms) and made the angular dependence more explicit. Many-body aspects, and the possible breakdown of the one-body description, have been discussed in Refs. [34][35][36], using the spin-spin analogy in simplified models. Liouville equations for neutrino distribution matrices are derived in [37]. An algebraic approach to the neutrino propagation in media is given in [38] and the evolution of the manybody problem is formulated as a coherent-state path integral. This allows one, in particular, to calculate corrections to the mean-field equations as a determinant coming from the path integral [38]. The algebraic based formulation is further employed in [39], where the neutrino Hamiltonian of the many-body system (with mixings and the neutrino-neutrino contribution but without the matter term) is put in connection with the (reduced) BCS pairing Hamiltonian describing superconductivity. It is pointed out that the corresponding constants of motion show the exact solvability of the problem.
The present work formulates the problem of the neutrino evolution, in terms of reduced density matrices, using the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy [40][41][42][43]. This replaces the Liouville Von-Neumann equation for the many-body density matrix describing the full many-body problem, by an (unclosed) set of coupled integro-differential equations for the reduced density matrices. Using this theoretical framework, we derive the equations of motion for a system of neutrinos, traveling in a medium made of ordinary matter, and interacting with each other. First we show that truncating the hierarchy at lowest order produces the mean-field equations for the one-body density matrix, commonly used in the literature, to implement the neutrino interaction with matter and with neutrinos. Next, we focus on the neutrino evolution description beyond the mean-field approximation and include for the first time neutrino-antineutrino pairing correlations to the two-body density matrix. Such bilinear products have been neglected so far, based upon the argument that their expectation value over free states typically oscillate fast around zero. However, since the neutrino evolution equations are often non-linear, it is worthwhile to investigate their possible impact on neutrino flavor evolution in a medium. We show that the inclusion of neutrinoantineutrino pairing correlations leads to extended timedependent mean-field equations both for the normal and for an abnormal neutrino density matrix. Finally we conclude by discussing the possible implications of these contributions for neutrino flavor conversion in the astrophysical and cosmological environments.
The manuscript is organised as follows. Section II presents the theoretical framework of the BBGKY hierarchy, its lowest order truncation that furnishes the mean-field approximation and the evolution equation for the two-body correlation function. In Section III we rederive the mean-field neutrino evolution equations including both the coupling to matter (the MSW contribution) and to neutrino themselves (the νν interaction term). We discuss the relationship with the equations commonly used in the literature. Section IV introduces the contribution from νν pairing correlations to the two-body correlation function. Our extended mean-field equations for the normal and abnormal neutrino density matrices, involving the normal and pairing mean-fields, are presented. Section V includes a discussion and our conclusions.

II. THE THEORETICAL FRAMEWORK
A. The Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy In numerous contexts one is interested in determining the dynamics of a system made up of N interacting particles, such a gas of weakly interacting neutrinos, or an ensemble of strongly interacting nucleons in a nucleus or in a collision among nuclei. The Hamiltonian for such a system of N-particles, interacting through a two-body interaction, readŝ comprising the one-bodyĤ 0 kinetic and the two-bodŷ V interaction terms. The k, k ′ indices run over singleparticle quantum states, identified by single-particle properties like momentum, flavor, helicity, isospin, etc...
The system's evolution is determined by solving the Schrödinger equation for the many-body quantum state |ψ(t) in case of a pure state or, more generally, the Liouville Von-Neumann equation for the many-body density matrixD: withĤ given by Eq.(1) (here we takeh = c = 1). In the BBGKY hierarchy 1 theoretical framework [40][41][42][43], one introduces the s-reduced density matrixρ 1...s defined aŝ tr s+1 indicating that we are tracing over the s+1 particle, and replaces Eq.(2) forD, by an unclosed equation for ρ 1...s : where 2 denotes the s-body density matrix components 3 . The a s and a † s correspond to the particle annihilation and creation operators for a particle in the quantum state s, respectively. In particular, the one-body and two-body matrix elements components are: In Eq. = k V (k, s + 1) with k=1. . . s. This equation is unclosed since the ρ 1...s evolution is coupled to the one of the (s+1)-reduced density 1 See e.g. Ref. [44]. 2 Note that we do not write explicitly "hat" over the creation and annihilation operators through the whole manuscript, not to overburden the text. 3 Note that we denote withρ 1...s the operators, while we indicate the density matrix components with ρ 1...s or ρ(1 . . . s). ρ 1...s+1 , via the two-body interaction. The BBGKY hierarchy (4) can easily be deduced by applying successive traces to Eq.(2) and using the property More explicitly Eq.(4) can be written as a hierarchy of equations of motion for the 1-body ρ 1 to the s-reduced ρ 1...s density matrix components 4 : (9) Solving Eq.(9) is completely equivalent to determining the exact evolution forD Eq.(2). The advantage of the BBGKY framework is that it furnishes a hierarchy of evolution equations for the reduced density matrices of increasing order, so that one can test different approximations, by going at a higher truncation level in the hierarchy.

The mean-field approximation for the evolution equations
Let us consider the first equation of the BBGKY hierarchy: One can separate the correlated from the uncorrelated 5 contribution of the two-body density matrix 6 : where c 12 is the two-body correlation function. If one deals with identical fermionic particles the uncorrelated contribution ρ 1 ρ 2 has to be replaced by ρ 1 ρ 2 (1 − P 12 ), with P 12 being the operator that exchanges particle 1 with particle 2, to properly account for the antisymmetrization. Inserting Eq. (11) in Eq.(10), one obtains Here, since no approximation is made, the dynamical equation for the one-body density is exact. Now, neglecting the correlated contribution to the twobody correlation function, one gets or, equivalently: with h 1 (ρ) = H 0 (1) + Γ 1 (ρ) and Γ 1 (ρ) = tr 2 (V (1, 2)ρ 2 ) being the mean-field acting on particle 1. This is the so-called mean-field approximation.
Writing such an equation more explicitly, it reads The mean-field potential is built up from a complete set of one-body density matrix components for particle 2 ρ 2,nm , each contributing with the matrix element 7 v (im,jn) = im|V 12 |jn , with jn (im) incoming (outgoing) single-particle states. From Eq.(16) one can see the dependence of our mean-field on the one-body density associated with particle 2, while in some cases the interaction itself might also have an explicit dependence on ρ. To solve Eq.(15) one has to assign the state of the many-body system at initial time, which can be either a correlated state, or a product of independent singleparticle states. In the latter case, the condition inherent to Eq.(15), i.e. c 12 = 0, ensures that it stays as such at any time.
It is worthwhile to mention that a first-order evolution equation for the one-body density matrix associated with a givenD can also be obtained by applying the Ehrenfest theorem: 7 Note that in case of identical particles the matrix elements are antisymmetrised, i.e.ṽ (im,jn) = im|V 12 |jn − im|V 12 |nj .
with ρ 1,ij = a † j a i andĤ given by Eq.(1) (in second quantization). In particular, Eq.(15) is recovered when neglecting the correlated contribution to the two-body density matrix [45].

Beyond the mean-field approximation
Our main goal will be to discuss contributions beyond the mean-field approximation given by Eq. (14), to the evolution equations for a system of relativistic neutrinos that interact among themselves and with matter. To this aim a useful reformulation of Eq.(9) is given by a hierarchy of evolution equations for the correlation functions (details of the demonstration can be found in Ref. [46]), where only linked terms are shown to remain. Such a reformulation has the advantage that higher-order contributions are expected to decrease with increasing rank [46]. In this context, the mean-field equation (14) is unchanged; while one gets for the two-body correlation function [46] with h 1 (ρ) (h 2 (ρ)) the mean-fields acting on particles 1 (2) respectively and P 13 (P 23 ) is the operator that exchange particle 1 (2) with 3 (see Appendix A for an explicit formulation of Eq. (18)). Such an equation contains three main contributions coming from two-body interactions. Retaining the second and third term on the r.h.s. of Eq. (18), and making the "molecular chaos" assumption, that the build up of correlations due to collisions is negligible, one gets a collision term with incoming and outgoing particles described by free particle states (see e.g. [47]). Such a term gives rise to a Boltzmann equation 8 . The (1 − ρ 1 )(1 − ρ 2 ) factor 9 ensures the appropri- 8 Note that a Boltzmann equation for a system of neutrinos and antineutrinos is derived in [29,30] and in [32]. Such evolution equations include a collision term as required in the context of the early Universe. In core-collapse supernova simulations, full transport equations for neutrinos (but without the inclusion of mixings) are usually employed in the dense region where neutrinos are trapped (see e.g. [49][50][51][52][53]). 9 This factor as well as the (1 − ρ 1 − ρ 2 ) one come from the linked contribution of the trace term over the third particle.
ate statistics (no contribution if the final single-particle states are already occupied). The fourth and fifth terms on the r.h.s. of Eq.(18) implements contributions from the correlated part of the two-body correlation function. It is on this term that we will focus later on. In particular we will consider cases involving both particles and antiparticles. In such systems, one might have possible contributions from the expectation values of the product operators of the type a k b l and a † k b † l . Such bilinear products include particle-antiparticle correlations that can be seen as correlations of the pairing type. These terms have been neglected so far, based on the argument their expectation value (over free states) typically oscillate fast around zero, a † ( p, t)b † ( p ′ , t) (see e.g. [30,37]). Finally, the three-body terms give the contribution from the two-body interaction among three particles, obtained by tracing over the third particle. We will neglect these higher order correlations here. Since, in this work, we focus on the inclusion of νν contributions to c 12 , we will not consider the collision and the three-body terms. Our evolution equation for the two-body correlation function is

B. The application to neutrinos
We are here mainly interested in discussing the neutrino evolution and flavor conversion when neutrinos propagate in an astrophysical environment, or in the early universe. While the BBGKY hierarchy is generally employed for systems of interacting particles without mixings, we here consider that the density matrices in Eq.(4) are associated with mixed particles. Note that a "matrix of densities", generalizing the usual occupation numbers, is commonly used in the literature (see e.g. [29,30]). More explicitly, in the three flavors case the neutrino density matrix reads The off-diagonal (or coherent) terms are non-zero to encode the presence of neutrino mixings. The neutrino occupation number for a given ν α,i flavor state is given by the diagonal element of the density matrix a † να,i a να,i , with N να = i a † να,i a να,i the total occupation number. The definition in Eq. (20) can easily be extended to the case of an arbitrary neutrino families, to account for the presence of both sterile and active neutrinos.
Since the systems we are interested in involve both particles and antiparticles, one introduces a density matrix ρ ν , in a way analogous to Eq. (20), but replacing the particle operators a † , a with the antiparticle b † , b ones [30]. For the sake of clarity concerning the convention used in the present work, we give its explicit expression: To implement the antiparticle degrees of freedom, we consider the usual expansion of the neutrino fields in the Schrödinger picture The summation is over the helicity h states, the integration over the momenta p of the (anti)particles and u p,h (v p,h ) are the usual Dirac spinors, E p being neutrino energy. Note that in such an expression the particle creation and annihilation operators are associated with states of a given mass for which the usual anticommutation rules are properly defined. In the present manuscript we employ 10 and The single-particle states associated with neutrino mass eigenstates are with | being the vacuum state defined by a m | = 0. The flavor eigenstates are related to the mass eigenstates 10 Similar relations hold for the antiparticle annihilation b( p, h) and through |ν α = i U * αi |ν i (i and α correspond to an arbitrary number of neutrino families), where U is the Maki-Nakagawa-Sakata-Pontecorvo (MNSP) unitary matrix. In three flavors, the MNSP matrix depends upon three neutrino mixing angles that have been measured, one Dirac and two Majorana unknown phases [3]. Since the appropriate anticommutation rules only hold for the mass eigenstates operators, Eqs. (23)(24) require special attention (see e.g. Ref. [37] for a discussion). This subtlety is sometimes avoided by writing the fields Eq. (22) for massless neutrinos, as done e.g. [30]. As discussed for example in Refs. [37,54], it is not possible to rigorously build up a Fock space for flavor states since the neutrino flavor creation and annihilation operators do not satisfy the canonical anticommutation rules Eqs. (23)(24). Indeed, as shown in Ref. [54], one can define an approximate Fock space by introducing neutrino "weak-particle states" that depend upon the specific weak process under consideration. However, in the limit of relativistic neutrinos, the anticommutation relations (23-24) also approximately hold for flavor states. Since this is a good approximation for our cases of interest (solar and supernova neutrinos, cosmological neutrinos at the epoch of big-bang nucleosynthesis), we make the assumption that an approximate Fock space for our flavor states can be built. It is worthwhile to mention that extending our results without making this approximation does not introduce extra conceptual difficulties. In much the same way as done in Ref. [54] for the muon decay case, but in the mean-fields expressions, one should retain an explicit dependence on the mixing matrix elements in the calculation of the interaction matrix elements, depending on the specific weak process under consideration 11 .
Finally, it is worthwhile to remind that, while the BBGKY theoretical framework is based upon the density matrix formalism, the evolution equations (9) can be formulated as a hierarchy for many-body Greens' functions as done in Ref. [48]. In particular, in the equal time limit, when only linked contributions are retained, these two formalisms are completely equivalent. Note that while historically BBGKY was developed to describe the evolution of nonrelativistic systems of N-particles, the hierarchy applies to a system of relativistic particles (as of interest here). In this case Eqs. (9) are replaced by an infinite set of equations.

III. NEUTRINOS EVOLVING IN A MEDIUM IN THE MEAN-FIELD APPROXIMATION
We now take the example of neutrinos interacting with the electrons, protons and neutrons composing a medium to show how the formalism presented in Section II can be used to rigorously derive well known neutrino evolution equations. We just sketch the main lines of the derivation 12 .
The assumption that is usually made is that our system of neutrinos interacting with a medium can be described at lowest order as a system of independent particles, so that implicitly the problem reduces to following the evolution of a single-particle at a time and calculating the evolution of the associated one-body density matrix. Let us consider the BBGKY hierarchy truncated in the mean-field approximation given by Eq. (14). The Hamiltonian Eq.(1) for our case of interest in the flavor basis reads where U is the MNSP unitary matrix relating the neutrino flavor basis to the mass eigenstate basis (with eigenenergies E i ). The first H m = diag(E i ) contribution is the propagation term, while the second one H int corresponds to the two-body interaction between a neutrino and another particle.
A. The mean-field associated with neutrino interaction with matter (MSW case) To follow the one-body density matrix evolution given by Eq. (17), one needs to determine the mean-field (16) created by the background particles and acting on the "test" neutrino. The interaction term H int corresponds to the charged-or neutral-current Hamiltonian describing the neutrino interaction with the medium. We take the example of the charged-current interaction on electrons, where H int is given by 13 12 A different derivation of the mean-field equations accounting for the neutrino interaction with a medium is given e.g. in [5]. 13 The effective low energy approximation is sufficient for the applications envisaged.
FIG. 1. Neutrino interaction with electrons and the corresponding mean-field associated with the electron background. In the evolution equations, the mean-field acting on a single neutrino state is build up from the summation of the electron single-particle states (see Eqs. (16), (28) and (32)).
where G F is the Fermi coupling constant and with the fields φ given by Eq. (22). This requires calculating the matrix elements for neutrino-electron scattering v νe,e im,jn = ν ei e m |H CC |ν ej e n with the interaction Hamiltonian (27). One gets The summation over the m, n single-particle states in Eq.(16) becomes here a sum over the electron helicity states and an integration over momenta ( Figure 1).
The key quantity to define is the expectation value of the a † e ( p, h)a e ( p ′ , h ′ ) operator, over the electron background 14 characteristic of the considered medium 15 : The assumption of a homogeneous and unpolarized medium corresponds to: Equation (29) constitutes for example a good approximation for the case of the Sun. Using Eqs. (29)(30) the total electron number is where the factor of 2 comes from the summation over the electron helicity states, and V is the quantization volume. In Eq.(28), by tracing over the spinors, implementing that neutrinos have only one possible helicity state, and assuming the electron background is homogeneous, unpolarized and isotropic, one gets for the meanfield Eq.(16) where n e = N e /V is the electron number density. In Eq.(32) the δ-function ensures that the momentum of the neutrino propagating in the electron medium is unchanged, as a consequence of homogeneity. The prefactor (2π) 3 2E k are present because of the chosen normalization of the (anti)particle anticommutation relations Eqs. (23)(24). Our result (32) is the mean-field contribution to the neutrino Hamiltonian, corresponding to neutrino interaction with matter, that is linear in the weak coupling constant and depends upon the number densities of the particles composing the medium. We find the well known low energy effective Hamiltonian (see e.g. [5]) that can give rise, if neutrinos propagate adiabatically in a medium, to a resonant 16 flavor conversion -the Mikheev-Smirnov-Wolfenstein effect [6,7]. The procedure just outlined can be applied to the neutral-current ν e scattering on electrons, to the chargedand neutral-currentν e scattering on electrons, of ν e or ν e on positrons, as well as neutral-current scattering on protons and neutrons, giving the expected results. In particular, it is immediate to show that, when dealing with antiparticles, the associated mean-field Γ depends upon matrix elements that involve the b † and b operators (instead of the particle operators), introducing a minus sign. For example, for the case ofν e evolving in an electron background, one obtains the expected result 16 The occurence of such a resonant phenomenon depends upon neutrino properties (energies, mixing angles, squared-mass differences value and sign) and the specific matter number density profile for the system under consideration.

B. The mean-field associated with neutrino self-interactions
The second case we are going to consider, in the meanfield approximation, is when the (anti)neutrino is evolving in a background of ν andν. This case is of interest since recent studies have shown the important role of the neutrino-neutrino interaction for instance in a corecollapse supernova, and the variety of new flavor conversion phenomena that can arise, when implementing this contribution in the neutrino Hamiltonian (see e.g. [14] for a review). A series of works have discussed the neutrino evolution equation in presence of such terms [8, 9, 30-33, 38, 39]. Note that the role of such terms was first pointed out in the context of the early universe [55]. Here we just sketch the derivation of these equations, following the same procedure as for the MSW case (more details are given in Appendix B).
The neutrino-neutrino mean-field depends this time on the effective low energy neutral-current Hamiltonian: with ν y = ν e , ν µ or ν τ . From Eq.(16) one gets for the mean-field Γ να,ν β (ρ ν ) with ν α , ν β = ν e , ν µ , ν τ (or the corresponding antineutrinos): and requiring the expectation value correspond to a homogeneous and unpolarized system: The key difference with the case of the electron background is that the quantity ρ ν β ,να p ′ h ′ , ph has diagonal and off-diagonal terms ( Figure 2). As a consequence the mean-field acting on a ν e for example has both diagonal Γ να,να (ρ ν ) and off-diagonal Γ να,ν β (ρ ν ) contributions. One recognizes in the mean-field contributions that the diagonal one is associated with forward scattering without flavor exchange, while the off-diagonal contribution is associated with forward scattering "with flavor exchange" (the Pantaleone diagram [8]). We emphasize that in our derivation we do not directly start with an effective Hamiltonian having these terms, as done e.g. in Ref. [38]. Here these contributions naturally arise when considering the weak interaction term, using the density matrix Eq.(20) and calculating Γ (or, in other words, closing the loop). By implementing that neutrinos have one helicity state only, tracing over the spinors, in case the neutrino background is homogeneous and anisotropic, one gets for the off-diagonal contribution the following expression: withˆ p = p/| p| andˆ k = k/| k|. A similar expression holds for the diagonal one if ν β = ν α , except for an extra factor of 2. In the case of an isotropic medium the angular term cosθ k p =ˆ k ·ˆ p averages out, and one recovers the same result as Eq. (32). If the background includes antineutrinos as well, one needs to add the contribution fromρ * , but with a minus sign. By adding up the two contributions the total neutrino-neutrino mean-field reads 17 : where here ρ να, p (ρ * να, p ) stand for the density matrix Eq.(20) (Eq. (21)), ν α refers to a neutrino that is initially born in the α flavor. In fact, one has to sum over all neutrino flavors present in the system. Note that in Eq.(39) the trace term, tr(ρ να, p −ρ * να, p ), has been subtracted.

C. Neutrino evolution equations in the mean-field approximation
With the results of Eqs. (32) and (39), the mean-field equations (15)(16) for the density matrix (20) describing the neutrino evolution in a medium becomes explicitly: with

0). A similar equation holds if an antineutrino is traveling instead of a neutrino:
17 Note that the contribution coming for the antineutrino density matrix isρ * if one employs the definitionρ ij = b † j b i (see Appendix B). If one defines such a quantity asρ ij = b † i b j then the contribution appearing in the evolution equations depends upon ρ, instead ofρ * .
Note that in Eq.(40-43) ρ andρ have two indices in flavor, but only one in momentum as a consequence of homogeneity. In case the quantum state at initial time is an independent particle state, the one-body density matrices (20)(21) can be replaced by single-particle onebody densities. In the neutrino case, the diagonal elements ρ i directly give the neutrino survival probabilities |ν i | 2 ν i being the neutrino amplitude for flavor i, while the off-diagonal ones ρ ij are the mixing terms ν i ν j (and similarly for antineutrinos). We conclude here by emphasizing that the evolution equations (40)(41)(42)(43) we find are in agreement with those of Refs. [30,31,38], commonly used in the investigation of solar, of supernova neutrinos and of the low energy neutrinos in accretion-disks black-hole scenarios.

A. The extended neutrino dynamical equations
Our main goal is to obtain the evolution equations, beyond the mean-field approximation, for a system of neutrinos and antineutrinos evolving in an environment, taking into account possible νν pairing correlations. The latter correspond to the following linked contribution to the two-body density matrix Eq. (11): where the quantities κ αβ and κ * α ′ β ′ are called the abnormal densities. Imposing that the pair products conserve individual lepton numbers, only product of operators associated with the same flavor, such as a να b να , and the corresponding hermitian conjugates, are admitted. Since neutrinos have mixings, one can have contributions from products involving neutrino-antineutrino 18 Note that, according to our definition forρ Eq.(21), antineutrinos do not transform the same way as neutrinos under U .
pairs with different flavors. More generally, e.g. if total lepton-number is not conserved one could have correlations associated with the pair product operators like a να a να (and similarly for antineutrinos). We note that the expectation values 19 a † a , b † b , ba and a † b † naturally appear as components of the field correlation function φφ . So far, neutrino-antineutrino correlations have been neglected 20 (see e.g. Refs. [29,30,37]). For example, in the formulation of Ref. [37], they correspond to the rapidly oscillating cross terms between the positive-and negative-frequency parts of the quantum density function iΓ lm ij = N ν l i (y)ν m j (z) , where the ν(y) andν(z) are the neutrino and antineutrino quantum field operators.
Two perspectives are possible to investigate the impact of the neutrino-antineutrino correlations on the neutrino evolution. The first is to assume that such correlations are non-zero at initial time, e.g. at the neutrino sphere in a supernova, from previous interactions among neutrinos in the dense supernova region, where they are trapped. The second possibility is to study whether such terms can be dynamically produced through the interactions. However this is more demanding, since one has to retain the collision term in the two-body correlation function Eq.(18) as well. In this manuscript we adopt the first perspective. We will see below that, according to our extended evolution equations if their contribution is zero at initial time it is zero at all times. Therefore, in this case, the mean-field approximation given by Eqs. (40)(41)(42)(43) is correct (if for the considered system, the collision term can also be neglected in Eq. (18)). This is the approximation that is usually implicitly made.
Let us now discuss how the neutrino evolution equation Eqs. (40)(41)(42)(43) have to be extended to implement pairing correlations between ν andν. In this case, the evolution of the system is determined by using the first two equations of the BBGKY hierarchy Eq.(9). By substituting Eq. (44) in Eq. (19) one obtains the evolution equation for the abnormal density (and its complex conjugate), while from the first equation of the BBGKY hierarchy and Eq.(11), we obtain the evolution equation for the normal density. We finally get the extended evolution equations: (45) where here the indices i, j stand for ν α , ν β ; k, l forν α ,ν β with α, β that vary over the different electron, muon and tau flavor states. For the sake of clarity, in Eqs. (45) we show explicitly the dependence on particle 1 and particle 2 of the quantities. Obviously, in our extended equations one has to determine the evolution of the normal density associated with a neutrino (particle 1) and an antineutrino (particle 2) to determine the neutrino-antineutrino pair evolution in κ and κ * . We emphasize that, the new evolution equations Eq. (45) can also be derived by using the Ehrenfest theorem and determining the first-order evolution equations for the abnormal density: and for the normal one using Eq. (17). The expression for the abnormal mean-field is and of its complex conjugate One can see that such expressions are analogous to the one for the mean-field Γ Eq. (16), but with the abnormal density κ replacing the normal density ρ, and by summing over the initial (or final) single-particle states instead of over a final and initial single-particle state. We note that Eqs. (45) can be cast in an elegant and compact matrix form: where we have introduced the generalized density R that depends upon both the normal densities for the neutrinos and the antineutrinos and the abnormal density.ν 3. The figure shows, in a pictorial way, the interaction of a neutrino with an antineutrino and the corresponding abnormal mean-field ∆ * Eq.(48) which is build up from the sum of single-particle neutrino-antineutrino states. Its expression is analogous to the one of the normal mean-field by replacing the normal density by the abnormal one κ and by summing over the final (initial in ∆) states. The contributions shown correspond to the diagonal component ∆ * να,να Eq.(56) (upper figures), and to the off-diagonal one ∆ * να,ν β Eq.(52) (lower figures). Note that the diagonal abnormal mean-field also receives a contribution from a term where neutrinos of the same flavor α run in the loop.
The generalized Hamiltonian H governing the evolution is given by It comprises the mean-field Hamiltonians for the neutrinos and the antineutrinos, as well as the abnormal meanfield.

B. The abnormal mean-field
Let us now compute the expression of the abnormal mean-field ∆. First, the expression for the mean-field h (h * ) acting on a neutrino (antineutrino) is given by Eqs. (41) and (43). To derive an explicit expression for ∆ Eq.(47) and ∆ * Eq. (48), one needs to calculate the matrix element v (ij,kl) associated with the neutral-current interaction of a neutrino with an antineutrino. One obtains the following expression for the off-diagonal components of the abnormal field (see Figure 3): The explicit expression of the abnormal field Eq.(52) depends on the properties of the background: Different assumptions can be made on the expectation value of the correlator a † να ( k, h α )b † ν β ( p, h β ) and its complex conjugate. In order to remain as general as possible, various options are considered in the following. We give the results for ∆ * να,ν β while similar expressions hold for ∆ να,ν β as well as for the diagonal contributions of the abnormal mean-field ∆ να,να and its complex conjugate ∆ * να,να (only the final results are presented here, while more details on their derivation can be found in Appendix C.). In particular, for the calculation of the latter, two contributions need to be added: as can be seen from Eq.(48). The first contribution ∆ * (eq) να,i,να,j comes from the two amplitudes for the process of neutrino-antineutrino scattering and annihilation for equal flavors (see Figure 3) that summed give a factor of 2 times an expression similar to the one shown in Eq.(52): The second contribution to ∆ * (uneq) να,να comes from different flavors than α running in the loop (Fig.3). Let us now perform an explicit calculation of the pairing mean-field, depending on the background properties.
The requirement that the background through which neutrinos are traveling is homogeneous, is fulfilled, if we impose that the neutrino and the antineutrino have opposite momentum in the expectation values of the neutrinoantineutrino pair operators. This corresponds to: and we take h α = −1 and h β = +1. Note that the quantities κ ναν β * p (or κ νανα * p ) are related to the neutrinoantineutrino pair number density: in a way analogous e.g. to the electron number densities (see Eq.(31)). From Eq.(52), by employing the procedure discussed in Section III.A and III.B for the normal fields, one obtains that the abnormal mean-field in spherical coordinates reads . Note that in case there is cylindrical symmetry only the last term remains: The expression Eq. (59) gives in cartesian coordinates is So far we have neglected the neutrino masses so that ν andν have definite helicities. Let us now consider the case that e.g. h α = h β = +1 or −1: While such a contribution is possible, we expect the corresponding abnormal field to be suppressed by a factor (m/E) 2 , E being the neutrino energy. For this specific derivation, since one is sensitive to the neutrino masses, the calculation has to be performed by working in the neutrino mass eigenstate basis. Following this procedure, one finds where here the notation E i ′ k ′ indicates the energy of a neutrino with mass m i ′ and momentum k ′ .

C. The three-flavor neutrino case
We conclude by giving the expressions of the extended dynamical equations for three neutrino flavors. In this case the abnormal density is defined as and its complex conjugate as Equation (45) gives to the following set of equations 21 :  (42)) for neutrinos (antineutrinos), while the κ ν and κ * ν are given by Eq.(66-67). In order to determine the neutrino flavor evolution, one needs to assign initial conditions for the normal and abnormal densities, i.e. ρ 0 = ρ(t = 0),ρ 0 =ρ(t = 0) and κ 0 = κ(t = 0) for the extended mean-field equations Eq. (49) with (50)(51). In the case of core-collapse supernova neutrinos, such an assignment could be done by extracting the relevant information from realistic simulations of the dense region where neutrinos are trapped that would include the relevant correlations. Otherwise, a guess that one can make is to assume that the system is approximately described at initial time by a "stationary state" of the extended mean-field Hamiltonian Eq.(51); while the system is driven out of such a solution at later times [59].
Finally, we would like to point out that, in the application of our equations to the case of supernova neutrinos, the conditions Eqs. (57) and (64) cannot be met, strictly speaking. In fact, in the region outside the neutrinosphere, where the density becomes low enough, all neutrinos and antineutrinos start free streaming. However, in the region just before the neutrino-sphere and at its boundary the presence of collisions might produce nonzero neutrino-antineutrino correlations as we consider here. It is this intermediate region, that we are interested in, where one goes from the regime where neutrinos are trapped requiring a Boltzmann treatment, and the one where only forward scattering becomes relevant and the mean-field approximation given by Eqs. (40)(41)(42)(43) constitutes a good approximation. The presence of ν −ν correlations in this transition region might indeed play a role from the point of view of the supernova dynamics, or eventually modify neutrino flavor conversion when neutrinos start free streaming.

V. DISCUSSION AND CONCLUSIONS
In the present work, we have derived the neutrino evolution equations in a medium composed of matter and eventually of neutrinos, having in mind environments that produce copious neutrino amounts such as the Sun, core-collapse supernovae, accretion disks around black holes or the early universe. In particular, observations of such neutrinos require a precise understanding of the corresponding neutrino number fluxes and spectra, their modifications through such media, and their specific signatures both in solar and supernova neutrino detectors, as well as on stellar, or primordial, nucleosynthesis abundances.
In the astrophysical environments the neutrino flavor evolution is essentially treated in the mean-field approximation. While this is a good description in the case of solar neutrinos, the investigation of neutrino flavor evolution in core-collapse supernovae might require going beyond. For example, in the transition between the dense region where neutrinos are trapped, and the neutrinosphere where they start free streaming, many-body correlations might have sizeable effects and impact flavor conversion. The physical argument to support the meanfield assumption is that the MSW resonances occur in the outer layers of the star, because of the large neutrino mass-squared differences, and that the MSW effect is well accounted for in the mean-field approximation. However, this picture -well separating the region where flavor conversion occurs and the very dense one where the mean-free path characterizes neutrino propagationhas been modified by recent theoretical studies. The inclusion of the neutrino-neutrino interaction in the treatment of the neutrino propagation, within the mean-field approximation has shown that significant flavor modification can occur in a region close to the neutrinosphere [14], with a possible impact on nucleosynthesis [61] and maybe on the explosion dynamics. Ref. [60] furnishes an example showing that a careful treatment of the transition region might be necessary, even from the point of view of the neutrino emission. Therefore, although the first investigations show that flavor conversion occurs out of the relevant region to influence the shock waves (see e.g. [62]), it might still be too early to draw definite conclusions.
In the context of core-collapse supernovae, a significant step forward beyond the mean-field approximation in the investigation of neutrino flavor conversion would require the solution of the Boltzmann equation implementing neutrino masses and mixings. Such numerical simulations are demanding, and appear still ahead also within realistic supernova simulations. The situation is different in the context of the early universe, since collisions are an essential ingredient bringing neutrino plasma to equilibration, while the mixings bring the system close to flavor equilibration. The neutrino history is usually determined by solving the Boltzmann equation for particles with mixings, or approximate versions of it.
Several works have addressed many-body aspects of the problem of neutrino flavor evolution focussing either on the theoretical formulation and the inherent symmetries [39], or on possible implications in schematic models (see e.g. [36]). In the present work we have adopted a novel perspective with respect to previous works and employed the BBGKY hierarchy as a theoretical framework to go beyond the mean-field approximation, in a realistic treatment of the neutrino interactions with matter and among themselves. In particular we have included contributions that have not been implemented so far. The BBGKY hierarchy formulation offers a natural truncation scheme with respect to the order of correlations. Such a formulation is equivalent to a description in terms of Green's function in the equal time limit. In particular, the BBGKY hierarchy furnishes a theoretical scheme to go from a many-body treatment to an effective one-body treatment. We have included, for the first time, novel νν correlations of the pairing type. We have obtained coupled non-linear evolution equations for the normal and abnormal densities, that depend upon the normal and pairing mean-fields. The explicit expression for such mean-fields has been obtained using the usual low energy limit for the charged-and neutral-current interactions.
The abnormal densities involve the expectation values of the bilinear products of the neutrino and antineutrino operators. Such pairs correspond to neutrinos with the same lepton number, while different individual lepton numbers can appear because of the presence of mixings, giving rise to off-diagonal contributions of the abnormal density. Further conditions are imposed to the neutrino-antineutrino pairs, in the calculation of the abnormal mean-fields, depending on the specific properties of the background that neutrinos are traversing. Requir-ing that the neutrino-antineutrino pairs have opposite momentum corresponds to a homogeneous medium. For this case we have derived a pairing mean-field with different helicity conditions for the pairs. In particular, if one considers contributions from states of positive (negative) helicity for the neutrinos (antineutrinos), the pairing mean-field turns out to depend on the ratio, of the neutrino mass over its energy, squared (as one expects). Although tiny, such off-diagonal contributions might give non-trivial effects. More generally one could consider the general case of an inhomogeneous background. The price to pay is that one should retain two indices in momentum in the evolution equations and all the normal and abnormal quantities involved.
An hypothesis made in the present work is that the abnormal densities are non-zero at initial time; while non-zero expectation values for these bilinear operator products can arise by implementing the collision term in a Boltzmann treatment. Such a collision term is not accounted for here. However we have been discussing that, in principle, following the BBGKY hierarchy one has the complete evolution equation for the two-body correlation function for our system of neutrinos and antineutrinos including both the collision term (also called the Born terms) and the terms dependent on the two-body correlation function (also called the PP terms) that have been assumed to depend upon the abnormal density only in the present work. The three-body correlation function contribution is expected to give higher order corrections. Therefore the framework discussed here could also be employed to investigate the relevance of the linked contributions of the two-body correlation function, that have made the object of the present work, for the case of cosmological neutrinos, at the epoch of Big-bang nucleosynthesis.
Obviously, the procedure we have been describing and employing for our calculations is very general. While it has been used to derive dynamical equations for a system of neutrinos propagating in a medium of ordinary matter and neutrinos, applying the procedure to more general cases is straightforward. For example, one can consider that matter is made up of other particles, or that it presents of non-standard interactions. Also, our results can be used for an arbitrary number of neutrino families and, in particular, to account for the presence of sterile neutrinos.
Finally, in the present manuscript, we have been focussing on the formal aspects implied by the two-body correlations of the neutrino-antineutrino type, and how to implement them in extended evolution equations.
Clearly numerical calculations are required to assess the impact of such correlations, or, more generally, of linked contributions of the two-body density matrix (or correlation function) on neutrino flavor conversion in a medium such as a core-collapse supernova, on nucleosynthesis and maybe the dynamics of these explosive phenomena, or on cosmological neutrinos at the epoch of primordial nucleosynthesis (within a Boltzmann treatment). Further investigations might tell us if these contributions can engender surprising features, or novel mechanisms, in this fascinating domain.
The above expression can be evaluated using the following definitions. We use the chiral representation. The four-component spinors can be expressed as | k| k 2 x + k 2 y k x k z − ik y | k| , k z k y + ik x | k| , −(k 2 x + k 2 y ) ,