Seasonal Modulation of the $^7$Be Solar Neutrino Rate in Borexino
M. Agostini K. Altenmuller S. Appel V. Atroshchenko D. Basilico G. Bellini J. Benziger D. Bick G. Bonfini L. Borodikhina D. Bravo B. Caccianiga F. Calaprice A. Caminata S. Caprioli M. Carlini P. Cavalcante A. Chepurnov K. Choi D. D'Angelo S. Davini A. Derbin X. F. Ding L. Di Noto I. Drachnev K. Fomenko D. Franco F. Froborg F. Gabriele C. Galbiati C. Ghiano M. Giammarchi M. Goeger-Neff A. Goretti M. Gromov C. Hagner T. Houdy E. Hungerford Aldo Ianni Andrea Ianni A. Jany D. Jeschke V. Kobychev D. Korablev G. Korga D. Kryn M. Laubenstein B. Lehnert E. Litvinovich F. Lombardi P. Lombardi L. Ludhova G. Lukyanchenko I. Machulin S. Manecki G. Manuzio S. Marcocci J. Martyn E. Meroni M. Meyer L. Miramonti M. Misiaszek M. Montuschi V. Muratova B. Neumair L. Oberauer B. Opitz F. Ortica M. Pallavicini L. Papp A. Pocar G. Ranucci A. Razeto A. Re A. Romani R. Roncin N. Rossi S. Schonert D. Semenov P. Shakina M. Skorokhvatov O. Smirnov A. Sotnikov L. F. F. Stokes Y. Suvorov R. Tartaglia G. Testera J. Thurn M. Toropova E. Unzhakov A. Vishneva R. B. Vogelaar F. von Feilitzsch H. Wang S. Weinz M. Wojcik M. Wurm Z. Yokley O. Zaimidoroga S. Zavatarelli K. Zuber G. Zuzel








GSSI]M. Agostini
Munchen]K. Altenmüller
Munchen]S. Appel
Kurchatov]V. Atroshchenko
Milano]D. Basilico
Milano]G. Bellini
PrincetonChemEng]J. Benziger
Hamburg]D. Bick
LNGS]G. Bonfini
Kurchatov]L. Borodikhina
Virginia,Milano]D. Bravo
Milano]B. Caccianiga
Princeton]F. Calaprice
Genova]A. Caminata
Milano]S. Caprioli
LNGS]M. Carlini
LNGS,Virginia]P. Cavalcante
Lomonosov]A. Chepurnov
Honolulu]K. Choi
Milano]D. D'Angelo
GSSI,Genova]S. Davini
Peters]A. Derbin
GSSI]X.F. Ding
Genova]L. Di Noto
GSSI,Peters]I. Drachnev
Dubna]K. Fomenko
APC]D. Franco
Princeton]F. Froborg
LNGS]F. Gabriele
Princeton]C. Galbiati
Genova]C. Ghiano
Milano]M. Giammarchi
Munchen]M. Goeger-Neff
Princeton]A. Goretti
Lomonosov]M. Gromov
Hamburg]C. Hagner
APC]T. Houdy
Huston]E. Hungerford
LNGS]Aldo IanniCanfranc
Princeton]Andrea Ianni
Krakow]A. Jany
Munchen]D. Jeschke
Kiev]V. Kobychev
Dubna]D. Korablev
Huston]G. Korga
APC]D. Kryn
LNGS]M. Laubenstein
Dresda]B. Lehnert
Kurchatov,Kurchatovb]E. Litvinovich
LNGS]F. LombardiDIEGO
Milano]P. Lombardi
Juelich,RWTH]L. Ludhova
Kurchatov]G. Lukyanchenko
Kurchatov,Kurchatovb]I. Machulin
Virginia]S. ManeckiQueens
Genova]G. Manuzio
GSSI]S. Marcocci
Mainz]J. Martyn
Milano]E. Meroni
Hamburg]M. Meyer
Milano]L. Miramonti
Krakow]M. Misiaszek
Ferrara]M. Montuschi
Peters]V. Muratova
Munchen]B. Neumair
Munchen]L. Oberauer
Hamburg]B. Opitz
Perugia]F. Ortica
Genova]M. Pallavicini
Munchen]L. Papp
UMass]A. Pocar
Milano]G. Ranucci
LNGS]A. Razeto
Milano]A. Re
Perugia]A. Romani
LNGS,APC]R. Roncin
LNGS]N. Rossi
Munchen]S. Schönert
Peters]D. Semenov
Peters]P. Shakina
Kurchatov,Kurchatovb]M. Skorokhvatov
Dubna]O. Smirnov
Dubna]A. Sotnikov
LNGS]L.F.F. Stokes
UCLA,Kurchatov]Y. Suvorov
LNGS]R. Tartaglia
Genova]G. Testera
Dresda]J. Thurn
Kurchatov]M. Toropova
Peters]E. Unzhakov
Dubna]A. Vishneva
Virginia]R.B. Vogelaar
Munchen]F. von Feilitzsch
UCLA]H. Wang
Mainz]S. Weinz
Krakow]M. Wojcik
Mainz]M. Wurm
Virginia]Z. Yokley
Dubna]O. Zaimidoroga
Genova]S. Zavatarelli
Dresda]K. Zuber
Krakow]G. Zuzel 

[cc]Corresponding author: spokeperson-borex@lngs.infn.it
[Canfranc]Also at: Laboratorio Subterráneo de Canfranc, Paseo de los Ayerbe S/N, 22880 Canfranc Estacion Huesca, Spain
[DIEGO]Present address: Physics Department, University of California, San Diego, CA 92093, USA
[Queens]Present address: Physics Department, Queen's University, Kingston ON K7L 3N6, Canada

