Neutrinophilic Axion-Like Dark Matter

The axion-like particles (ALPs) are very good candidates of the cosmological dark matter, which can exist in many extensions of the standard model (SM). The mass range of the ALPs as the dark matter can extend from a sub-eV scale to almost $10^{-22}~{\rm eV}$. On the other hand, the neutrinos are found to be massive and the SM must be extended to explain the sub-eV neutrino masses. It becomes very interesting to consider an exclusive coupling between these two low scale frontiers that are both beyond the SM. The propagation of neutrinos inside the Milky Way would undergo the coherent forward scattering effect with the ALP background, and the neutrino oscillation behavior can be modified by the ALP-induced potential. Assuming a derivative coupling between the ALP and the three generations of active neutrinos, possible impacts on the neutrino oscillation experiments have been explored in this paper. In particular, we have numerically studied the sensitivity of the Deep Underground Neutrino Experiment (DUNE). The astrophysical consequences of such coupling have also been investigated systematically.


I. INTRODUCTION
As one of the most promising dark matter candidates, the ALP with a possible mass spanning from sub-eV to O(10 −22 ) eV 1 is drawing more and more attention (see [5][6][7][8] for recent reviews). The typical QCD axion [9,10] was first identified as the pseudo-Nambu-Goldstone (pNG) boson in the Peccei-Quinn (PQ) mechanism [11,12] to solve the strong CP problem of QCD. It was later realized [13][14][15][16][17] to be a very good candidate of the cold dark matter which can acquire an effective anomalous mass by interacting with the gluons. The mass of the QCD axion is settled by the QCD phase transition scale Λ QCD and the PQ scale f a , e.g. m a ≈ Λ 2 QCD /f a ≈ 10 −5 eV (10 12 GeV/f a ) with Λ QCD ≈ O(10 2 ) MeV and f a ≈ 10 9 -10 12 GeV. Moreover, the axion can interact with the standard model (SM) fermions through a derivative type of coupling whose strength is proportional to its mass m a . On the other hand, the ALPs with very similar properties to the QCD axion are well motivated in many other extensions of the SM. For example, generating the tiny neutrino mass by spontaneously breaking the lepton number [18][19][20][21][22][23][24] predicts the existence of a Nambu-Goldstone (NG) boson, the Majoron, which inevitably couples with the neutrinos. To spontaneously break the family symmetries will lead to the familons [25,26]. In the string theory framework, many different ALPs can appear naturally [27][28][29][30][31][32], and one of them could just be the QCD axion. In some unified models [33][34][35][36][37][38][39], the pseudoscalar particle can even * huanggy@ihep.ac.cn † newton@ihep.ac.cn 1 The ultralight dark matter with a mass ∼ 10 −22 eV is often called "fuzzy dark matter" [1]. Its macroscopic wavelength can suppress the power of the galactic clustering and resolve the small scale tension of the cold dark matter paradigm [2]. A lower bound on the ALP mass can be set by the observation of the Lyman-α forest at O(10 −22 ) eV [3,4].
play multiple roles among the QCD axion, the majoron and the familon at the same time.
The mass spectrum and the coupling strength with the SM particles of the generic ALPs are not limited as in the QCD axion model, greatly enriching their phenomenology and the experimental efforts to hunt them. The ALPs can be non-thermally produced in the early Universe by the misalignment mechanism or the topological defect decay. After the ALPs form the dark matter halo of our Milky Way, their de Broglie wavelength λ dB can be as long as several kpc for an ALP mass m a ≈ 10 −22 eV, almost comparable to the scale of the Milky Way ∼ O(10) kpc. The local number density of the ultralight particles is n a = ρ a /m a ≈ 10 30 cm −3 if they saturate the dark matter energy density ρ ≈ 0.3 GeV · cm −3 . The huge particle number within the de Broglie wavelength makes them oscillate coherently as a single classical field with a(t, x) = a 0 cos (m a t − p · x), where a 0 = √ 2ρ/m a is the amplitude of the ALP field, p is the momentum of the ALP, and m a denotes the mass of the ALP. The ALP field maintains its coherence during the expansion of the universe. However, after the formation of the Milky Way, it develops a slight stochastic decoherence due to the gravitational clustering effect. The typical velocity dispersion of the ALP dark matter is σ v ∼ 10 −3 c [40], corresponding to a coherence length of λ coh ≈ λ a × 10 3 with λ a = 2π/m a being its Compton wavelength. In our work, we consider an ALP-neutrino interacting Lagrangian which preserves the shift symmetry as −L int = g αβ ∂ µ a ν α γ µ γ 5 ν β , where g αβ with α, β = (e, µ, τ ) is the dimensional coupling strength in the neutrino flavor basis with g αβ = g * βα , and ν α(β) represents the neutrino flavor eigenstate. In the QCD axion model, the coupling strength g αβ is suppressed by the PQ scale f a , which is not necessarily the case in the generic ALPs context.

