Neutrino flavor oscillations in stochastic gravitational waves

We study neutrino flavor oscillations in a plane gravitational wave (GW) with circular polarization. For this purpose we use the solution of the Hamilton-Jacobi equation to get the contribution of GW to the effective Hamiltonian for the neutrino mass eigenstates. Then, considering stochastic GWs, we derive the equation for the density matrix for flavor neutrinos and analytically solve it in the two flavors approximation. The equation for the density matrix for the three neutrino flavors is also derived and solved numerically. In both cases of two and three neutrino flavors, we predict the ratios of fluxes of different flavors at a detector for cosmic neutrinos with relatively low energies owing to the interaction with such a GW background. The obtained results are compared with the recent observation of the flavor content of the astrophysical neutrino fluxes.


I. INTRODUCTION
After the recent success in experimental studies of solar [1], atmospheric [2], and accelerator neutrinos [3], these particles are believed to be massive and have nonzero mixing between different flavors. These neutrino properties result in transformations of the flavor content of a neutrino beam, called neutrino flavor oscillations [4].
Neutrino flavor oscillations can happen in vacuum. However, external fields, such as a magnetic field [5] or the electroweak interaction with background matter [6], can significantly modify the process of neutrino oscillations. The neutrino interaction with a gravitational field, in spite of its weakness, was reported in Refs. [7,8] to contribute the propagation and oscillations of flavor neutrinos. Recently, the method to account for the contribution of gravity, developed in Refs. [7,8], was further used in Refs. [9,10] to study neutrino flavor oscillations in static gravitational backgrounds such as Schwarzchild and Kerr metrics.
It is interesting to analyze how neutrino flavor oscillations proceed in a time dependent gravitational field, such as a gravitational wave (GW). This interest is inspired by the recent detection of GW emitted by merging binary black holes (BHs), reported in Ref. [11]. Later on, multiple GW signals, including that of the neutron stars coalescence, were observed. The summary of these events is presented in Ref. [12]. These phenomena are a strong evidence of the validity of the general relativity. Moreover, presently significant efforts are made to detect GW and astrophysical neutrinos, emitted by the same source (see, e.g., Ref. [13]). It would open a window for the multimessenger astronomy.
The influence of GW on neutrino oscillations was studied in Ref. [14], where spin oscillations are considered. In that situation, transitions between left and right polarized particles, belonging to the same neutrino type, were discussed. Neutrino spin oscillations driven by GW, studied in Ref. [14], are analyzed using the formalism for the description of the neutrino spin evolution in external fields in curved space-time developed in Refs. [15,16]. Note that the evolution of a spinning particle in GW was also considered in Ref. [17].
This paper, where we continue our study of neutrino oscillations in GW, is organized in the following way. First, in Sec. II, we adapt the formalism for the description of neutrino flavor oscillations, developed in Refs. [7][8][9][10], to describe the neutrino interaction with a time dependent gravitational field, such as GW. Using the solution of the Hamilton-Jacobi equation for a massive particle in GW, previously obtained in Ref. [18], we derive the effective Hamiltonian for neutrino flavor oscillations in GW. Then, considering a stochastic GW background, we obtain the equation for the neutrino density matrix and analytically solve it for the two neutrinos system. The equation for the density matrix of three flavor neutrinos is derived and solved numerically. In Sec. III, we consider a possible astrophysical application consisting in the interaction of cosmic neutrinos with random GWs emitted by coalescing supermassive BHs both in the two flavors approximation and in the general situation of the three neutrino flavors. We predict a specific flavor content of astrophysical neutrinos in a detector owing to the interaction with such GWs. In Sec. IV, we consider other random factors which can influence flavor oscillations of cosmic neutrinos on large distances. Our results are summarized in Sec. V. In Appendix A, we calculate the averaged phase of neutrino oscillations induced by stochastic GWs.