The Borexino Collaboration


[APC]AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/IRFU, Observatoire de Paris, Sorbonne Paris Cité, 75205 Paris Cedex 13, France
[Dubna]Joint Institute for Nuclear Research, 141980 Dubna, Russia
[Genova]Dipartimento di Fisica, Università degli Studi e INFN, 16146 Genova, Italy
[Hamburg]Institut für Experimentalphysik, Universität Hamburg, 22761 Hamburg, Germany
[Heidelberg]Max-Planck-Institut für Kernphysik, 69117 Heidelberg, Germany
[Krakow]M. Smoluchowski Institute of Physics, Jagiellonian University, 30059 Krakow, Poland
[Kiev]Kiev Institute for Nuclear Research, 03680 Kiev, Ukraine
[Kurchatov]National Research Centre Kurchatov Institute, 123182 Moscow, Russia
[Kurchatovb] National Research Nuclear University MEPhI (Moscow Engineering Physics Institute), 115409 Moscow, Russia
[LNGS]INFN Laboratori Nazionali del Gran Sasso, 67010 Assergi (AQ), Italy
[Milano]Dipartimento di Fisica, Università degli Studi e INFN, 20133 Milano, Italy
[Perugia]Dipartimento di Chimica, Biologia e Biotecnologie, Università degli Studi e INFN, 06123 Perugia, Italy
[Peters]St. Petersburg Nuclear Physics Institute NRC Kurchatov Institute, 188350 Gatchina, Russia
[Princeton]Physics Department, Princeton University, Princeton, NJ 08544, USA
[PrincetonChemEng]Chemical Engineering Department, Princeton University, Princeton, NJ 08544, USA
[UMass]Amherst Center for Fundamental Interactions and Physics Department, University of Massachusetts, Amherst, MA 01003, USA
[Virginia]Physics Department, Virginia Polytechnic Institute and State University, Blacksburg, VA 24061, USA
[Ferrara]Dipartimento di Fisica e Scienze della Terra  Università degli Studi di Ferrara e INFN,  Via Saragat 1-44122, Ferrara, Italy
[Munchen]Physik-Department and Excellence Cluster Universe, Technische Universität  München, 85748 Garching, Germany
[Lomonosov] Lomonosov Moscow State University Skobeltsyn Institute of Nuclear Physics, 119234 Moscow, Russia
[GSSI] Gran Sasso Science Institute (INFN), 67100 L'Aquila, Italy
[Huston]Department of Physics, University of Houston, Houston, TX 77204, USA
[Dresda]Department of Physics, Technische Universität Dresden, 01062 Dresden, Germany
[UCLA]Physics and Astronomy Department, University of California Los Angeles (UCLA), Los Angeles, California 90095, USA
[Mainz]Institute of Physics and Excellence Cluster PRISMA, Johannes Gutenberg-Universität Mainz, 55099 Mainz, Germany
[Honolulu]Department of Physics and Astronomy, University of Hawaii, Honolulu, HI 96822, USA
[Juelich]IKP-2 Forschungzentrum Jülich, 52428 Jülich, Germany
[RWTH]RWTH Aachen University, 52062 Aachen, Germany









We present the evidence for the seasonal modulation of the ^7Be neutrino interaction rate with the Borexino detector at the Laboratori Nazionali del Gran Sasso in Italy. 
The period, amplitude, and phase of the observed time evolution of the signal are 
consistent with its solar origin, and the absence of 
an annual modulation is rejected at 99.99% C.L. The data are analyzed using three methods: the analytical fit to event rate, the Lomb-Scargle and the Empirical Mode Decomposition techniques, which all yield results in excellent agreement.





Solar neutrinos; neutrino oscillations; liquid scintillators detectors; low background detectors. 










§ INTRODUCTION



Since 2007 Borexino <cit.> has measured the fluxes of low-energy neutrinos, most notably those emitted in nuclear fusion reactions and β decays along the pp-chain in the Sun. Borexino was the first experiment to make spectroscopic and real-time measurements of solar neutrinos with energy <3 MeV, i.e. below the endpoint energy of long-lived, natural β radioactivity: ^40K and the ^232Th and ^238U decay chains.
The detector has made first direct observations of ^7Be <cit.>, pep <cit.>, and pp <cit.> solar neutrinos, lowered the detection threshold for ^8B solar neutrinos <cit.>. These measurements deepen our understanding of Solar Standard Model <cit.> and support the MSW-LMA mechanism of neutrino oscillations.  In addition Borexino has detected anti-neutrinos from the Earth and distant nuclear reactors <cit.> and has set a new upper limit for a hypothetical solar anti-neutrinos flux <cit.>.

Borexino, located deep underground (3,800 m water equivalent) in Hall C of the Gran Sasso Laboratory (Italy), measures solar neutrinos via their interactions with a target of 278 ton organic liquid scintillator.
The ultrapure liquid scintillator (pseudocumene (1,2,4-trimethylbenzene (PC)) solvent with 1.5 g/l 2,5-diphenyloxazole (PPO) scintillating solute) is contained inside a thin transparent spherical nylon vessel of 8.5 m diameter. Solar neutrinos are detected by measuring the energy and position of electrons scattered by neutrino-electron elastic interactions. The scintillator promptly converts the kinetic energy of electrons by emitting photons, which are detected and converted into electronic signals (photoelectrons (p.e.)) by 2,212 photomultipliers (PMT) mounted on a concentric 13.7 m-diameter stainless steel sphere (SSS).