arXiv:1809.01111v1 [hep-ph] 4 Sep 2018
The neutrino propagating in the Milky Way can coherently interact with the neutrinophilic ALP field, and the mixings among the neutrino flavors will be modified effectively. In this note, we will study systematically the impacts of the ALP-neutrino coupling on the neutrino oscillation experiments as well as the astrophysical phenomenology. In particular, we will take DUNE as a typical example, and numerically study its sensitivity. The influences of the possible ultralight dark matter field on the neutrino phenomenology have been discussed from various perspectives in the literature [41][42][43][44][45][46][47]. However, it is still very worthwhile to conduct a work like this one. In the existing works, the ultralight scalar [41][42][43][44][45][46][47], vector [42,45] and tensor [45] DMs have been considered for the neutrino oscillation experiments. The pseudoscalar ALP with a derivative coupling is now being investigated in this work. We will see the derivative coupling in Eq. (2) can have very different laboratory and astrophysical consequences. The impacts on the neutrino oscillation experiments greatly depend on the duration of one ALP oscillation cycle. The time-dependent perturbation analysis should be adopted when the ultralight DM oscillates a number of cycles during the neutrino propagation. If the ALP field loses its coherence within a single neutrino flight, the forward scattering effect which is coherently enhanced should vanish stochastically. Various astrophysical constraints on the coupling are considered systematically in this work. The free-streaming constraint of the cosmic microwave background (CMB) for the neutrino decay process should be the most stringent one. The influences on the propagation of supernova neutrinos and ultrahigh-energy (UHE) neutrinos are also discussed in detail.
This paper is organized as follows. In Section II, the influence of the ALP dark matter on the neutrino oscillation behavior have been explored in detail, and the impact on DUNE has been elaborated with a numerical sensitivity study. In Section III, we will generally discuss the existing constraints from astrophysical observations. We make our conclusion in Section IV.