II. NEUTRINO FLAVOR OSCILLATIONS IN A GRAVITATIONAL WAVE
In this section, we study flavor oscillations of neutrinos in a plain GW. Then we assume that there is a stochastic background of GWs. The equation for the density matrix of mixed neutrinos, accounting for both the vacuum term and the interaction with GW, is derived and solved analytically. Then we also consider the general case of three flavor neutrinos and present the evolution equation for the density matrix in this situation.
We start with the discussion of the system of flavor neutrinos in the two flavors approximation, (ν e , ν x ), with x = µ, τ . These neutrino flavor eigenstates participate in the electroweak interaction with other leptons. However, these particles do not have definite masses. To diagonalize the mass matrix we introduce the neutrino mass eigenstates ψ a , a = 1, 2, which are related to ν λ by means of the matrix transformation ν λ = U λa ψ a , where U λa are the components of the mixing matrix. If we consider two flavor neutrinos, (U λa ) has the form, where θ is the vacuum mixing angle. The neutrino mass eigenstate with the mass m a was found in Refs. [7,8] to evolve in a gravitational field as where S a (x, t) is the action for this particle, which obeys the Hamilton-Jacobi equation [9,10], where g µν is the metric tensor. Instead of dealing with the neutrino wave functions in Eq. (2.2), we can define the contribution to the effective Hamiltonian H m for the mass eigenstates as where we take that neutrinos are ultrarelativistic particles. Equation (2.4) means that ψ a obeys the equation, iψ a = (H m ) aa ψ a . One can check that, in case of two neutrino eigenstates in vacuum, Eq. (2.4) results in the correct vacuum oscillations phase Φ vac = ∆m 2 /4E, where ∆m 2 = m 2 2 − m 2 1 > 0 and E is the mean neutrino energy. As a background gravitational field, we consider a plane GW with the circular polarization propagating along the z axis. Choosing the transverse-traceless gauge, we get that the metric has the form [19], where h is the dimensionless amplitude of the wave, φ = (ωt − kz) is the phase of the wave, ω is frequency of the wave, and k is the wave vector. In Eq. (2.5), we use the Cartesian coordinates x µ = (t, x, y, z).
The solution of Eq. (2.3) for a plane GW of the arbitrary form, not necessarily a monochromatic one as in Eq. (2.5), was found in Ref. [18]. If we define g ⊥ µν = g µν , at µ, ν = 1, 2, and this solution takes the form, where E = p 0 + p 3 and p µ ⊥ = (0, p 1 , p 2 , 0) are the integrals of motion of Eq. (2.3), x µ ⊥ = (0, x 1 , x 2 , 0), and v = t + z. One can see in Eq. (2.7), in case of a neutrino propagating along GW, i.e. when p µ ⊥ = 0, this background gravitational field does not affect the neutrino motion. If a neutrino interacts with a plane electromagnetic wave, there is an effect on neutrino oscillations in case when particles move along the wave [20]. The fact that GW cannot induce a spin flip of a neutrino propagating along the wave, i.e. it does not directly influence neutrino spin oscillations, was revealed in Ref. [14].
If we study the neutrino motion along GW and, then, adiabatically turn off the gravitational field, the action in Eq. (2.7) takes the form, S a = p 0 x 0 + p 3 x 3 , where p 0 = m 2 a + p 2 3 . It means that Eq. (2.7) has a correct vacuum limit.
Since typically h ≪ 1 in Eq. (2.5), we can decompose S a in Eq. (2.7) in a series, and keep only the terms linear in h. The symbol ". . . " in the arguments of S (i) a (. . . ) incorporates other parameters except the contribution of the gravitational interaction. Since the amplitude of GW has been explicitly written in Eq. (2.8), i.e. the effects of a nontrivial geometry have been taken into account, we can conclude that the additional parameters, marked by the "· · · " symbol, obey relations in a flat space-time. For example, S We suppose that a neutrino propagates arbitrarily with respect to GW. Hence the components of the neutrino momentum have the form, p 1 = p cos ϕ sin ϑ, p 2 = p sin ϕ sin ϑ, and p 3 = p cos ϑ, where ϕ and ϑ are the spherical angles. Using Eqs. (2.4)-(2.7), we get the diagonal entries of H m , which contain the linear contribution of GW, in the form, where E a = m 2 a + p 2 is the neutrino energy, φ a = ωt(1 − β a cos ϑ) is the phase of GW accounting for the fact that a neutrino moves on a certain trajectory, which is a straight line approximately, and β a = p/E a is the neutrino velocity. Note that (H  We have taken into account the contribution of GW linear in h to the diagonal elements of H m . However, besides GW, there are usual vacuum contributions to these elements, which have the form (H ) aa , as well as using Eqs. (2.1) and (2.9), we obtain that H f has the form, m ) aa is given in Eq. (2.9), and σ are the Pauli matrices.
We mentioned above that a neutrino, propagating along GW, is not affected by such GW. Therefore we consider the interaction of a neutrino with a stochastic GW background. In this situation, following Ref. [21], it is more convenient to deal with the density matrix ρ. Let us define We should average ρ I over the directions of the GW propagation and its amplitude. Then we consider the δ-correlated where τ is the correlation time. The evolution equation for ρ I , obtained in Ref. [21], has the form, where In Eq. (2.13), we averaged over the angles ϕ and ϑ and accounted for the fact that neutrinos are ultrarelativistic.
First, we should supply Eq. (2.12) with the initial condition. Taking into account that U 0 (0) = 1 in Eq. (2.11), we get that ρ I (0) = ρ(0). Then, we take that ρ I (0) 11 = F e is the initial probability (∼flux) of ν e , ρ I (0) 22 = F x = 1 − F e is the initial probability (∼flux) of ν x , and ρ I (0) 12 = ρ I (0) 21 = 0, that implies that there are no correlations between the initial fluxes of different flavors. Now, Eq. (2.12) can be solved analytically. The components of ρ I (t) have the form where is the parameter describing the relaxation of the density matrix. The expression for ρ (t) = U 0 ρ I (t)U † 0 is quite cumbersome in general case. We present its diagonal elements in the limit Γt ≫ 1, One can see in Eq. (2.16) that ρ 11 + ρ 22 = 1, as it should be. Of course, this relation holds true for arbitrary t.
Despite the two flavors approximation, adopted above, allows the analytical solution of Eq. (2.12), it cannot correspond to a realistic situation because the mixing between ν µ and ν τ is close to maximal [22]. That is why we should generalize our treatment of neutrino oscillations in stochastic GWs for three flavors.
We suppose that the condition ωL|β a − β b | ≪ 1, where a, b = 1, . . . , 3 and L is the neutrino propagation distance, established in Appendix A, is valid. Then, performing similar calculations, as above, and using the results of Ref. [21], we obtain the following equation for the averaged density matrix ρ I for the three neutrino flavors: Here ∆m 2 ab = m 2 a − m 2 b is the standard definition for the mass squared differences. Unlike Eq. (2.1), in three flavors case, the mixing matrix U in Eq. (2.18) can be parametrized in the form [4,22] where c ab = cos θ ab , s ab = sin θ ab , θ ab are the corresponding vacuum mixing angles, and δ CP is the CP violating phase. Note that Eq. (2.18) cannot be solved analytically. We present its numerical solutions in Sec. III when we discuss some possible astrophysical applications.