The volume between the nylon vessel and the SSS is filled with 889 ton of ultra pure, non scintillating fluid and acts as a radiation shield for external gamma rays and neutrons.
A second, larger nylon sphere (11.5 m diameter) prevents radon and other radioactive contaminants from the PMTs and SSS from diffusing into the central sensitive volume of the detector. The SSS is immersed in a 2,100 ton water Čerenkov detector meant to detect residual cosmic muons <cit.>.

Radioactive decays within the scintillator form a background that can mimic neutrino signals. During detector design and construction, a significant effort was made to minimize the radioactive contamination of the scintillator and of all detector components in contact with it. A record low scintillator contamination of <10^-18 g/g was achieved for ^238U and ^232Th.

The identification of different components of the solar neutrino flux relies on fitting the recorded energy spectrum with a combination of identified radioactive background components and of solar neutrino-induced electron recoil spectra. The neutrino-induced spectra are derived from Standard Solar Model neutrino energy distributions (SSM <cit.>) and include the effect of neutrino oscillation.
The solar origin of the detected neutrinos is determined by the identification of crisp spectral signatures as predicted by the SSM. Exemplary is the Compton-like energy spectrum of electrons scattered by the mono-energetic ^7Be solar neutrinos. Remarkably, the ^7Be-induced Compton 'shoulder' was clearly identified with just one month of data <cit.>, thanks to the extremely low radioactive background in the scintillator.

In contrast with water Čerenkov detectors, Borexino cannot retain directional information of individual events due to the isotropic emission of scintillation light; direct solar imaging with neutrinos is thus not possible. The eccentricity of the Earth's orbit, however, induces a modulation of the detected solar neutrino interaction rate proportional in amplitude to the solid angle subtended by the Earth with respect to the Sun (neglecting neutrino oscillation effects). The effect appears as a 6.7% peak-to-peak seasonal amplitude modulation, with a maximum at the perihelion. Evidence for such a yearly modulation of the ^7Be signal was already observed with Borexino Phase-I data (collected from May 2007 to May 2010) <cit.>.
The period and phase were found to be consistent with a solar origin of the signal. 

Yearly modulation searches have also been carried out by other solar neutrino experiments: in particular SNO <cit.> and Super-Kamiokande <cit.> found evidence for an annual flux modulation in their time series datasets.
Similar analyses were also performed aiming to search for time-dependencies of solar neutrino rates with periods other than one year. An apparent anti-correlation with solar cycles was suggested by data from the Homestake chlorine experiment <cit.>, and claims of such a periodicity were reported for Super-Kamiokande-I  <cit.>.
The SNO <cit.>, Super-Kamiokande <cit.>, and Gallex/GNO <cit.> collaborations looked for these time variations, but found none in their data.

Here we report an improved measurement of time periodicities of the ^7Be solar neutrino rate based on 4 years of Borexino Phase-II data, acquired between December 2011 and December 2015.
Borexino Phase-II began immediately after an extensive period of scintillator purification. 
Borexino Phase-II,  in addition to higher statistics, lower background levels and an improved rejection of alpha-decay background, is characterized by the absence of major scintillator handling and thus displays a high degree of stability of the detector, crucially important for identifying time dependent signals. 
In the Borexino Phase-I analysis we based our annual modulation search on the well-established Lomb-Scargle approach as well as on the more recent Empirical Mode Decomposition (EMD) technique. The virtue of the latter technique is its sensitivity to transient modulations embedded in time series, emerging from analyzing data features with more than just standard reference sinusoidal functions.

The analysis reported here analyzes the Borexino Phase-II dataset, described in Sec. <ref>, by employing both the Lomb-Scargle and an updated version of the EMD techniques. Two independent sections of this paper describe the methods of each approach and their respective results (Sec. <ref> and Sec. <ref>).
For completeness, we have also carried out a search of the annual modulation directly in the time domain, using a straightforward analytical fit (Sec. <ref>).
All analysis methods clearly confirm the presence of an annual modulation of the ^7Be solar neutrino interaction rate in Borexino and show no signs of other periodic time variations. 




§ THE DATA SET





The data of Borexino Phase-II are used for this analysis (1456 astronomical days of data). 
Compared to Borexino Phase-I, background levels have been substantially reduced by an extensive purification campaign that took place during 2010 and 2011. 
Of particular importance for this study is the reduction of the ^85Kr and ^210Bi concentrations, both backgrounds in ^7Be region.
Data taking has seen only occasional, minor interruptions due to detector maintenance.





 §.§ Event selection




A set of cuts described in <cit.> has been applied on an event-by-event basis to
remove backgrounds and non physical events. In particular, muons and spallation events within 300 ms of parent muons, time-correlated events (^214Bi-^214Po), and noise events are identified and removed. In addition, events featuring vertices reconstructed outside a Fiducial Volume (FV) are rejected.
Recoil electrons from the elastic scattering of ^7Be-ν's are selected by restricting the analysis to the energy region 
∼215-715 keV (115-380 N_pe). In this range, the major backgrounds are the α decays of ^210Po and the β decays of ^210Bi and ^85Kr.
The 5.3 MeV α's appear as a peak at ∼450 keV (after quenching) in the energy spectrum (red line in Fig. <ref>). The β's
define a continuous spectrum beneath the ^7Be recoil spectrum (blue line in Fig. <ref>).
The time stability of the background was studied to factor out any influence on the annual modulation search.
Two major changes were implemented for this search from that with Borexino Phase-I data and described below: the FV (Sec. <ref>) was redefined and an enhanced method for the rejection of ^210Po α background was developed (Sec. <ref>).







  §.§.§ Fiducial Volume Selection



We define a FV of 98.6 ton by combining a spherical cut of R=3 m radius at the center of the detector with two paraboloidal cuts at the nylon vessel poles to reject γ-rays from the Inner Vessel end-cap support hardware and plumbing. 