II. IMPACTS ON NEUTRINO OSCILLATIONS
When neutrinos propagate inside our galaxy, which is always true for the current ground-based oscillation experiment, they will inevitably scatter with the ambient ALP background under our consideration. The coherent forward scattering process which can be analogous to the Mikhyev-Smirnov-Wofenstein (MSW) [49][50][51] matter effect is dominant. Under the interaction form in Eq. (2), the dispersion relation of the three generations of neutrinos will be modified as whereÎ is the 3×3 identity matrix, g m stands for the 3×3 Hermitian matrix of the coupling constant g αβ , m ν stands for the mass matrix of neutrinos in the flavor basis, E ν and p ν are the energy and momentum of the neutrino, respectively. The modified neutrino oscillation behavior can be conveniently described by the following effective Hamiltonian: where U stands for the flavor mixing matrix of the three generations of neutrinos in vacuum, ∆m 2 21 and ∆m 2 31 are the neutrino mass-squared differences in vacuum. Here V is a potential characterizing the MSW effect with the possible ordinary matter around the neutrino, and ξ αβ for α, β = (e, µ, τ ) is the potential contributed by the ALP background. ξ αβ can be derived from Eq. (3), which at the leading-order reads where p ν /| p ν | stands for the direction of the neutrino flux. Providing the velocity of the ALP field in the Milky Way is of the order O(10 −3 )c, the temporal term should dominate the potential in Eq. (5). The further approximation can be made to give with √ 2ρ ≈ 2.15 × 10 −3 eV 2 ρ/(0.3 GeV · cm −3 ). As has been mentioned before, the ALP potential oscillates with a period t a ≈ 2π/m a ≈ 4.8 day (10 −20 eV/m a ). The impact on the neutrino oscillation strongly depends on the magnitude relation between the ALP period and several time scales of the experiment. The relevant experimental time scales include: t opr , the operation time of the experiment, typically over several years for an oscillation experiment; t rsv , the minimal period that the experiment can resolve for a periodic modulation effect; t flt , the flight time of each single neutrino from the source to the detector, ∼ 0.33 ms for a baseline length of 100 km. There are four different cases that should be addressed: (i) Case I: t opr t a . In this case, neutrinos will experience an approximately constant ALP-induced potential throughout the operation time of the experiment. The influence to neutrino oscillations is very similar to the non-standard interactions (NSIs) of neutrinos 2 with a constant matter density profile. But the effect of the ALP field is irreducible for all oscillation experiments in our galaxy, even with no ordinary matter surrounding the neutrino. For an ALP mass m a ≈ 10 −22 eV which might be the minimal feasible mass to form the dark matter, the corresponding period is t a ≈ 1.3 yr. However, the operation time of most oscillation experiments is longer than 1.3 yr. For them, the ALP-induced potential is no longer constant during the operation. (ii) Case II: t rsv t a t opr . The operation time is longer than the period of the ALP field, so it can cover a number of ALP cycles during the operation. To detect the modulation effect of the oscillating potential, the duration of each ALP cycle must be longer than the experimental resolution limit of the periodicity. There are several experiments that have looked for the periodic modulation effect of the solar neutrinos, e.g. Super-Kamiokande [52][53][54][55] and SNO [55][56][57][58]. The period of possible periodic signals have been scanned from ∼ 10 min to ∼ 10 yr. No anomalous modulation effect is found in the data sets of Super-Kamiokande and SNO using various statistical methods. Based on this, Ref. [41] has set a strong constraint on the ultralight scalar coupling. The minimal recognizable period t rsv of an experiment depends on many factors: the time binning of the data, the events number, the statistical approach etc.. In principle, t rsv can be as small as the precision of the time measurement, e.g. ∼ 100 ns for SNO [57]. The statistical analysis in the solar case can be similarly applied to other types of oscillation experiments, to locate the ALP-induced signal. Moreover, in this case one can also choose to integrate over the data-taking time to get a smaller averaged effect [42][43][44]. (iii) Case III: t flt t a t rsv . The ALP field oscillates so fast that the experiment can no longer resolve the periodic modulation effect. But during each single neutrino flight from source to detector, the ALP field is still approximately constant. There is no other way but to integrate over the data-taking time and obtain the averaged distortion effect. (iv) Case IV: t a t flt t a . Throughout a single neutrino flight, the ALP field has been oscillating for a number of cycles. In this case, the neutrino flavor evolution is driven by a varying potential. The time-dependent perturbation theory is required for the perturbative expansion analysis. However, if the neutrino flight time is even longer than the decoherence time of the ALP t a ≈ 1000 × t a , i.e. t flt > t a , the coherent potential induced by the ALP field should vanish stochastically.

A. Perturbative Expansion
We will work under the two-neutrino-flavor scheme to see the influence of the ALP potential analytically. Let us assume the ALP-induced potential is very small, such that one can perturbatively expand the effective mixing parameters and the oscillation probability around the standard ones. Taking ν e and ν µ as the two neutrino fla-vors under our consideration, the effective Hamiltonian in Eq. (4) is then reduced to 3 where ∆m 2 (t) is the effective mass-squared difference which has been shifted by the ALP background, andŨ (t) is the effective mixing matrix that diagonalizes the total Hamiltonian. Setting θ andθ as the mixing angles of the matrices U andŨ (t) respectively, one can formally expand the effective quantities to the second-order as Here δ 1 and δ 2 represent the first-order and the secondorder ALP perturbations respectively, which read as follows with ξ ≡ ξ µµ − ξ ee , and ξ