III. ASTROPHYSICAL APPLICATIONS
In this section, we apply the results of Sec. II for the description of neutrino flavor oscillations in GWs emitted by random merging binaries.
To study neutrino flavor oscillations in stochastic GWs, we should estimate h 2 and τ for this gravitational background. For this purpose, it is convenient to use the spectral function [23], where f = ω/2π, G = M −2 Pl is the Newton constant, M Pl = 1.2 × 10 19 GeV is the Planck mass, ρ c = 3H 2 0 /8πG = 0.53 × 10 −5 GeV · cm −3 is the critical energy density of the Universe, and S h (f ) is the spectral density. The root mean square of the strain h is then wheref is the typical frequency of stochastic GWs. Despite merging BHs with several solar masses or neutron stars are more abundant sources of GWs, we study coalescing supermassive BHs, with masses up to 10 10 M ⊙ . Such sources of stochastic GWs have lower characteristic frequenciesf . Thus, h 2 is greater for such sources. In this situation, we can take Ω GW (f ) = 10 −8 atf = 10 −6 Hz [24] to get the upper limit for h 2 . Using Eq. (3.3), we obtain that h 2 = 1.6 × 10 −32 . To evaluate the correlation time we take that τ ∼f −1 = 10 6 s.
Using the above estimates, we get the parameter, which describes the rate of the relaxation of the probabilities in Eqs. (2.14) and (2.15), while a relativistic neutrino passes the distance L = t, in the form
Let us suppose that the ratios of fluxes at a source are (F e : F µ ) S = (1 : 2) and (F e : F τ ) S = (1 : 0). This situation is consistent with the suggestion of Ref. [25], (F e : F µ : F τ ) S = (1 : 2 : 0). If κ νe→νx > 1 in Eqs. (3.5) and (3.6), then, using Eq. (2.16), we predict the ratios of the fluxes in a terrestrial detector in the form, where we take the mixing angles from Ref. [22]. Note that Eq. (3.7) is based on two flavors approximation ν e → ν µ and ν e → ν τ , i.e. we do not take into account ν µ ↔ ν τ transitions. Nevertheless, we can see that (F e : Before we proceed with the studies of the three flavors evolution, it is necessary to repeat the analytical result in Eq. (3.7) using the direct numerical solution of Eqs. (3.8) and (3.9), rewritten for two flavors. The reduction to the two flavors case is straightforward. This numerical solution is represented in Fig. 1. In Fig. 1(a), we show the fluxes of ν e and ν µ versus t ′ . This solution is based on the fluxes at a source (the initial condition) in the form, F e (0) = 1/3 and F µ (0) = 2/3. The characteristics of neutrinos and GW are taken as in Eq. (3.5). One can see in Fig. 1(b) that the asymptotic value of the fluxes ratio is (F e /F µ ) ⊕ = 0.9, which is in the agreement with Eq. (3.7). Performing analogous simulations, we show that (F e /F τ ) ⊕ = 21.8 in Eq. (3.7) can be also reproduced.
After reproducing the evolution of the two neutrino flavors with the numerical code, we can study the most general situation of three flavor neutrinos. The fluxes of ν e,µ,τ , based on the numerical solution of Eqs. (3.8) and (3.9), are shown in Fig. 2. To build Fig. 2 we utilize the values of ∆m 2 ab , θ ab , and δ CP from Ref. [22]. We consider both normal and inverted mass hierarchies in our simulations.
In Fig. 2, we take L = 1 Gpc and E = 0.5 MeV for the fluxes to reach their asymptotic values. This neutrino energy is between the values used in Eqs. (3.5) and (3.6). The propagation length, taken in Fig. 2, is comparable with the size of the visible Universe [26]. Figure 2 is based on the initial condition (at a source) (F e : F µ : F τ ) S = (1 : 2 : 0). We can see in Fig. 2(a) that, for the normal ordering, the asymptotic fluxes (at the Earth) are F e⊕ = 0.3127, F µ⊕ = 0.3504, and F τ ⊕ = 0.3369. For the inverted ordering, one has F e⊕ = 0.3154, F µ⊕ = 0.3497, and F τ ⊕ = 0.3349 in Fig. 2(b). It means that, at the Earth, the predicted fluxes are close to the case (F e : F µ : F τ ) ⊕ = (1 : 1 : 1). However, there is a small deviation from this prediction of Ref. [25] for both normal and inverted mass orderings. Moreover, one can see that there is a small dependence of our results on the hierarchy of the neutrino masses. Eventually, we can see that the deviation of the ratios of the cosmic neutrino fluxes, caused by the interaction with stochastic GWs, from the values (F e : F µ : F τ ) ⊕ = (1 : 1 : 1), predicted in Ref. [25], exists in the three flavors case. However, the magnitude of such a deviation is smaller than that in the two flavors approximation, studied above. The difference between the two and three flavors cases is explained by accounting for θ 23 , which is close to π/4. Hence, in the three flavors situation, ν µ ↔ ν τ oscillations are more intense and the total neutrino flux becomes more uniform.
It is important that we study the situation of the fixed distance between a neutrino source and a neutrino detector. The only random influence on neutrino oscillations is caused by the interaction with stochastic GWs. Some other random factors, which can influence neutrino flavor oscillations, are considered in Sec. IV.
The recent measurement of the flavor content of cosmic neutrinos in Ref. [27] excludes the following cases: (F e : F µ : F τ ) ⊕ = (1 : 0 : 0) and (F e : F µ : F τ ) ⊕ = (0 : 1 : 0). Our prediction of the neutrino fluxes at a source in Fig. 2 is in the region not excluded in Ref. [27]. Of course, neutrino energies in Ref. [27], E > 35 TeV, are much higher than these considered in our work, E = (10 2 keV ÷ several MeV); cf. Eqs. (3.5) and (3.6), as well as Figs. 1 and 2. Nevertheless there are prospects to detect cosmic neutrinos even with lower energies (see, e.g., Ref. [28]). Perhaps, the proposal for the observation of the diffused flux of cosmic neutrinos in Ref. [29] would make it possible to study the influence of GWs on neutrino oscillations.
At the end of this section, we discuss the approximation made to obtain Eqs. (2.14) and (2.15). In Appendix A, we find that the expressions for the neutrino fluxes in Eq. (2.16) are valid in the limit when λ = ωt|β 1 − β 2 | ≪ 1. If a neutrino is a relativistic particle, it means that t = L and Let us consider ν e → ν µ oscillations channel with ∆m 2 ⊙ = 7.4 × 10 −5 eV 2 . Taking the neutrino energy as in Eq. (3.5), E = 10 2 keV, and ω ∼ 10 −6 s −1 [24], we get that L crit = 6.5 × 10 2 Gpc. One can see in Eq. (3.5) that, if L ∼ 1 Gpc (the size of the Universe), then this L is much less than L crit . Therefore, the constraint in Eq. (3.10) is satisfied with a large margin. Analogously one can check that Eq. (3.10) is fulfilled for ν e → ν τ oscillations. Hence, the expression for the fluxes ν e and ν x , proportional to ρ 11 and ρ 22 in Eq. (2.16), are valid.

IV. OTHER RANDOM FACTORS CONTRIBUTING NEUTRINO FLAVOR OSCILLATIONS
In this section, we study other possible random factors which, along with stochastic GWs, can contribute flavor oscillations of cosmic neutrinos.
First, we mention that neutrino interaction with randomly distributed matter can affect flavor oscillations [21,30]. On the large distances L ∼ Gpc, studied in Sec. III, the contribution of this factor to H f in Eq. (2.10) does not exceed the following quantity: where G F = 1.17 GeV −2 is the Fermi constant, Ω B = 0.042 is the barions contribution to the total energy of the Universe, and m p is the proton mass. The contribution of GW to H f is V GW ∼ hΦ vac ∼ 10 −25 eV, where we take h = 10 −16 , E = 10 2 keV, and consider ν e → ν µ channel. We can see that V m ≪ V GW , i.e., the random matter contribution is negligible for flavor oscillations of such neutrinos. A terrestrial detector can record a neutrino flux emitted by randomly distributed sources. We should analyze this factor in studying of flavor oscillations of cosmic neutrinos and compare it with the contribution of stochastic GWs. For this purpose, we should study neutrino flavor oscillations in vacuum and average the probabilities over the distance of the neutrino beam propagation. We analyze this case in the two flavors approximation since it allows the analytical solution.
Using Eq. (2.10) and omitting H 1 there, we get that probabilities to observe ν e,x in a detector at the distance L from a source are P e,x (L) = cos 2 (Φ vac L) + cos 2 2θ sin 2 (Φ vac L) P e,x (0) + sin 2 2θ sin 2 (Φ vac L)P x,e (0), (4.2) where P e,x (0) are the emission probabilities. To obtain Eq. (4.2) we assume that ν e (0)|ν x (0) = 0. Now, we average Eq. (4.2) over L by setting cos 2 (Φ vac L) = sin 2 (Φ vac L) = 1/2. Finally, we get the corresponding probabilities which formally coincide with these in Eq. (2.16). This fact is not surprising. It is a consequence of the ergodic theorem [31]. Nevertheless, Eqs. (2.16) and (4.3) correspond to completely different physical situations. In Sec. II, we study the neutrino propagation between a detector and a source which are at the fixed points. The random influence on neutrino oscillations by stochastic GWs is between these points. To derive Eq. (4.3), we consider flavor oscillations in vacuum of neutrinos emitted by randomly distributed sources. No other stochastic influence on neutrino system is assumed now.
To differentiate between these cases, we can consider a neutrino source with adiabatically changing luminosities of different neutrino flavors. In this situation, the neutrino fluxes at the source in Eq. (2.16) are slowly varying functions of time F e,x (t). On the contrary, one cannot expect that P e,x (0) change simultaneously in all sources, especially when these sources are causally disconnected.

V. CONCLUSION
In this work, we have studied neutrino flavor oscillations under the influence of a plain GW with the circular polarization for the first time. In Sec. II, we have analyzed the evolution of the mass eigenstates in the quasiclassical approximation. Using the expression for the action for a massive particle, interacting with GW, obtained in Ref. [18], we have derived the contribution of GW to the effective Hamiltonian for the neutrino mass eigenstates. We have revealed that, in case of the neutrino propagation along GW, GW does not influence neutrino flavor oscillations.
Then, we have assumed that we deal with stochastic GWs emitted by randomly distributed sources. In this situation, we have derived the equation for the density matrix of neutrinos using the approach in Ref. [21]. This equation has been solved analytically in the two flavors system. The asymptotic expressions for the diagonal elements of the density matrix, which the fluxes of flavor neutrinos are proportional to, have been presented in Eq. (2.16). The equation for the density matrix evolution for the three flavor neutrinos has been also derived in Sec. II; cf. Eqs. (2.17)- (2.19). However, this equation can be solved only numerically.
Then, in Sec. III, we have considered an astrophysical application of the obtained result. We have supposed that GWs are emitted by merging binary supermassive BHs. Using the two flavors approximation, we have obtained that the probabilities to detect relatively low energy neutrinos, with E = (10 2 keV ÷ several MeV), can reach the asymptotic values in Eq. (2.16) if the propagation distance is comparable with the size of the Universe L ∼ Gpc [26]. Recently, we revealed in Ref. [14] that spin oscillations can be significantly affected by GW only if the neutrino energy is low.
Note that the predicted fluxes of different flavors are not equal at a detector, contrary to the finding of Ref. [25]. In the two flavors approximation, they depend on the fluxes at a source and the vacuum mixing angle. The fluxes of flavor neutrinos at the detector in Eq. (3.7) are in a region not excluded by the observations in Ref. [27].
In Sec. III, we have also studied neutrino flavor oscillations in stochastic GWs, emitted by merging BHs, in the most general case of the three flavors. The deviation of the fluxes of flavor neutrinos at a detector from the prediction of Ref. [25], (F e : F µ : F τ ) ⊕ = (1 : 1 : 1), owing to the neutrino interaction with stochastic GWs, has been confirmed in the three neutrinos situation by means of the numerical simulations; see Fig. 2. This deviation turns out to be smaller than that in the two flavors approximation.
The prediction of the fluxes at a detector in Sec. III is based on the assumption that the distance between a source and a detector is fixed. The stochastic influence of external fields (GWs in our case) is applied between the emission and detection points. We have found in Sec. IV that the same asymptotic fluxes are obtained if one averages vacuum probabilities over the neutrino propagation distances. It is the consequence of the ergodic theorem [31]. In Sec. IV, we have have also pointed out how to pick out these completely different physical situations.
Finally, we mention that the results obtained in the present work have some uncertainty since no stochastic GW background from supermassive BHs has been observed yet. It is mainly related to the determination of the correlation time τ . The fact that GWs, emitted by merging BHs with several solar masses, have been observed [12], imposes strong constraints on the parameters of such GWs. Nevertheless there are efforts to detect stochastic GWs with 10 −9 Hz < f < 10 −6 Hz (see, e.g., Ref. [32]).

ACKNOWLEDGMENTS
I am thankful to G. Sigl and S. V. Troitsky for useful comments. This work is performed within the government assignment of IZMIRAN. I am also thankful to RFBR (Grant No. 18-02-00149a) and DAAD for a partial support.