The excluded paraboloids have different dimensions to remove the local background. The paraboloids are defined as R(θ)=d/cos^n θ, where θ is the angle with z-axis and d is the distance from the detector center to the paraboloid vertex. The top paraboloid is defined by d=250 cm and n=12 whitch corresponds to an aperture of 54 cm of radius;  the bottom one by d=-240 cm and n=4 which corresponds to a larger aperture of 91 cm of radius.








  §.§.§ ^210Po Rejection





^210Po in the scintillator constitutes a background for the search of time-varying signals because of its decay half-life of 138 days.
In general α-backgrounds and β-events in a liquid scintillator can be efficiently separated exploiting the largely different shapes of the scintillation pulses <cit.>.
A novel pulse-shape method based on MultiLayer Perceptron (MLP) machine learning algorithm was applied to distinguish between the scintillation pulses of α and β particles with high efficiency.
This multivariate method uses a neural network based on 13 α/β 
discriminating input variables, that are computed for each event from the time distribution of reconstructed PMT hits. Clean samples of α and β events were obtained from the radon daughters ^214Po and ^214Bi to train the neural network. 
The resulting mlp parameter assumes values mostly between 0 (α) and 1 (β). 
Figure <ref> shows the distributions of the mlp parameters for the 
^214Po and ^214Bi event samples.

The MLP provides excellent α-β discrimination: with the mlp parameter threshold set at 0.9 to retain β's, the α rejection efficiency is >99.98% for ^214Po candidate events (7.7 MeV). 
The discrimination technique is based upon scintillation pulse shape, therefore we expect a reduced performance for the lower energy ^210Po α's (5.3 MeV) due to lower photoelectron statistics.  
In this case, for a clean β-like electron-recoil sample, we select events with mlp >0.98. 
Fig. <ref> shows the energy spectrum with and without α subtraction (blue and red lines). The small residual ^210Po events and the unaffected β spectrum illustrate the efficacy of the discrimination.




 §.§ Residual Background



There are two main sources of background for this analysis: the residual ^210Po activity, and the stability of ^210Bi and ^85Kr β-decays in the FV. 



  §.§.§ Residual ^210Po



At the beginning of Borexino Phase-II (Dec. 2011), the count rate of ^210Po was  
∼ 1400 cpd/100 ton. Estimating an mlp α-β efficiency of ≃ 99%, the residual α contamination of the β spectrum is R_α∼ 14 cpd/100 ton, 
comparable to an average β count rate (ν-signal and background) R_β∼ 40 cpd/100 ton distributed over the entire analysis energy region.
We estimated the efficiency of the MLP cut by looking for any exponentially
decaying ^210Po residual still present in the dataset. The residual amount of R_α has been subtracted for a given mlp cut in each time bin R(t):


    R_β(t) =R(t)-ξ_mlp· R_α(t),