R(I) eµ
being the real (imaginary) part of ξ eµ . One can observe that ∆m 2 and θ have been shifted by amounts of ξE ν and ξE ν /∆m 2 respectively, with ξ denoting the general order of magnitude of ξ αβ . Assuming the ALP potential is constant within t flt (for Case I, Case II or Case III), the survival probability of ν µ → ν e reads where φ ≡ ∆m 2 L/(4E) is a dimensionless phase factor and the corrections are +2 sin 2 φ sin 4θ × δ 1 (θ), When the ALP field oscillates fast enough like in Case II and Case III, the first-order term δ 1 (P µe ) ∝ sin m a t can be averaged to 0. However, the second-order term after being averaged, δ 2 (P µe ) ∝ sin 2 m a t = 1/2, will lead to the distortion effect on the final probability. It should be emphasized that the second-order corrections in Eq. (8) are as important as the first-order ones when we derive δ 2 (P µe ) in Eq. (15). We notice that δ 1 (P µe )/P µe ≈ ξE ν /∆m 2 and δ 2 ( The above discussion is under the assumption that the ALP field is approximately constant during a single neutrino flight. For Case IV, the above results do not hold any more. In the case of the time-dependent potential, the transition amplitude M ≡ ν e |ν µ t can be formally expanded as M = M 0 + M 1 + M 2 + · · · with M i being the ith-order result. Using the time-dependent perturbation theory, one can work out the following first-order amplitude correction, where the subscripts i, j run over (1,2). Here t 1 and t 2 denote the initial time and final time of the neutrino flight, respectively,Ĥ 0,11 = 0 andĤ 0,22 = ∆m 2 /(2E ν ) are the eigenvalues for the Hamiltonian in vacuum H 0 , and we define with η, γ run over (e, µ). One can notice a singularity at the point m a = ∆m 2 /(2E ν ) in the denominator. However, the singularity can be cancelled by the numerator and reduce to a factor ∼ t flt . The baseline of an oscillation experiment is usually selected around the oscillation maximum, i.e. t flt ≡ (t 2 −t 1 ) ≈ 2πE ν /∆m 2 . For Case IV with t a ≡ 2π/m a < t flt , we should have ∆m 2 /E ν m a . One might as well further set ∆m 2 /(2E ν m a ) O(0.1), such that Eq. (16) can be approximated as The final survival probability of ν µ → ν e is given by As has been mentioned before, the flight time for a baseline length of L = 100 km is t flt ≈ 0.3 ms. For Case IV under the discussion, the oscillating period of ALP should be in the range of (0.3 µs, 0.3 ms). If the experiment can resolve the periodic fluctuation in this time range (that would require a very refined scanning in the experimental data set), the modulation pattern with an amplitude ξ/m a can be identified. Regardless of the practical periodicity resolving power, one can always choose to average over the observation time, obtaining the probability correction at the order of O(ξ/m a ) 2 . Note that the cross term Re(M o M * 1 ) vanishes for the averaged probability. In conclusion, the ALP field will induce a potential ξ ≈ g √ 2ρ sin m a t in the neutrino effective Hamiltonian, with g denoting the general order of magnitude of the coupling strength g αβ . If the ALP field is approximately constant for a single neutrino flight and its periodicity is also resolvable for the concerned experiment, the potential will shift ∆m 2 and θ by E ν ξ and ξE ν /∆m 2 respectively. No matter whether the ALP periodicity is resolvable for the experiment, one can always average over the observation time. In this case, the shifts of ∆m 2 and θ are (E ν ξ) 2 and (ξE ν /∆m 2 ) 2 respectively. If the ALP field oscillates very rapidly along the neutrino course but still remains in coherence, the oscillation probability will also be modified. The probability can either fluctuate with the amplitude ξ/m a or averagely get distorted by the magnitude (ξ/m a ) 2 .

B. General Sensitivities
In this section, we will roughly estimate the sensitivities (or constraints) of the oscillation experiments based on the analytical results of the last section. Suppose that an oscillation experiment can measure ∆m 2 and θ with the 1σ experimental errors of σ(∆m 2 ) and σ(θ). The experiment should be sensitive to the ALP coupling strength whose corrections to ∆m 2 and θ are bigger than the experimental errors. To obtain the rough sensitivity, one can demand that the ALP-induced corrections should be smaller than errors of the relevant ∆m 2 and θ. This argument should be good enough as an orderof-magnitude estimate which is the main purpose of this section. Setting σ ≡ min σ(∆m 2 )/∆m 2 , σ(θ) , the argument can be translated into the following conclusive results for two different scenarios: when the ALP potential is either approximately constant or periodically resolvable; and when the oscillating potential is averaged to the distortion effect. We have used a factor max (2π, m a t flt ) and the condition t flt ≈ 2πE ν /∆m 2 to combine Cases I, II, III and Case IV, which is good enough as an order-ofmagnitude estimate. With higher energy the oscillation experiment can have stronger sensitivity. One can observe that the sensitivity (constraint) can be made by an experiment depends on: (1) its baseline; (2) its beam energy; (3) the resolution power of the periodicity; (4) the statistics. We list in Table 1 [63], which correspond to σ(θ 12 ) ≈ 0.18% and σ(∆m 2 21 )/∆m 2 21 ≈ 0.24%. With an exposure of 150 kt · MW · yr for DUNE, the resolutions of sin 2 2θ 13 , sin 2 θ 23 and ∆m 2 31 can reach 0.007, 0.012 and 1.6 × 10 −5 eV 2 [64], respectively, corresponding to σ(θ 13 ) ≈ 0.63%, σ(θ 23 ) ≈ 1.2% and σ(∆m 2 31 )/∆m 2 31 ≈ 0.64%. For the atmospheric neutrino experiment, we will take its baseline and energy as 7 × 10 3 km and 10 GeV as an example. It is able to push the 1σ errors of sin 2 θ 23 and ∆m 2 31 down to 0.05 and 0.1 × 10 −3 eV 2 respectively [65,67], corresponding to σ(θ 23 ) ≈ 5% and σ(∆m 2 31 )/∆m 2 31 ≈ 4%. In Fig. 1, the sensitivities or constraints of different oscillation experiments are demonstrated. For comparison, the existing astrophysical bounds which will be discussed in Section III have also been given. In the left panel, the sensitivities and constraints are given based on Eq. (19) for the averaged case, while the right panel is based on Eq. (18) for the non-averaged case. The red, blue and green curves demonstrate the future sensitivity of DUNE, JUNO, and the atmospheric neutrino experiment respectively. They can exceed the most stringent astrophysical bound given by the CMB free-streaming argument. One should keep in mind that the sensitivity study here is based on the rough arguments, while the exact one can only be done with the dedicated experimental simulation. We will take DUNE as an example with detailed numerical simulations and analysis in the next section. We will find the rough estimation in this section is consistent in the order of magnitude with the result given by the detailed simulation.

C. Impact on DUNE
In this section, we will perform a detailed study of the ALP-neutrino coupling parameters considering DUNE both at the probability as well as the χ 2 level. DUNE is a proposed long baseline superbeam neutrino oscillation experiment at Fermilab [64,68]. It will have two neutrino detectors. The near detector will be placed near the source at Fermilab, while the far one will be installed 1300 km away from the neutrino source at Sanford Underground Research Facility (SURF) in Lead, South Dakota. The far detector will utilize four 10 kton liquid argon time projection chambers (LArTPCs). Also, it will use the existing Neutrinos at the Main Injector (NuMI) beamline design at Fermilab as the neutrino source.
For the numerical analysis of DUNE, we use the GLoBES package [69,70] along with the required auxiliary files given in Ref. [68]. Throughout this work we perform our simulation considering the 40 kton liquid argon detector. The flux corresponding to an 1.07 megawatt beam power gives 1.47×10 21 protons on target (POT) per year due to an 80 GeV proton beam energy has been considered. The remaining experimental details like the signal and background normalization uncertainties for both the appearance and disappearance channels have been taken from DUNE CDR [68]. In addition, we consider a 3.5 years running time for each of the neutrino and antineutrino modes unless otherwise stated.
To study the impact of the ALP-neutrino coupling on DUNE, we use the GLoBES extension file snu.c as has been presented in Refs [71,72]. In our analysis, we have modified the NSI matter potential in snu.c with the potential induced by the ALP field. We simulate the fake data of DUNE using the standard oscillation parameters sin 2 θ 12 = 0.307, sin 2 θ 13 = 0.022, sin 2 θ 23 = 0.535, δ = −π/2, ∆m 2 21 = 7.40 × 10 −5 eV 2 and ∆m 2 31 = 2.50 × 10 −3 eV 2 which are compatible with the latest global-fit results [73][74][75]. We marginalize the standard parameters θ 23 , δ and the mass hierarchy over their current 3σ ranges in the test. As the remaining standard oscillation parameters have been measured with very high precisions, we keep them fixed in the analysis. On the other hand, for the ALP-neutrino coupling parameters, we fix their true values as zero during the statistical analysis. For the fit we switch on only one of them at one time for numerical simplicity. While performing the sensitivity study of the diagonal ALP-neutrino coupling parameters, we marginalize over the standard parameters in the fit. Moreover, for the off-diagonal parameters, we also marginalize their phase in the range (0 → 2π). All along this work, we perform our analysis considering two different scenarios for the ALP-induced potential, namely non-averaged and averaged cases. For the non-averaged case, we take the ALP-induced potential ξ αβ as a constant over time for simplicity. A more careful treatment would require the statistical analysis of the modulation on the time-binned data, which might be interesting for a future work. For the averaged case, we have modified the snu.c file by averaging the oscillation probability over the time points within one ALP cycle.
In Fig. 2, we show the oscillation probabilities of the appearance channels (ν µ −→ ν e and ν µ −→ ν e ) for DUNE. For illustration, we perform this study with a single non-zero representative value (g ee = 5.6×10 −11 eV −1 ) of the ALP-neutrino coupling term. The green (dashed and dotted) curves show the case where the ALP-induced potential is not averaged over the observation time (see the figure legend for more details). Whereas, the solid blue curves describe the scenario where the ALP-induced potential is averaged over time. Besides this, the stan-

Parameters
Non-averaged Averaged  dard case without the ALP-induced potential is shown as the gray solid curves for comparison. For the nonaveraged case, one can observe significant deviations of the probabilities around the peak energy ∼ 2.5 GeV compared to the standard case. When the ALP-induced potential oscillates as g ee √ 2ρ sin m a t, the associated probability will vary roughly in between the two green curves. After averaging over the time, one would obtain the distortion effect as the blue curve in each panel, which is smaller than the non-averaged case. As our intention at the probability level is for demonstration, hence we illustrate this analysis with a single non-zero ALP-neutrino coupling term for the appearance channel. Nevertheless, we will perform a detailed sensitivity study considering both the appearance and disappearance channels for all the coupling terms at the χ 2 level.
In Fig. 3, we show the sensitivity of DUNE to the six ALP-neutrino coupling parameters g αβ in the χ 2g αβ plane. The top and bottom panels represent our results for the non-diagonal and the diagonal parameters, respectively. The green solid curves represent the non-averaged case where a constant potential has been taken during the numerical simulation. The blue dotted curves signify the averaged case where the ALP-induced potential is averaged over time in the final probability. Note that the black dotted and black dashed horizontal lines stand for the χ 2 values corresponding to the 2σ and the 3σ C.L., respectively. In Fig. 3, one can clearly observe that the averaged case has looser sensitivity than the non-averaged case. This can be easily understood from the analytical results before and from the probability demonstration in Fig. 2. It can be ascribed as the loss of information after averaging. From the first plot of the bottom row, we notice that for the green curve there is a local minimum around g ee = −10.1 × 10 −11 eV −1 other than a global minimum at zero. This shows a degenerate solution for g ee and it arises from the unknown hierarchy (as we have marginalized over the mass hierarchy in the fit). In our careful analysis by fixing the hierarchy to the normal one both in data and theory we find no degenerate solution for g ee . This tells that the lack of knowledge of the hierarchy may lead to a degenerate solution. Thus, for all these different cases we perform our numerical analysis over the marginalized hierarchy. Also, comparing the top and the bottom row, we notice that overall sensitivity of the off-diagonal parameters is better than that of the diagonal ones. Table II summarizes the expected sensitivities for all the ALP-neutrino coupling parameters g αβ for DUNE, where the 2σ and 3σ intervals are presented respectively. Note that the second and third panels of the table show the sensitivity limit for the non-averaged and averaged cases, respectively.

A. Early Universe
The cosmological evolution will be more or less modified, if neutrinos have strong interactions with an ultralight degree of freedom. During the expansion of the Universe, the thermal ALPs can be produced via the processes like ν + ν ↔ a + a and ν i ↔ ν j + a 4 , where ν i and ν j are the neutrino mass eigenstates with m i > m j . If the ALPs are fully thermalized before the neutrino decoupling around 1 MeV, they will contribute to the extra effective neutrino number by an amount of ∆N eff ≈ 0.57 [76]. The extra N eff is strongly constrained by the observations of BBN and CMB. The current constraint of BBN is ∆N eff 1 at 95% C.L. [77], while CMB can give a stronger constraint with the latest Planck 2018 result [78]: ∆N eff 0.52 (95% C.L., P lanck TT+lowE). On the other hand, the strong interaction during the CMB epoch will also suppress the free-streaming length of neutrinos [79,80], such that the evolution of CMB perturbations is severely disturbed. This can lead us to a very strong constraint on the coupling strength of the secret neutrino interactions between neutrinos and ALPs.
First, let us focus on the possible influence on ∆N eff . We will find the binary process ν + ν ↔ a + a is the dominant one for generating ALPs. The reaction rate of the binary process can be approximated as Γ ≈ g 4 m 2 ν ·T 3 with T being the temperature of the neutrino plasma 5 . During the epoch of radiation domination, the Hubble expansion rate is given by H ≈ 1. 66 √ g * T 2 /M Pl , with g * denoting the effective number of relativistic degrees of freedom at T (refer to [82,83] for its values at various temperatures) and M Pl 1.221 × 10 19 GeV being the Planck mass. Because of the ratio Γ/H ∝ T , the thermal ALPs should first be copiously produced at high temperatures in the very beginning of the Universe. They should eventually freeze out at some low temperature T fo to evade the CMB constraint. The extra N eff contributed by the thermal ALPs is found to be for T fo > T ν , where T ν is the neutrino decoupling temperature around 1 MeV, and g * s ≈ g * holds before T ν . To satisfy ∆N eff 0.52, one can find the freeze-out temperature of thermal ALPs should be T fo 20 MeV [84]. 4 We assume the neutrino mass is much larger than the ultralight ALP such that the process like ν + ν ↔ a is kinematically forbidden. 5 The derivative coupling in Eq. (2) is equivalent to the pseudoscalar coupling h αβ aν α γ 5 ν β with a dimensionless coupling constant h αβ = ij U αi U * βj (m i + m j )g αβ , only if each neutrino line in the Feynman diagram is attached to only one NG boson line [79,81]. However, this is not the case for the process ν + ν ↔ a + a. The reaction rate for the pseudoscalar coupling case is proportional to T instead of T 3 .

Demanding Γ
H at 20 MeV, we obtain the following constraint The ∆N eff generated by the freeze-in process ν i → ν j + a yields a negligible bound g 10 −2 eV −1 , therefore not elaborated here.
During the CMB era, the process ν + ν ↔ a + a is suppressed by the low plasma temperature. There are mainly two processes that can reduce the neutrino freestreaming length: ν+ν ↔ ν+ν and ν i ↔ ν j +a. For these two processes, the derivative coupling is equivalent to the pseudoscalar coupling. By requiring Γ(ν + ν ↔ ν + ν) H at the photon decoupling temperature T γ ≈ 0.256 eV [34,79], one can obtain the constraint The decay process ν i ↔ ν j + a is only relevant for the off-diagonal coupling g ij with i = j in the mass basis. The decay rate differs for different choices of neutrino mass patterns. If the three masses of neutrinos are very hierarchical, the decay rate reads Γ ij ≈ g 2 ij m 4 i /(48πT ) with g ij ≡ ij U αi U * βj g αβ and m i being the mass of the heavier neutrino ν i . To make sure the free propagation of neutrinos is sufficiently disrupted, one should include an angular factor (m i /3T ) 2 [79]. Requiring Γ t ≡ Γ ij (m i /3T ) 2 H at T γ , the following bound can be derived [79]: However, if the neutrino masses are degenerate, one should instead have which is in the same order of magnitude as the the former one in Eq. (23). Overall, the most stringent constraint can be placed on the off-diagonal couplings in the mass basis g ij 10 −10 eV −1 based on the free-streaming of CMB. The ∆N eff limit of CMB can place bounds on all elements of g αβ as g 10 −8 eV −1 . These early Universe bounds have been included in Fig. 1 as the dash-dotted gray curves.

B. SN1987A Explosion
The ALPs can also play an important role in the evolution of the supernovae. The constraints on the ultralight scalar coupling with neutrinos from the supernova neutrino observation SN1987A [85,86] have been extensively studied in the literature, see Ref. [80,81] and the references therein for more details. The presence of the ALP-neutrino coupling can modify the standard evolution of supernova explosion by the following processes: (1) The ALPs can be thermally generated inside the core, carrying energy away from the core-collapse supernovae; (2) The ALPs can transport energy inside the supernova core and reduce the cooling time or the duration of the neutrino burst signal; (3) The spectrum of the observed neutrino flux might be distorted; (4) If the neutrino is the Majorana particle, the ALP-neutrino coupling can cause the the excessive deleptonization which will disable the supernova explosion. However, if the coupling is so strong that the ALPs is severely trapped inside the core, no constraint can be made on the contrary. The most stringent upper bound can be placed on the coupling g ee from the energy loss argument as [81] g ee 10 −5 eV −1 0.05 eV m ν , which is much weaker than the limits of CMB. We have used Eq. (25) to represent the constraining power of the SN1987A explosion in Fig. 1.

C. Galactic Propagation of Neutrinos from SN1987A and Blazar TXS 0506+056
The propagation of astrophysical neutrino fluxes in the dark matter halo of our galaxy will be affected when the ALP-neutrino interaction is very strong [42,47,87]. These neutrinos will suffer the energy loss and the direction change, due to the frequent scattering with the ALPs which forms the dark matter. The observed neutrino flux of SN1987A agrees well with the standard supernova neutrino theory, which should in turn give us a constraint on the ALP-neutrino coupling. On the other hand, the recent UHE neutrino event IceCube-170922A observed by IceCube coincides with a flaring blazar TXS 0506+056 with ∼ 3σ level [88]. The estimated neutrino luminosities are similar to that of the associated γ-rays, consistent with the prediction of blazar models [89,90]. However, if there is a large attenuation effect for the UHE neutrino flux by scattering with ALPs, the original flux from TXS 0506+056 must be much stronger than the estimated one, which is inconsistent with the blazar observation. For the diffuse UHE neutrinos, the observed flux is around the Waxman-Bahcall (WB) bound [91]. The large attenuation effect of ALP would require a diffuse flux to be much larger than the WB bound, which can in turn set a limit on the ALP-neutrino coupling. However, the original diffuse UHE neutrino flux might be well above the WB bound, since the complete sources of the UHE neutrino still remain unknown for us.
A very simple argument would be that the mean free path of neutrinos λ ν should be smaller than the radius of the Milky Way R MW ∼ 26 kpc. However, this is not true in general because the energy loss for each ALP-neutrino collision is unclear, which might be negligible compared with the initial neutrino energy in some case. For the scattering process ν i + a → ν j + a, the average relative energy loss and angular change under one collision are estimated as where s ≈ 2m a E ν + m 2 i + m 2 a is the square of the centerof-mass energy. ∆E ν /E ν and ∆θ are defined as the halves of their possible maximum values for the collision. One can easily find ∆θ ∆E ν /E ν for our case, so we focus on the energy loss of neutrinos. Assuming 2m a E ν m 2 ν such that s ≈ 2m a E ν , we can have ∆E ν /E ν ≈ O(1). In this case, the energy loss of one single collision is significant, so one can adopt the argument λ ν R MW . However, in the case of 2m a E ν m 2 ν with i = j, one would have a negligible energy loss for each collision. It is difficult to transfer energy to a static ALP when the ALP mass is too small. Integrating over the energy loss along the neutrino course in our galaxy and in order not to violate the observations of SN1987A and blazar TXS 0506+056, a general requirement would be where λ ν = (σn a ) −1 . Here σ ∼ g 4 m 2 ν is the approximated cross section for ν + a ↔ ν + a, which should be good enough for an order-of-magnitude estimation. The general bound can be formally written as For an ALP mass m a m 2 ν /(2E ν ), one can obtain the following constraint . Note that the constraint given by IceCube is almost two orders of magnitude stronger than that of the supernova under the same argument. This is simply due to the higher energy of the IceCube event (∼ 290 TeV) compared with the supernova one (∼ 30 MeV). The constraint with m ν ≈ 0.05 eV over the whole ALP mass range for IceCube can be found in Fig. 1 as the magenta curve. The corresponding supernova constraint along with the explosion one in the last section are shown as the purple curves in Fig. 1.

IV. CONCLUSION
Assuming the derivative coupling between the ALP and the neutrinos, we have investigated the impacts of the ALP dark matter on neutrino oscillation experiments. Depending on the mass of the ALP, there are two different scenarios regarding the data analysis: the modulation effect induced by the oscillating ALP field can be resolved for the experiment (Non-Averaged ); the modulation effect is simply averaged to a distortion effect (Averaged ). Based on the simple argument, we find that the existing experiment like T2K can already exclude g 3 × 10 −10 eV −1 for 10 −22 eV m a 10 −11 eV. The projected experiments like DUNE are sensitive to the coupling strength g 10 −12 eV −1 in the ALP mass range 10 −22 eV m a 10 −12 eV for the non-averaged case. The 1σ sensitivity is reduced by a factor of O(10) for the averaged case due to the cancellation of the oscillating ALP-induced potential. Using the GLoBES package, we have numerically simulated the data of DUNE. The sensitivity results on the six coupling parameters g αβ for α, β = (e, µ, τ ) have been yielded and summarized in Table II. The numerical results agree well in the order of magnitude with the simple estimation for DUNE. The impact of such coupling on the evolution of the early Universe has been discussed. A very stringent bound from the free-streaming of CMB can be made as g ij 10 −10 eV −1 . The propagation of neutrinos from SN1987A can put a constraint g 4 × 10 −8 eV −1 , while the IceCube observation can put a much stronger one g 7 × 10 −10 eV −1 . The next generation of neutrino experiments can probe the parameter range two orders of magnitude beyond the astrophysical limits.