Search for three body pion decays π^+{\to}l^+νX
A. Aguilar-Arevalo, M. Aoki, M. Blecher, D.I. Britton, D. vom Bruch, D.A. Bryman, S. Chen, J. Comfort, S. Cuen-Rochin, L. Doria, P. Gumplinger, A. Hussein, Y. Igarashi, S. Ito, S. Kettell, L. Kurchaninov, L.S. Littenberg, C. Malbrunot, R.E. Mischke, T. Numao, D. Protopopescu, A. Sher, T. Sullivan, D. Vavilov
aa r X i v : . [ h e p - e x ] F e b Search for three body pion decays π + → l + νX A. Aguilar-Arevalo, M. Aoki, M. Blecher, D.I. Britton, D. vom Bruch, ∗ D.A. Bryman,
5, 6
S. Chen, J. Comfort, S. Cuen-Rochin,
6, 9
L. Doria,
6, 10
P. Gumplinger, A. Hussein,
6, 11
Y. Igarashi, S. Ito, † S. Kettell, L. Kurchaninov, L.S. Littenberg, C. Malbrunot, ‡ R.E. Mischke, T. Numao, D. Protopopescu, A. Sher, T. Sullivan, § and D. Vavilov (PIENU Collaboration) Instituto de Ciencias Nucleares, Universidad Nacional Aut´onoma de M´exico, CDMX 04510, M´exico Department of Physics, Graduate School of Science,Osaka University, Toyonaka, Osaka, 560-0043, Japan Virginia Tech., Blacksburg, Virginia 24061, USA SUPA - School of Physics and Astronomy, University of Glasgow, Glasgow, G12-8QQ, United Kingdom Department of Physics and Astronomy, University of British Columbia, Vancouver, British Columbia V6T 1Z1, Canada TRIUMF, 4004 Wesbrook Mall, Vancouver, British Columbia V6T 2A3, Canada Department of Engineering Physics, Tsinghua University, Beijing, 100084, China Physics Department, Arizona State University, Tempe, AZ 85287, USA Universidad Aut´onoma de Sinaloa, Culiac´an, M´exico PRISMA + Cluster of Excellence and Institut f¨ur Kernphysik,Johannes Gutenberg-Universit¨at Mainz, Johann-Joachim-Becher-Weg 45, D 55128 Mainz, Germany University of Northern British Columbia, Prince George, British Columbia V2N 4Z9, Canada KEK, 1-1 Oho, Tsukuba-shi, Ibaraki, 300-3256, Japan Physics Department, Osaka University, Toyonaka, Osaka, 560-0043, Japan Brookhaven National Laboratory, Upton, NY, 11973-5000, USA (Dated: February 10, 2021)The three body pion decays π + → l + νX ( l = e, µ ), where X is a weakly interacting neutral bo-son, were searched for using the full data set from the PIENU experiment. An improved limiton Γ( π + → e + νX ) / Γ( π + → µ + ν µ ) in the mass range 0 < m X <
120 MeV/ c and a first resultfor Γ( π + → µ + νX ) / Γ( π + → µ + ν µ ) in the region 0 < m X < . c were obtained. TheMajoron-neutrino coupling model was also constrained using the current experimental result ofthe π + → e + ν e ( γ ) branching ratio. I. INTRODUCTION
The existence of massive or massless weakly interactingneutral particles ( X ) has been suggested to augment thestandard model with motivations that include providingdark matter candidates [1], explaining baryogenesis [2],revealing the origin of neutrino masses [3], and findingsolutions to the strong CP problem [4, 5] involving theaxion [6–10]. Pion and kaon decays are potential sourcesof X particles as discussed by Altmannshofer, Gori, andRobinson [11] who investigated a model with axionlikeparticles involved in pion decay π + → e + νX . Batell et al. [12] studied a model of thermal dark matter emitted inthree body meson decay π + ( K + ) → l + χφ where χ and φ are assumed to be sterile neutrinos. Light vector bosonsemitted in π + ( K + ) → l + νX decay have been discussed byDror [13]. ∗ Present address: LPNHE, Sorbonne Universit´e, Universit´e ParisDiderot, CNRS/IN2P3, Paris, France. † Corresponding author ([email protected]).Present address: Faculty of Science, Okayama University,Okayama, 700-8530, Japan. ‡ Present address: Experimental Physics Department, CERN,Gen`eve 23, CH-1211, Switzerland. § Present address: Department of Physics, University of Victoria,Victoria BC V8P 5C2, Canada.
A Nambu-Goldstone boson, the “Majoron” proposedby Gelmini and Roncadelli [14], is also a candidate ofinterest. It arises in gauge models that have a sponta-neous breaking of the baryon and lepton numbers ( B − L )global symmetry [14, 15]. In the Majoron models, neu-trino masses arise from the vacuum expectation valueof a weak isotriplet scalar Higgs boson. Barger, Ke-ung, and Pakvasa extended the Majoron model to thedecay processes of pions and kaons π + ( K + ) → l + νX viaMajoron-neutrino couplings [16]. Other related processesand models have been discussed in Refs. [17–20].Three body pion decays π + → l + νX can be investi-gated using the decay lepton energy spectra in pion de-cays. Figure 1 shows the total and kinetic energy spectraof π + → e + νX and π + → µ + νX decays assuming the de-cay products of X are invisible or have very long life-times allowing undetected escape. The signal shapeswere obtained from Eq. (12) in Ref. [12]. A previoussearch for the decay π + → e + νX was performed by Pic-ciotto et al. [21] as a byproduct of the branching ratiomeasurement R π = Γ[ π + → e + ν e ( γ )] / Γ[ π + → µ + ν µ ( γ )],where ( γ ) indicates the inclusion of radiative decays, us-ing stopped pions in an active target [22]. The upperlimit on the branching ratio was found to be R πeνX =Γ( π + → e + νX ) / Γ( π + → µ + ν µ ) . × − in the mass range m X from 0 to 125 MeV/ c . The sensitivity was limitedby statistics and the remaining background originated (MeV) e E N o r m a li z e d c oun t s X ν + e → + π (a) =0 MeV/c X m =40 MeV/c X m =80 MeV/c X m (MeV) µ T N o r m a li z e d c oun t s X ν + µ→ + π (b) =5 MeV/c X m =15 MeV/c X m =25 MeV/c X m FIG. 1. Total energy spectra of π + → e + νX and kinetic en-ergy spectra of π + → µ + νX decays. (a) π + → e + νX decay withmass m X of 0 MeV/ c (solid black), 40 MeV/ c (dotted red),and 80 MeV/ c (dashed blue). (b) π + → µ + νX decay withmass m X of 5 MeV/ c (solid black), 15 MeV/ c (dotted red),and 25 MeV/ c (dashed blue). from pion decay-in-flight ( π DIF) events. For π + → µ + νX decay, no comparable studies have been performed.In the present work, the decays π + → e + νX and π + → µ + νX were sought using the full data set of thePIENU experiment [23] corresponding to two orders ofmagnitude larger statistics than the previous experiment[21]. The analyses were based on the searches for heavyneutrinos ν H in π + → e + ν H decay [24] and π + → µ + ν H de-cay [25], and the decays π + → e + ν e ν ¯ ν and π + → µ + ν µ ν ¯ ν [26]. II. EXPERIMENT
The PIENU detector [27] shown schematically in Fig.2 was designed to measure the pion branching ratio R π = Γ[ π + → e + ν e ( γ )] / Γ[ π + → µ + ν µ ( γ )]. The decaypositron in π + → e + ν e decay has total energy E e = 69 . π + → µ + ν µ decay followed by µ + → e + ν e ¯ ν µ de-cay ( π + → µ + → e + decay chain), the decay muon has ki-netic energy T µ = 4 . CFBN π + $T*$T* /B*(cid:1)(cid:9) l (cid:10) (cid:18)(cid:17)(cid:1)DN FIG. 2. Schematic of the PIENU detector [27]. tillator of about 1 mm; the total energy of the positronin the subsequent muon decay µ + → e + ν e ¯ ν µ ranges from E e = 0 . ± c providedby the TRIUMF M13 beam line [28] was tracked by twomultiwire proportional chambers (WC1 and WC2) andtwo sets of silicon strip detectors (S1 and S2). FollowingWC2, the beam was degraded by two thin plastic scin-tillators (B1 and B2) to measure time and energy lossfor particle identification. After S2, pions stopped anddecayed at rest in the center of an 8 mm thick plasticscintillator target (B3). The pion stopping rate in B3was 5 × π + / s.Positrons from pion or muon decay were detected byanother silicon strip detector (S3) and a multiwire pro-portional chamber (WC3) located downstream of B3 toreconstruct tracks and define the acceptance. Two thinplastic scintillators (T1 and T2) were used to measurethe positron time, and its energy was measured by a48 cm (dia.) ×
48 cm (length) single crystal NaI(T ℓ )calorimeter surrounded by 97 pure CsI crystals to detectshower leakage. The energy resolution of the calorimeterfor positrons was 2.2% (FWHM) at 70 MeV.The pion and positron signals were defined by a coin-cidence of B1, B2, and B3, and a coincidence of T1 andT2, respectively. A coincidence of the pion and positronsignals within a time window of −
300 ns to 540 ns withrespect to the pion signal was the basis of the main trig-ger condition. This was prescaled by a factor of 16 toform an unbiased trigger (Prescaled trigger). π + → e + ν e event collection was enhanced by an early time triggerselecting all events occurring between 6 and 46 ns afterthe arrival of the pion (Early trigger). The typical triggerrate including calibration triggers was about 600 s − .To extract the energy and time information, plas-tic scintillators, silicon strip detectors and CsI crystals,and the NaI(T ℓ ) crystal were read out by 500 MHz, 60MHz, and 30 MHz digitizers, respectively. The wirechambers and trigger signals were read by multi-hittime − to − digital converters with 0.625 ns resolution [27]. III. π + → e + νX DECAYA. Event selection
The decay π + → e + νX was searched for by fitting the π + → e + ν e energy spectra after π + → µ + → e + backgroundsuppression. The cuts used for the pion selection, therejection of the extra activity in scintillators, and thesuppression of π + → µ + → e + backgrounds were the sameas for the analysis of π + → e + ν e ν ¯ ν decay [26]. Pions wereidentified using the energy loss information in B1 andB2. Events with extra activity in B1, B2, T1 or T2 wererejected. Since the calibration system for the CsI crystalswas not available before November 1, 2010, the data weredivided into two sets (dataset 1, before, and dataset 2,after November 1, 2010). A 15% solid angle cut was usedfor the dataset 2, and a tighter cut (10%) was applied tothe dataset 1 to minimize the effects of electromagneticshower leakage.The π + → µ + → e + backgrounds were suppressed usingdecay time, energy in the target, and tracking infor-mation provided by WC1, WC2, S1, and S2 [25, 26].Events were first selected by the Early trigger and adecay time cut t = 7 −
35 ns after the pion stop wasapplied. The energy loss information in B3 was usedbecause π + → µ + → e + backgrounds deposit larger energyin B3 than π + → e + ν e decays due to the presence of thedecay muon ( T µ = 4 . π + → e + ν e events (mostly, π DIF eventsbefore B3) [27]. Figure 3 shows the decay positron en-ergy spectra of π + → e + ν e decays after π + → µ + → e + back-ground suppression cuts ((a) dataset 1 and (c) dataset 2).The bumps in the positron energy spectra at about 58MeV are due to photo-nuclear reactions in the NaI(T ℓ )[30]. The total number of π + → e + ν e events was 1 . × (5 × in dataset 1 and 8 × in dataset 2). B. Energy spectrum fit
The energy spectrum was fitted with a combinationof background terms and a shape to represent the sig-nal. The background component due to the remaining π + → µ + → e + events was obtained from the data by re-quiring a late time region t >
200 ns. The shape of thelow energy π + → e + ν e tail was obtained by Monte Carlo(MC) simulation [31] including the detector responsewhich was measured using a mono-energetic positronbeam [27, 30]. Because the solid angle cut was reducedand the CsI was not used for dataset 1, the shapes of thelow energy π + → e + ν e tails are slightly different for thetwo datasets. Another background came from the decays-in-flight of muons ( µ DIF) following π + → µ + ν µ decays inB3 that has a similar time distribution to π + → e + ν e de-cay. The shape of the µ DIF event spectrum was obtainedby MC simulation. The signal shapes as shown in Fig.1 (a) were produced with mass range m X from 0 to 120MeV/ c in 5 MeV/ c steps by MC simulation includingthe detector response. These shapes were normalized to 1and used for the fit to search for the signals. To combinethe two data sets, simultaneous fitting with a commonbranching ratio as a free parameter was performed. Thefit in the range of E e = 5 −
56 MeV without any signalresulted in χ / d.o.f.=1.04 (d.o.f.=402). The addition ofthe signals did not change the fit result. C. Results
Figure 3 (b) and (d) show the residual plots with-out any signal in datasets 1 and 2; hypothetical signalsassuming m X = 80 MeV/ c with the branching ratio R πeνX = 2 . × − are also shown. No significant ex-cess above the statistical uncertainty was observed. Forexample, the branching ratio with m X = 0 MeV/ c ob-tained by the fit was R πeνX = (0 . ± . × − . Figure 4shows the 90% confidence level (C.L.) upper limits for thebranching ratio π + → e + νX in the mass region from 0 to120 MeV/ c calculated using the Feldman and Cousins(FC) approach [32]. Since the signal shape at a massof 55 MeV/ c is similar to the π + → µ + → e + energy spec-trum, the sensitivity was worse than for other masses dueto the strong correlation; R πeνX = ( − . ± . × − .The statistical uncertainty dominates because the sys-tematic uncertainties and the acceptance effects are ap-proximately canceled out by taking the ratio of the num-ber of signal events obtained by the fit to the number ofpion decays. The acceptance effect due to the cuts wasexamined by generating positrons in B3 isotropically withan energy range of E e = 0 −
70 MeV using the MC simula-tion and the systematic uncertainty was estimated to be < (MeV) e E10 20 30 40 50 60 70 C oun t s / . M e V
10 (a) DataFitting function e ν + e → + π + e → + µ→ + π DIF µ (MeV) e E10 20 30 40 50 D a t a - F i t -100-50050100 (b) (MeV) e E10 20 30 40 50 60 70 C oun t s / . M e V
10 (c) DataFitting function e ν + e → + π + e → + µ→ + π DIF µ (MeV) e E10 20 30 40 50 D a t a - F i t -150-100-50050100150 (d) FIG. 3. First and third panels from the top: the E e spectra of π + → e + ν e decay after π + → µ + → e + suppression cuts for datasets1 (a) and 2 (c). The black crosses with the statistical uncertainties show the data. Background components illustrated by thedashed and dotted green line, dotted blue line, dashed gray line, and solid red line represent π + → µ + → e + decays, low energy π + → e + ν e tail, µ DIF events, and the sum of those three components, respectively (see text). Second and fourth panels fromthe top: the residual plots shown by the black circles with statistical error bars and hypothetical signals (solid red lines) witha mass of m X = 80 MeV/ c and a branching ratio R πeνX = 2 . × − from datasets 1 (b) and 2 (d) (the branching ratioobtained by the fit at this mass was R πeνX = ( − . ± . × − ). IV. π + → µ + νX DECAY
The decay π + → µ + νX can be sought by a measure-ment of the muon kinetic energy in π + → µ + ν decay (fol-lowed by µ + → e + ν e ¯ ν µ decay) in the target (B3). In the π + → µ + → e + decay chain, three hits are expected in B3:the first signal is from the beam pion, the second is fromthe decay muon, and the third is from the decay positron.Thus, the second of three pulses in B3 would be due tothe muon kinetic energy. However, the pulse detectionlogic could not efficiently identify pulses below 1.2 MeV[24]. Therefore, the search was divided into two muon en-ergy regions, above and below 1.2 MeV. The number of Prescaled trigger events used for the analysis was 4 × .The analysis strategy and event selection cuts were basedon the massive neutrino [25] and three neutrino decay [26]searches, briefly described in the following sections. A. Analysis of the region above 1.2 MeV
As described in Sec. III A, pions were identified usingB1 and B2 and events with extra hits in B1, B2, T1, orT2 were rejected. A solid angle acceptance of about 20%for the decay positron was used. To ensure the selectedevents were from π + → µ + → e + decays, a late positron de- ) (MeV/c M Mass m0 20 40 60 80 100 120 ( % C . L . ) X ν e π U pp e r li m i t o f b r a n c h i ng r a t i o R -8 -7 -6 -5
10 Picciotto et al. [21]This work
FIG. 4. Results of the 90% C.L. upper limit branching ratio R πeνX . Dashed black line: previous TRIUMF results [21].Solid red line with filled circles: results from this work. (MeV) µ T C oun t s / . M e V (a) DataSum of the functionsGaussian function γ µ ν + µ→ + π (MeV) µ T D a t a - F i t -20-15-10-505101520 (b) FIG. 5. (a) The T µ spectra of π + → µ + → e + decay. The blackcrosses with the statistical uncertainties show the data. Thedotted green line, dashed blue line, and solid red line repre-sent a Gaussian distribution centered at 4.1 MeV, π + → µ + ν µ γ decay, and the sum of those two functions, respectively. (b)Residual plots shown by the black circles with statistical errorbars in the range T µ =1.3 to 3.4 MeV. The solid red line rep-resents a hypothetical signal with mass of m X = 15 MeV/ c and the branching ratio R πµνX = 6 . × − ; the branchingratio obtained by the fit was R πµνX = ( − . ± . × − . cay time t >
200 ns after the pion stop and the positronenergy in the NaI(T ℓ ) calorimeter E e <
55 MeV wererequired. Then, the events with three clearly separatedpulses in the target (B3) were selected and the secondpulse information was extracted and assigned to the de-cay muon [24]. The muon kinetic energy ( T µ ) spectrumafter the event selection cuts is shown in Fig. 5 (a).As described above, the drop below 1.2 MeV was dueto the inefficiency of the pulse detection logic [24]. Themain background below 3.4 MeV was due to the radia-tive pion decay π + → µ + ν µ γ (branching fraction 2 × − [29]). The total number of π + → µ + → e + events availablewas 9.1 × .The decay π + → µ + νX was searched for by fitting the T µ energy spectrum of π + → µ + → e + decays. The fit wasperformed using a Gaussian peak centered at 4.1 MeV(energy resolution σ = 0 .
16 MeV), the π + → µ + ν µ γ de-cay spectrum obtained by MC simulation [31], and thenormalized signal spectra including the energy resolutionin B3. The signal spectra as shown in Fig. 1 (b) weregenerated with the mass range 0 < m X <
26 MeV/ c with 1 MeV/ c steps using MC including detector res-olution. The fit for T µ from 1.3 to 4.2 MeV withoutany π + → µ + νX signal introduced gave χ / d.o.f.=1.27(d.o.f.=53) and the residuals of the fit for the signal sen-sitive region are shown in Fig. 5 (b). The addition ofsignal components did not change the fit result.No significant signal beyond the statistical uncer-tainty was observed. For example, the branching ra-tios for the signals with mass m X = 0 MeV/ c and 26 MeV/ c obtained by the fit were R πµνX =Γ( π + → µ + νX ) / Γ( π + → µ + ν µ ) = ( − . ± . × − and( − . ± . × − , respectively. Systematic uncertain-ties and acceptance effects were approximately canceledby taking the ratio of amplitudes for the signal and π + → µ + ν µ decays. The systematic uncertainties and ac-ceptance effects due to the cuts were examined by gen-erating decay muons in the target with several kineticenergies in the range T µ = 0 − . < R πµνX in thisenergy region calculated using the FC method. B. Analysis of the region below 1.2 MeV
For T µ < . ℓ ) calorimeter were all the same as in the analysisin the energy region T µ > . π DIFevents, the same tracking cut by WC1, WC2, S1, and S2used in Sec. III A was also applied. After these basiccuts, the energies observed in B3 in a wide time win-dow (700 ns) including pion and positron energies wereobtained. To cleanly subtract the positron contributionfrom the integrated energy, events with late positron de- ) (MeV/c X Mass m0 5 10 15 20 25 30 35 ( % C . L . ) X ν µ π U pp e r li m i t o f b r a n c h i ng r a t i o R -6 -5 -4 >1.2 MeV µ Analysis for T <1.2 MeV µ Analysis for T
FIG. 6. Summary of the 90% C.L. upper limit branching ratio R πµνX in this work. The black circles show the result of thesearch in the energy region T µ > . T µ < . (MeV)) µ Total energy after subtraction of 17 MeV (T-6 -4 -2 0 2 4 6 C oun t s / . M e V DataSum of the functionsMain peak at 4.1 MeVQuadratic function (a) (MeV)) µ Total energy after subtraction of 17 MeV (T -1.5 -1 -0.5 0 0.5 1 1.5 D a t a - F i t -150-100-50050100150200 (b) FIG. 7. (a) The total energy in the target due to the pionand muon after subtracting 17 MeV. The black crosses withstatistical uncertainties show the data. The dotted green line,dashed blue line, and solid red line represent the main peak at4.1 MeV, quadratic background due to π DIF events, and thesum of those two functions, respectively. (b) Residual plotsshown by the black circles with the statistical error bars in thesignal region T µ =-1.8 to 1.8 MeV. The solid red line representsa hypothetical signal with mass of m X = 33 . c andthe branching ratio R πµνX = 3 . × − . cay t >
300 ns were selected and the isolated positronenergy was subtracted. After that, the contribution ofthe averaged pion kinetic energy ( ∼
17 MeV) was sub-tracted from the total energy (due to the pion and themuon). Figure 7 (a) shows the total energy (correspond-ing to T µ ) after subtracting 17 MeV. The backgroundbelow T µ < π DIFevents. The number of π + → µ + → e + events available forthe analysis is 1 . × .There are two background shapes, the 4.1 MeV peakand the π DIF events. A quadratic function was usedfor the π DIF events. To search for π + → µ + νX decay,the width of the signal shape was scaled using that atthe 4.1 MeV peak. Figure 7 (b) shows the residual plotsin the signal region from -1.8 to 1.8 MeV without anysignal shape and a hypothetical signal shape assuminga mass of m X = 33 . c with the branching ratio R πµνX = 3 . × − . The branching ratio obtained by thefit was (1 . ± . × − . The fit was performed from -4.0to 4.1 MeV and the fitting range of -4.0 to 2.0 MeV (sig-nal region) resulted in χ / d.o.f.=1.03 (d.o.f.=115); thereis some small deviation above 2 MeV due to a small mis-match due to the kinetic energy distribution of the beampion.The signals of π + → µ + νX decay were searched for inthe mass range of m X = 26 to 33.9 MeV/ c , but nosignificant excess beyond the statistical uncertainty wasobserved. The red squares in Fig. 6 represent the resultof the 90% C.L. upper limit branching ratio R πµνX inthis energy region calculated using the FC approach. V. CONSTRAINTS ON THE MAJORONMODEL
The Majoron model can be constrained using the ex-perimental value of the pion branching ratio R π . Thepredicted branching ratio including the massless Majoron X and a light neutral Higgs H ′ ( . c ) can bewritten asΓ( π → eL ) / Γ( π → µL )Γ( π → eν e ) / Γ( π → µν µ ) = 1 + 157 . g (1)where L is the final state ν , νX , and νH ′ , and g is the Majoron-neutrino coupling constant [16]. Theupper limit of the ratio R π exp /R π SM at 90% C.L. us-ing the current averaged experimental value R π exp =(1 . ± . × − [33] is R π exp R π SM < . . (2)Using this limit, the 90% C.L. upper limit of the couplingconstant can be found to be g < × − , (3)which was improved by a factor of three over the previousexperiment [22]. VI. CONCLUSION
No evidence of the three body pion decays π + → e + νX or π + → µ + νX was found and new upper limits were set.The limits on the branching ratio π + → e + νX were im-proved by an order of magnitude over the previous ex-periment. For π + → µ + νX decay, the limits obtained arethe first available results. The Majoron model was alsoconstrained using the pion branching ratio R π . ACKNOWLEDGMENTS
This work was supported by the Natural Sciences andEngineering Research Council of Canada (NSERC, No. SAPPJ-2017-00033), and by the Research Fund for theDoctoral Program of Higher Education of China, byCONACYT doctoral fellowship from Mexico, and byJSPS KAKENHI Grant No. 18540274, No. 21340059,No. 24224006, and No. 19K03888 in Japan. We aregrateful to Brookhaven National Laboratory for the loanof the crystals, and to the TRIUMF operations, detec-tor, electronics and DAQ groups for their engineering andtechnical support. [1] G. Bertone, D. Hooper, and J. Silk, Phys. Rep. , 279(2005).[2] A.D. Dolgov, arXiv:hep-ph/9707419; V.A. Rubakov andM.E. Shaposhnikov, Phys. Usp. , 461 (1996).[3] Y. Fukuda et al ., Phys. Rev. Lett. , (1998) 1562.[4] R.D. Peccei and H.R. Quinn, Phys. Rev. Lett. (1977)1440.[5] R.D. Peccei and H.R. Quinn, Phys. Rev. D (1977)1791.[6] F. Wilczek, Phys. Rev. Lett. , 1549 (1982); see alsoA. Davidson and K. C. Wali, Phys. Rev. Lett. 48, 11(1982).[7] J. Jaeckel and A. Ringwald, Annu. Rev. Nucl. Part. Sci. , 405 (2010).[8] P. Agrawal and K. Howe, J. High Energy Phys. (2018)029.[9] D.S. M. Alves and N. Weiner, J. High Energy Phys. (2018) 092.[10] K.S. Jeong, T.H. Jung, and C.S. Shin, Phys. Rev. D ,035009 (2020).[11] W. Altmannshofer, S. Gori, and D.J. Robinson, Phys.Rev. D , 075002 (2020).[12] B. Batell, T. Han, D. McKeen, and B.S.E. Haghi, Phys.Rev. D , 075016 (2018).[13] J.A. Dror, Phys. Rev. D , 411(1981); see also G.B. Gelmini, S. Nussinov, and M. Ron-cadelli, Nucl. Phys. B209 (1982) 157-173.[15] Y. Chikashige, R. N. Mohapatra, and R. D. Peccei, Phys.Lett. , 265 (1981).[16] V. Barger, W.Y. Keung, and S. Pakvasa, Phys. Rev. D , 907 (1982). [17] A. Masiero, J.W.F. Valle, Phys. Lett. B251 , 273-278(1990).[18] A.P. Lessa and O.L.G. Peres, Phys. Rev. D , 094001(2007).[19] M. Hirsch, A. Vicente, J. Meyer, and W. Porod, Phys.Rev. D , 055023 (2009).[20] X. Garcia i Tormo, D. Bryman, A. Czarnecki, and M.Dowling, Phys. Rev. D , 113010 (2011).[21] C.E. Picciotto et al ., Phys. Rev. D , 1131 (1988).[22] D.I. Britton et al ., Phys. Rev. Lett. , 3000 (1992) andPhys. Rev. D , 28 (1994).[23] A. Aguilar-Arevalo et al ., Phys. Rev. Lett. , 071801(2015).[24] M. Aoki et al ., Phys. Rev. D , 052002 (2011) and A.Aguilar-Arevalo et al ., Phys. Rev. D , 072012 (2018).[25] A. Aguilar-Arevalo et al ., Phys. Lett. B , 134980(2019).[26] A. Aguilar-Arevalo et al ., Phys. Rev. D , 012001(2020).[27] A. Aguilar-Arevalo et al ., Nucl. Instrum. Methods Phys.Res., Sect. A , 38 (2015).[28] A. Aguilar-Arevalo et al ., Nucl. Instrum. Methods Phys.Res., Sect. A , 102 (2009).[29] G. Bressi, G. Carugno, S. Cerdonio, E. Conti,A.T. Meneguzzo, and D. Zanello, Nucl. Phys. B (1998) 555.[30] A. Aguilar-Arevalo et al. , Nucl. Instrum. Methods Phys.Res., Sect. A , 188 (2010).[31] S. Agostinelli et al . (GEANT4 Collaboration), Nucl. In-strum. Methods Phys. Res., Sect. A , 250 (2003);http://geant4.cern.ch.[32] G.J. Feldman and R.D. Cousins, Phys. Rev. D , 3873(1998).[33] M. Tanabashi et al. (Particle Data Group), Phys. Rev.D98