where ξ_mlp is the `inefficiency' parameter.



For ξ_mlp = 1%  the exponential component due to the residual alphas become negligible in the overall time series of the dataset, leaving the remaining β's rates with a constant average value in time.




  §.§.§ Background stability


The β-decays of ^210Bi and ^85Kr cannot be distinguished from recoil electrons of the same energies induced by neutrinos.  To study the stability of the background rate over time, we compared the spectral fits to the data divided in short periods.
The fit procedure
is the same as in the ^7Be analysis <cit.>.
No appreciable variation of the background rate is observed within uncertainties.




 §.§ Detector Stability

The stability of the detector response also needs to be characterized, in particular of energy and position reconstruction and fiducial mass.



  §.§.§ Energy and Position Reconstruction

The stability of the energy scale over time was checked by comparing the number of events in the selected energy window and in the FV with those expected by Monte Carlo. A detailed simulation that includes the run per run detector performance is used. The stability of the energy scale over the period of interest was proven to be better than 1%, adequate for our purposes. 



  §.§.§ Fiducial Mass


The liquid scintillator density varies with temperature as:  
ρ_PC=((0.89179± 0.00003)-(8.015 ± 0.009)10^-4× T) g/cm^3,
where T is the temperature in degrees Celsius <cit.>.
The temperature is monitored at various positions inside the detector. The volume closest to the IV where temperature is recorded is the concentric Outer Buffer, where the thermal stability is measured to be better than 1^∘ C. 
In the FV, the maximum scintillator mass excursion corresponding to temperature variations is 0.1 ton, ∼0.1% of the FV mass. 
A Lomb-Scargle analysis (Sec. <ref>) on the temperature data was performed. The largest amplitude corresponded to a frequency of ∼0.6 year^-1, reflecting a significant real trend which anyhow cannot mimic  the annual modulation. 






§ MODULATION ANALYSIS




We have implemented three alternative analysis approaches to identify the seasonal modulation. The first is a simple fit to the data in the time domain (Sec. <ref>). The second is the Lomb-Scargle method (Sec. <ref>) <cit.>, an extension of the Fourier Transform approach. The third method is the Empirical Mode decomposition (EMD) (Sec. <ref>) <cit.>.

For each approach we define a set of time bins of equal length t_k and their
corresponding event rate R(t_k), obtained as the ratio of the number of selected events and the corrected life time (subtracted of the muon veto dead time and any down-time between consecutive runs).  



The time bins are too short to allow extracting a value of the ^7Be neutrino interaction rate via a spectral fit. 
We use the raw β-event rate instead, which include background contributions.

                                                                                                                         




 §.§ Fit to the Event Rate



Due to Earth's orbital eccentricity (ϵ = 0.0167), the total count rate is expected to vary as


    R(t) = R_0 + R[1 + ϵcos2 π/T(t- ϕ)]^2

where T is the period (one year), ϕ is the phase relative to the perihelion, R is the average neutrino interaction rate and R_0 is the time independent background rate.  This formalism is consistent with the MSW solution in which are no additional time modulations, at the ^7Be energies <cit.>.

In this approach, the event rate as a function of the time is fit with the function defined in equation <ref>.
Figure <ref> shows the folded, monthly event rate relative to the average rate measured in Borexino, with t=0,  365 representing perihelia. Data from the same months in successive years are added into the same bin.
Having normalized to 1 the overall mean value, the data are compared with Eq. <ref> and show good agreement with a yearly modulation with the expected amplitude and phase.
The no modulation hypothesis is excluded at  3.91 σ (99.99% C.L.) by comparing the χ^2 obtained with and without an annual periodicity.  



To extract the modulation parameters, we perform a χ^2 fit of the data with 30.43-day bins, without folding multiple years on top of each other.
Figure <ref> shows the event rate (in cpd/100 ton) along with the best fit.  
From <cit.>, the expected neutrino average rate in this energy range is ∼32  cpd/100 ton. 
The fit returns an average neutrino rate of R=33 ± 3 (cpd/100 ton), within 1σ of the expected one (χ^2/ndof = 0.68, ndof=42).
The best-fit eccentricity is ϵ = 0.0174± 0.0045, which corresponds to an amplitude of the modulation
of (7.1 ± 1.9)%, and the best-fit period is T=367± 10 days. Both values are in agreement with the expected values of 6.7% and of T=365.25 days. 
The fit returns a phase of ϕ=-18 ±24 days.  The robustness of the fit has been studied by varying the
bin size between 7 and 30 days, by shifting the energy range for selected events, and with and without α-β mlp inefficiency. 
Fit results are found not to vary greatly and are all in agreement with the expected modulation due to the Earth's orbit 
eccentricity. The resulting systematic uncertainty on the eccentricity is 10%.










 §.§ The Lomb-Scargle method



The second approach uses the Lomb-Scargle method. This extension of the Fourier Transform is well suited for our conditions since it can treat data sets that are not evenly distributed in time. 
In the Lomb-Scargle formalism, the Normalized Spectral Power Density, P(f), also known as the Lomb-Scargle periodogram and derived for N data points (R_1 … R_j … R_N) at specific times t_j, is 
evaluated and plotted for each frequency f as:


    P(f)=12σ^2{[ Σ_j (R_j - R) cosω
    (t_j-τ) ] ^2 Σ_j cos^2ω (t_j-τ).


    .+[Σ_j (R_j -R)sinω (t_j-τ)]^2Σ_j sin^2ω (t_j-τ)}





    R = R_1 + R_2 + R_3 + ...+R_NN = 1/N∑_j=1^N R_j 
    σ^2 = 1/N-1∑_j=1^N ( R_j - R)^2 
    tan 2 ωτ = ∑_j sin 2 ω t_j∑_j cos 2 ω t_j

where ω=2 π f.
After finding the frequency f_0 corresponding to the maximum of the Lomb-Scargle Power distribution <cit.>, the sine wave that best describes the time-series, in the case of 
a pure signal, is:


    R(t)=A cos ω_0 t + B sin ω_0 t


where, for ω_0=2 π f_0 and

    A=12σ^2[Σ_j R_jcosω_0 (t_j-τ)]^2Σ_j cos^2ω_0 (t_j-τ)


    B=12σ^2[Σ_j R_jsinω_0 (t_j-τ)]^2Σ_j sin^2ω_0 (t_j-τ)




The modulation amplitude is the peak-to-peak variation of the curve resulting from Eq. (<ref>).

For this analysis the data are grouped, after selection cuts, into 7-day bins as shown in Fig. <ref>.
The Spectral Power Density P(f) is calculated using the corresponding normalized event rate R(t_k)
 and it is shown in Fig.<ref>. 



The maximum of the periodogram is at f=1 year^-1 and corresponds to  a P(f)  value of 7.9. A zoom-in is shown in Fig. <ref>. 





 

 


Following  <cit.>, we have evaluated the significance of the largest peak found in the periodogram of our experimental data set with a toy Monte Carlo simulation assuming a realistic signal-to-background ratio and a  time interval of 4 years. 
Figure <ref> displays the P(f), at f = 1 year^-1, distribution (red filled area) obtained applying the Lomb-Scargle analysis to 10^4 simulations of a constant rate signal corresponding to the null
hypothesis (absence of modulation). This distribution is exponential as expected
for the power at a given frequency of the standard Lomb-Scargle periodogram of a pure white noise time series, Prob(P(f)>z)= e^-z <cit.>. 
In the plot, the vertical lines mark the 1σ (solid), 2σ (dashed) and 3σ (dotted) sensitivity to the null hypothesis. The blue distribution is obtained from 10^4 simulations of an expected yearly modulated signal plus constant backgrounds and its most probable value is P(f)=9.9 with rms of 4. 




The Spectral Power Density P(f) of 7.9 for f = 1 year^-1, obtained from the data, is within the range expected from Monte Carlo and corresponds to >3.5σ significance with respect to the null hypothesis. 

In addition we have estimated via Monte Carlo the significance of the two 4.5 high peaks
in the L-S periodogram.  Missing any a-priori information about the presence of periodicities
other than the annual one, the significance of these two peaks must be evaluated as
global significance, which takes into account the so called Look Elsewhere Effect,
i.e. the blind search over a frequency range <cit.>. Basically, one performs a Monte Carlo evaluation of the distribution of the highest peak induced by a pure noise time series over the
searched frequency interval. The significance (or p-value) is computed comparing the obtained distribution with the Power value of the highest peak detected in the Lomb-Scargle periodogram of the data. In this way  we determined for  the two 4.5 high peaks the p-value of 85%. Hence these two peaks are fully compatible with being pure noise induced fluctuations in the spectrum.

Finally, a sinusoidal function is constructed via Eq. (<ref>) for f_0=1 year^-1 and overlaid to the time-binned data in Fig. <ref> (red curve).
The peak-to peak amplitude is ∼ 5.7%, slightly less than that expected from the eccentricity of the Earth's orbit, because the Lomb-Scargle method cannot disentangle the background from neutrino signal.
The same analysis using data selected with slightly different cuts and without applying the rate correction for MLP inefficiency (see Sec. <ref>), returns consistent results. The resulting total uncertainty for the period is 4%, and for the amplitude 7%. No phase information is available with this technique.








 §.§ Empirical Mode Decomposition



The third method, the “Empirical Mode Decomposition” (EMD) <cit.>, has been
designed to work with non periodical signal, in order to extract the main parameters
from a time series as instantaneous frequency, phase and amplitude. 
The algorithm does not make any assumption about the functional form of the signal, in contrast to the 
Fourier analysis, and can therefore extract any time variation embedded in the data set.   

The EMD is a methodology developed to perform time-spectral analysis based on a empirical and iterative algorithm called
sifting, able to decompose an initial signal in a set of complete, but not orthogonal, oscillation mode 
functions called "Intrinsic Mode Function" or IMF <cit.>. 

Here we adopt a new technique for the noise assisted method called “Complete Ensemble Empirical Mode Decomposition with Adaptive Noise” 
(CEEMDAN) <cit.> showing a greater efficiency and stability on the final results than the EEMD method <cit.>.
The algorithm is more capable to separate the signals of interest from background because it removes the residual 
noise present in the final IMFs together with the spurious oscillation modes <cit.>.




  §.§.§ Standard Algorithm


The sifting algorithm (Sec. <ref>) requires a large number of points for a best performance. 
To maximize this number we chose bins of 1 day. As a consequence, statistical fluctuations dominate the dataset time-series
(red points in Fig. <ref>). However, the intrinsic dyadic filter <cit.>, 
removes all high frequency components created by the Poisson statistical noise.

The intrinsic mode functions, IMFs, are extracted from the original function through an iterative
procedure: the sifting algorithm. The basic idea is to interpolate at each step the local maxima and
minima of the initial signal, calculate the mean value of these interpolating functions, and
subtract it from the initial signal. The same procedure is then repeated on the residual subtracted signal
until suitable stopping criteria are satisfied. 

These are numerical conditions, which slightly differ in literature according to the approach followed (see e.g. <cit.>).
They aim at making sure that the IMFs obey two features inherited from harmonic functions: first, the number of extrema (local maxima
and minima) has to match the number of zero crossing points or differ from it at most by one;
second, the mean value of each IMF must be zero. 

The i-th IMF obtained by the k-iteration is given by:


    IMF_i(t) = x_i(t) - ∑_j=1^km_ij


where x_i(t) is the residual signal when all “i-1” IMF’s have been subtracted from the
original signal R(t), x_0(t) = R(t), and the m_ij are the average function of the max and min
envelopes at each j-th iteration. Following the results from a detailed simulation, 
we fixed the number of sifting iterations to 20. This number guarantees a good
symmetry of the IMF with respect to its mean value, preserving the dyadic-filter property of the
method (i.e., each IMF has an average frequency that is half of the previous one
<cit.>). Thus we obtain all i IMFs down to the last one called
“trend”, that is a monotonic IMF. 



The EMD approach features two potential issues: on one hand, the method is strongly dependent on small changes of
the initial conditions; on the other, mode mixtures could occur for a physical component present in the
data set especially when the ratio between signal and noise[In this case noise
means the statistical fluctuations of the rate with respect to the amplitude of the seasonal
modulation signal.] is low (about S/N=0.2, in our case). In order to account for these problems, a noise-assisted 
technique has been adopted. A random white noise signal (dithering) was added several times to
the data set under study and the average of all the IMFs taken.

As for the Borexino Phase-I analysis <cit.>, we repeat the single extraction of the IMF
1000 times, adding to the data a white noise component with an average value μ_wn = 0.0 and 
σ_wn = √(N_bin), where N_bin is the rate of the single bin (Poisson's error).
The main difference with respect to the Borexino Phase-I analysis is the use of the noise-assisted approach, 
called CEEMDAN. 

The final decomposition of our data set is shown in Fig. <ref>, where the lower frequency components identified by the algorithm become visible in the higher IMFs.
The ones shown are the resulting IMFs averaged over the 1000 extractions with different regenerations of 
white noise.

In particular, Fig. <ref>c shows the grey band corresponding to 1000 noise regenerated
IMF-7 containing the seasonal modulation. The resulting average
function is shown as black solid line, while the red-dashed curve corresponds to the expected seasonal modulation.



  §.§.§ Modulation Parameters Estimation

Here we can only provide a short account of the procedures to calculate the modulation parameters. A more
detailed and formal description of the numerical calculations and theoretical explanations are
reported in <cit.>.

The frequency and the amplitude values of a periodic function (as the seasonal modulation) are
constant in time. We therefore expect that in the IMF7 (Fig. <ref>) where a 
modulation of 1-year period is visible, these parameters will be constant in time, the average curve peaking on the 
expected values.
Naturally, due to the numerical procedure with which the “signal” has been obtained, some
small fluctuations of the frequency and of the amplitude are expected.



The IMF functions extracted by the sifting algorithm are not based on an analytical function.
Therefore, in order to extract information on frequency, phase and amplitude, it is necessary to build a
complex function z(t) by means of a Hilbert transform of the initial signal  <cit.>:


    z(t)=a(t)+ib(t)=A(t) e^- iθ(t),

in which the real part a(t) is the IMF and the imaginary part b(t) is the 
Hilbert transform of the real function:


    b(t)= 1/πP∫_t'a(t')/(t-t') dt'

where P is the Cauchy principal value. In Eq.(<ref>), A(t) is defined as 


    A(t)=√(a^2(t)+b^2(t)).

A(t) is also called the amplitude modulation function (AM), while


    θ(t)=arctan( b(t)/a(t))

defines the phase of the carrier function or frequency modulation (FM) function.
This method provides a function of the phase of the time that we can use to define the
instantaneous frequency (IF) as simple time derivative of the phase θ(t).
Unfortunately a direct calculation of the IF, starting from the signal, 
gives unphysical results with negative values for the frequencies. In order to solve this problem, 
an additional numerical procedure is required: the “Normalized Hilbert Transform” (NHT) 
<cit.>. 
Performing the NHT we obtain a normalized carrier function over all the time series.
Building the z(t) function, we are able to calculate a reliable instantaneous frequency 
function with a real physical meaning as follows:


    f(t)= dθ(t) dt.


We calculate the IF f(t) and the amplitude A(t) for all the IMFs extracted from each noise
regeneration and take the distribution of their average in time. A
Gaussian fit is applied to the resulting distribution to obtain f(t), A(t) and their respective errors.
  
In Fig. <ref> we compare IMFs obtained from the real dataset (Fig. <ref>c) with simulated data sets from a toy Monte Carlo with/without the sinusoidal signal expected for the seasonal modulation (Fig. <ref>a and <ref>b respectively).

For both real and MC data set, the resulting IMF average shows a very good agreement with the expected seasonal modulation function, while in the case of the null hypothesis (Fig. <ref>b) the 
amplitudes of the resulting IMFs are substantially smaller while frequencies and phases are varying randomly. 



A power spectrum is defined based on the average in time of the square amplitudes 
(⟨ A^2(t)⟩) (<ref>)  for each frequency ω (t).
Fig. <ref> shows the relative power spectra for the simulations with and without modulation
(Fig. <ref>a and  <ref>b) and the real data set, respectively (Fig. <ref>c).

The colored histograms are the Power Spectra from the last 4 IMFs, while the dark green 
are the full spectra of the whole set of IMFs (full dataset spectrum).

As expected in the presence of the seasonal modulation signal (Fig. <ref>a and <ref>c), we observe a
narrow peak centered on the expected frequency (f=1/T  = 2.73 × 10^-3 day^-1), while in the case of the null hypothesis this spectral component remains 
almost flat, featuring an amplitude comparable with other background IMFs that are present at higher 
frequencies. The power is an order of magnitude lower than the signal case (Fig. <ref>b).  



Applying equation <ref>, we compute the average parameters shown in Tab.<ref> for the simulated and real data. The results are in agreement with the expected seasonal modulation. 

                                                                                           


Based on the comparison of the power spectrum and the parameters resulting from the zero-modulation MC data sets we conclude the presence of a seasonal modulation.

We have calculated a χ^2-map varying both the phase and modulation amplitude of the sinusoidal function with respect to the average IMF obtained over the complete 1000 noise regenerations. The χ^2-contours are displayed in Fig. <ref>, where we assumed the standard deviation of the IMFs from the average curve to equal 1σ-uncertainties divided by the number of time bins minus one.






§ SUMMARY




Four years of Borexino Phase-II data have been analyzed searching for the expected annual modulation of the ^7Be solar neutrino interaction rate induced by the eccentricity of the Earth's orbit around the Sun.

Both the detector and the data have shown remarkable stability throughout the entire Phase-II period, allowing for the clear emergence of the annual periodicity of the signal.

Three analysis methods were employed: an analytical fit to event rate, a Lomb-Scargle periodogram and an Empirical Mode Decomposition analysis. Results obtained with all three methods are consistent with the presence of an annual modulation of the detected ^7Be solar neutrino interaction rate.  Amplitude and phase of the modulation are consistent with that expected from the eccentric revolution of the Earth around the Sun, proving the solar origin of the low energy neutrinos detected in Borexino.
The absence of an annual modulation is rejected with a 99.99% C.L.. The direct fit to the event rate yields an eccentricity of ϵ = (1.74 ± 0.45)%, while the Lomb-Scargle method identifies a clear spectral maximum at the period T=1 year. The EMD method provides a powerful and independent confirmation of these results. 
 



§ ACKNOWLEDGEMENTS


The Borexino program is made possible by funding from
INFN (Italy), NSF (USA), BMBF, DFG, HGF and MPG
(Germany), RFBR (Grants 16-02-01026 A, 15-02-02117 A, 16-29-13014 ofi-m,
17-02-00305 A) (Russia), and NCN Poland (Grant No. UMO-2013/10/E/ST2/00180). We acknowledge the generous hospitality and support of the Laboratory Nazionali del Gran Sasso (Italy). 




00


§ REFERENCES


bib:detector 
G. Alimonti et al. (Borexino Collaboration), The Borexino detector at the Laboratori Nazionali del Gran Sasso, Nucl. Instr. and Methods A 600 (2009) 568.
bib:be7 G. Bellini et al., (Borexino Collaboration), Precision Measurement of the 7Be Solar Neutrino Interaction Rate in Borexino, Phys. Rev. Lett. 107 (2011) 141302.
bib:pep Bellini, G. et al. (Borexino Collaboration),
First Evidence of pep Solar Neutrinos by Direct Detection in Borexino, Phys. Rev. Lett. 108, (2012) 051302.
bib:bxpp 
G. Bellini et al. (Borexino Collaboration), Neutrinos from the primary  proton-proton  fusion  process  in  the  Sun,  Nature  512 (2014) 383.
bib:b8 G. Bellini, et al. (Borexino Collaboration), Measurement of the solar 8B neutrino rate with a liquid scintillator target and 3 MeV energy threshold in the Borexino detector, Phys. Rev. D 82, (2010) 033006.
bib:carlos J. Bergstrm, M.C. Gonzalez-Garcia, M. Maltoni et al, Updated determination of the solar neutrino fluxes from solar neutrino data, J. High Energ. Phys. (2016) 132.
bib:geo  M. Agostini et al. (Borexino Collaboration),  Spectroscopy of geoneutrinos from 2056 days of Borexino data, Phys. Rev. D 92, 031101(2015) 031101.
bib:anti  G. Bellini, G. et al. (Borexino Collaboration), Study of solar and other unknown anti-neutrino fluxes with Borexino at LNGS, Phys. Lett. B 696, (2011) 191.
bib:muon G.  Bellini  et  al.  (Borexino  Collaboration),  Muon  and  Cosmogenic  Neutron  Detection  in  Borexino,  J.  Instrum.  6  (2011) P05005.
bib:serenelli16  A. Serenelli, Alive and well: a short review about standard solar models, 
arXiv: hep-ph/1601.07179 (2016)
bib:be2007 Bellini, G. et al. (Borexino Collaboration), First real time detection of 7Be solar neutrinos by Borexino, Phys. Lett. B 658  (2008) 101.
bib:long 
G.  Bellini  et  al.  (Borexino  Collaboration),  Final  results  of Borexino  Phase-I  on  low-energy  solar  neutrino  spectroscopy, Phys. Rev. D 89 (2014) 112007.

bib:SNO B. Aharmim (SNO Collaboration), Search for periodicities in the B8 solar neutrino flux measured by the Sudbury Neutrino Observatory, Phys. Rev. D 72 (2005) 052010.
bib:SK J. Hosaka et al. (Super-KKamiokande Collaboration), Solar neutrino measurements in Super-Kamiokande-I, Phys. Rev. D 73, (2006) 112001. 

bib:homestake R. Davis Jr.,  A review of measurements of the solar neutrino flux and their variation, Nucl. Phys. B - Proceedings Supplements, 48, (1996) 284.
bib:mils A.Milsztajin, A Search for Periodicity in the Super-Kamiokande Solar Neutrino Flux Data
arXiv:hep-ph/0301252 (2003).
bib:stu D.O. Caldwell and P.A. Sturrock, Further Evidence for Neutrino Flux Variability from Super-Kamiokande Data, arXiv:hep-ph/0305303 (2003)
bib:stu2 D.O. Caldwell and P.A. Sturrock, Evidence for solar neutrino flux variability and its implications, Astropart. Phys., 23, (2005) 543.
bib:SK1 J.Yoo et al. (SK Collaboration), Search for periodic modulations of the solar neutrino flux in Super-Kamiokande-I, Phys.Rev. D68, (2003) 092002.
bib:GNO L. Pandola, Search for time modulations in the Gallex/GNO solar neutrino data
Astropart. Phys. 22 (2004) 219.




bib:Lomb N.R. Lomb, Least-squares frequency analysis of unequally spaced data, Astrophysics and Space Science, 39 (1976) 447.
bib:Scargle J.D. Scargle, Studies in astronomical time series analysis. II - Statistical aspects of spectral analysis of unevenly spaced data, Astrophysical Journal 263 (1982) 835.

bib:Huang98 N. E. Huang, Z. Shen, S. R. Long et al., The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc. R. Soc. of Lond. A,  454 (1998) 903.

bib:MSW G. Bellini, G. et al. (Borexino Collaboration), Absence of a day-night asymmetry in the 7Be solar neutrino rate in Borexino, Phys. Lett. B 707, (2012) 22.



bib:horneJ. H. Horne and S. L. Baliunas, A prescription for period analysis of unevenly sampled time series, Astrophys. J. 302 (1986) 757.
bib:gio G. Ranucci., Likelihood scan of the Super-Kamiokande I time series data,
Phys. Rev. D 73  (2006) 103003.


bib:WuHuang2010 Z. Wu, N.E. Huang, On the filtering properties of the empirical mode decomposition, Advances in Adaptive Data Analysis, 2, No. 04, (2010) 397.
bib:Torres2011
M.E. Torres, M.A. Colominas, G. Schlotthauer, P. Flandrin, A complete ensemble empirical mode decomposition with adaptive noise, Proc. of 36th IEEE  Intern. Conference on Acoustic, Speech and Signal Processing, Prague, Czech Republic, (2011) 4144. 
bib:Colominas2014 M.A. Colominas, G. Schlotthauer G., M.E. Torres, Improved complete ensemble EMD: A suitable tool for biomedical signal processing, Biomedical Signal Processing and Control 14 (2014) 19.
bib:WuHuangEEMD2009
Z. Wu and N. E. Huang, Ensemble empirical mode decomposition: a noise-assisted data analysis method, Advances in Adaptive Data Analysis, 1, No.1 (2009) 1.


bib:WuHuang2004 
Z.Wu and N.E.Huang, A study of the characteristics of white noise using the empirical mode decomposition method, Proc. R. Soc. Lond. A 460 (2004) 1597.


bib:Rilling2003 G. Rilling, P. Flandrin, P. Goncalves, On empirical mode decomposition and its algorithms, Proc. of IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing NSIP-03, (2003) 177.

bib:onIF N.E.Huang et al., On instantaneous frequency,
Advanced in Adaptive Data Analysis,1, No.2, (2009) 177



 



