Solar neutrino measurements using the full data period of Super-Kamiokande-IV

,


I. INTRODUCTION
The existence of neutrino oscillations [1,2], which is a consequence of neutrino masses and mixing, is experimental evidence of elementary particle physics 'beyond the standard model'.Observations of solar neutrinos were first made by the Homestake experiment [3] using a radiochemical method, and then followed by real-time measurement with KAMIOKANDE-II [4] and other radiochemical experiments using gallium by SAGE and GALLEX/GNO [5][6][7].An initial indication of solar neutrino oscillations was obtained from the difference between the 8 B solar neutrino fluxes as measured in the elastic-scattering channel at Super-Kamiokande (SK) and the charged-current channel at the Sudbury Neutrino Observatory (SNO) in 2001 [8,9].Solar neutrino oscilla- tion was subsequently established by including neutralcurrent measurements from SNO [10].Solar neutrino oscillations were confirmed using reactor anti-neutrinos by KamLAND [11].Since these discoveries, Borexino and KamLAND experiments have measured the neutrino fluxes from different solar nuclear fusion processes, such as pp, pep, 7 Be, and carbon-nitrogen-oxygen (CNO) cycle [12][13][14][15].All measurements to date are naturally explained by neutrino flavor change due to neutrino oscillations with matter effects predicted by Mikheyev, Smirnov, and Wolfenstein [16,17], termed the MSW effect: higher energy neutrinos undergo adiabatic conversion from the electron flavor state to the second mass eigenstate.While neutrino oscillations and MSW effect is consistent with all current solar neutrino measurements, two distinctive predictions are yet to be observed: the characteristic energy dependence of the solar neutrino electron-flavor survival probability P ee (E ν ) distortion due to the MSW effect in the Sun and the day/night flux asymmetry induced by the matter effect in the Earth [18][19][20][21][22][23].One of the interests of solar neutrino experiments is to determine the neutrino oscillation parameters of ∆m 2  21 and sin 2 θ 12 .Independent of solar neutrino measurements, the KamLAND experiment used reactor anti-neutrinos to measure the same oscillation parameters, assuming CPT symmetry holds [24].Our previous paper on solar neutrino measurements [25] reported consistency in sin 2 θ 12 while observing about 2σ tension in ∆m 2  21 between the solar global fit and the KamLAND result.Further precise measurements of neutrino oscillation parameters are required to test the framework of three-neutrino oscillation as well as the conservation of CPT in the neutrino sector [26].
In this article, the results of solar neutrino analysis from the full observation period of the fourth phase of SK (SK-IV) are described.Moreover, the combined results together with the earlier phases of SK (SK-I, SK-II, and SK-III) are also presented.This article is organized as follows: Sec.II provides an overview of the SK detector and summarises the simulations.In Sec.III, improvements to the event reconstruction methods are explained.In Sec.IV, the calibration methods and the detector performance are presented.In Sec.V, the data set of SK-IV, data reduction, and its systematic uncertainties are described.In Sec.VI, observed solar neutrino results from SK are presented.In Sec.VII, VIII, and IX, the oscillation analysis from SK and other experiments are discussed.In the final section, we conclude this study and give future prospects.

II. SUPER-KAMIOKANDE-IV DETECTOR
The Super-Kamiokande detector consists of about 50, 000 tonnes of purified water in a stainless steel cylindrical water tank with photomultiplier tubes [27].The inner detector (ID) holds 32, 000 tonnes of water as physics target and its standard fiducial volume for solar neutrino analysis is 22, 500 tonnes.The outer detector (OD) is optically separated from the ID and used to veto cosmicray muons penetrating the mountain above [28].The ID is viewed by 20-inch photomultiplier tubes (PMTs) and the OD is viewed by 8-inch PMTs.SK was originally started in April 1996, and observation in phase IV was finished in May 2018 to refurbish the detector in preparation for dissolution of gadolinium into the water [29].Table I shows a summary of the experimental phases of SK.Starting with SK-II, an acrylic cover and fiber-reinforced plastic (FRP) case were installed around the ID PMTs to avoid a chain reaction of implosions [30].
The start of SK-IV corresponds to the installation of new 'QBEE' front-end electronics [31,32].The new electronics is capable of high-speed signal processing and records every hit of all the PMTs including hits from PMTs' dark current.
The coordinate of the SK detector is defined the same way as in the previous publication [25].The origin is at the center of the detectors, with X-and Y -axes lying in the horizontal plane and the Z-axis vertically upwards.We categorize the cylindrical ID surfaces into "top", "barrel", and "bottom" regions.

A. Water system and water quality
Water purification is extremely important for the experiment, not only to improve the attenuation length of Cherenkov photons but also to reduce the radioactive backgrounds due to dissolved isotopes.For these purposes, the SK-IV water is continuously recirculated through a water purification system and returned to the detector [27,33].The water quality, which affects the absorption and scattering of photons during their propagation, is monitored by using a laser system [28], and the water transparency (water attenuation length) is measured using cosmic-ray muon data.
The electrons (or positrons) from decays of stopped muons are used to monitor the variation over time of the water attenuation length.The light intensity observed by PMTs depends on the distance from the position of the decay electrons; light reaching more distant PMTs will be more attenuated.Using many events, the correlation between light observed and distance-to-PMT is used to measure the water transparency in the water tank [34].Figure 1 shows the time variation of water attenuation length over the data sample described in this article.Purified water is supplied from the bottom of the detector, and drained from the top of the detector.The temperature of the supply water is carefully controlled as this affects the amount of convection in the detector.Normally, convection inside the detector occurs below Z = −11 m [28].This water flow results in a small asymmetry in the water transparency as a function of vertical position Z, and is responsible for a top-bottom asymmetry (TBA).In order to model the Z-dependence of photon absorption, we introduce the parameter α TBA , which is defined as where N top , N bottom , and N barrel are the averages of the hit probabilities of top, bottom, and barrel of the ID.
The TBA parameter, α TBA , is measured by two calibration devices.One is the auto-Xenon lamp, which continuously injects light into the detector every 1 s.The other is a Ni-Cf calibration source, which is inserted into the detector to take calibration data about once per month, as shown in Fig. 2. Due to the water temperature dis-  The red dots (blue upward-pointing triangle) show the TBA value measured by the auto-Xenon calibration system (Ni-Cf calibration source) [28].The black vertical line shows when the temperature of the supply water was changed, as described in the main text.Bottom: The difference between two TBA measurements, which is normally within ±0.5% during the SK-IV data set.
tribution in the detector, a TBA of about −5% has been observed by the auto-Xenon lamp, indicating a higher hit rate in the bottom region.The TBA measured by the Ni-Cf calibration source has a similar value to that measured by the auto-Xenon lamp.The difference over the whole data set is around ±0.5% level.
In February 2018, we changed the temperature of the supply water injected at the bottom of the ID region, from +13.06 • C to +13.52 • C.This lowered the density of the water in the bottom region and evoked large scale convection.The water temperature measured at various places in the detector became consistent after about two weeks of re-circulation, this indicated that water is fully mixed at that time and hence water quality should be uniform across the entire detector volume.After achieving this uniform water condition, the TBA value became significantly smaller, but a residual 1-2% asymmetry remained.This remaining asymmetry is interpreted as the asymmetry of photon detection efficiency of the PMTs, and is corrected for in the energy reconstruction as described in Sec.III B. Another convection test was done in May 2018.A consequence of these 'water convection tests' was increased Rn contamination (see Sec. V A 1).

B. Radioactive contamination in the water
A large part of the low-energy intrinsic background of the SK-IV detector comes from radon (hereafter Rn) remaining in the purified water.To prevent Rn gas from the mine entering into the water tank, the detector is tightly sealed and a low-Rn buffer gas is continuously supplied to the air layer between the top of the tank and the water surface.To keep the detector clean, the radon concentrations in the buffer gas, as well as the experimental site of the SK detector, are monitored by several kinds of Rn detectors [35,36].The measured Rn concentrations in the buffer gas supply is 0.08 ± 0.07 mBq/m 3 while that in the air layer is 28.8 ± 1.7 mBq/m 3 .This indicates Rn sources exist inside the detector.
In addition to the buffer gas monitoring system, another Rn measurement system was developed in 2013 [33] which evaluates the Rn concentration in purified water in the detector directly.Since then, water has been sampled from the various positions in the detector and its Rn concentrations have been measured continuously.The Rn concentrations in the center and bottom regions of the SK-IV detector in 2015 are < 0.23 mBq/m 3 and 2.63±0.22mBq/m 3 , respectively.Details of the Rn study will be given in a future publication [37].

C. Time variation of PMT gain and dark rate
During SK-IV, we observed a time variation of the gain and dark rate of PMTs, as shown in Fig. 3 and Fig. 4.
The average PMT gain increased by around +10%∼ + 15% during SK-IV, depending on the production year of PMTs.Since the hardware threshold at which a PMT's signal is recorded was fixed during the time depending on the PMT production year, as shown in Fig. 4. In our previous report, based on the first 1664 days of SK-IV [25], the effect from these variations was small and within the energy scale systematic uncertainty.However, since the full SK-IV is substantially longer than earlier SK phases, the effect cannot be ignored in the full SK-IV period.A correction factor is therefore included in the SK-IV solar neutrino energy reconstruction to account for this time variation.Some improvement in gain over time is often seen in large PMTs.We are not sure why our PMTs continue to improve after such a long period of operation, but we believe that our empirical correction factor will account for it within our stated systematic uncertainties.This is described in more detail in Sec.III B and Appendix A.

D. SK-IV detector simulation
A Monte Carlo (MC) simulation that reproduces physics events is used to understand the SK-IV detector performance.This simulation is based on Geant3 package [38], and customized for the SK-IV detector.
In order to reproduce the data more accurately, the following parameters are evaluated, typically day-by-day, then used in the solar neutrino event simulations in SK-IV: water transparency, TBA of the water transparency, PMT dark rate, and PMT gain.For the 8 B solar neutrino event generation, the Winter spectrum [39] is used for the initial 8 B solar neutrino energy.For the hep solar neutrinos, the Bahcall spectrum [40] is used.We consider the actual detector live time when simulating the direction of solar neutrinos.Ejection of a recoil electron by neutrino-electron elastic scattering is simulated using the cross section including radiative correction [41].The expected 8 B and hep solar neutrino event rates in the whole SK ID volume (32.5 kt) without neutrino oscillation are 294.7 events/day and 0.6375 events/day, respectively.Flux values used in this calculation are based on the neutral current results from SNO [42] for 8 B solar neutrino ((5.25 ± 0.20) × 10 6 cm −2 s −1 ), and BP2004 [43] standard solar model (SSM) for hep solar neutrino (7.88 × 10 3 cm −2 s −1 ).We generate about 10 5 simulation events/day each for 8 B and hep solar neutrino fluxes.

A. Vertex and Direction
The event vertex and direction are reconstructed using the same method as the previous analysis [25].The electrons scattered by the solar neutrinos can travel only a few cm in water and thus the location of the Cherenkov emission is recognized as a point.Under this assumption, the vertex position is reconstructed with a maximum likelihood fit of the photon arrival time on the PMTs in the solar neutrino analysis [44].The log-likelihood function is defined as, where x is the vertex position, being tested, N hit is the number of PMTs that have received light (hereafter hit-PMT), t i is the hit time of i-th hit-PMT, P (∆t) is the probability density function of the timing residual (t i − t tof − t 0 = ∆t) for a single photo-electron signal, t tof is the timing after subtracting the time-offlight (hereafter TOF), and t 0 is the initial time of the event.The probability density function for the timing residual is extracted from LINAC calibration data [45], and is shown in Fig. 6.Although small secondary peaks corresponding to the PMTs' after pulses can be seen in the residual density function, the likelihood function is strongly peaked at ∆t = 0. Figure 7 shows the vertex resolution of an electron event in the standard fiducial volume.Here, the vertex resolution is defined as the root mean square of the residual difference between the generated position and the reconstructed vertex position.To evaluate this, electrons are generated uniformly across the whole ID volume with a random direction.Then events reconstructed within the standard fiducial volume are selected and the difference between the vertex positions is calculated.Typical vertex resolutions are 101.5 cm, 64.0 cm, 49.3 cm, and 39.1 cm at 3.49 MeV, 6.49 MeV, 9.49 MeV, and 14.49 MeV of electron kinetic energy, respectively.
The direction of the event is reconstructed with a maximum likelihood method to find a Cherenkov ring that best matches the position of the hit PMTs.The loglikelihood function is defined as where N 30 is the number of hit-PMTs with residual time (t i −t tof ) within 30 ns; θ dir,i is the angle between the electron direction and the vector from the reconstructed vertex to the position of i-th hit-PMT; E is the reconstructed energy of the electron and f (cos θ dir,i , E) is the probability of simulated electron events with θ dir,i and E.
The logarithm is weighted by an acceptance term based on θ i , which is the angle between the vector from the reconstructed vertex to the i-th hit-PMT position and the direction that the i-th hit-PMT is facing.The parameter a(θ i ) is the typical acceptance of a PMT, which is shown in Fig. 8.The use of the acrylic cover and FRP cases since SK-II has little effect on the acceptance of photons at small angles to the tube facing, but the decrease of acceptance at large angles is noticeable, so this effect is corrected in this solar analysis.Figure 9 shows typical distributions of the function of f (cos θ dir,i , E).Because of the multiple scattering of electrons in water, the dis- tribution peaks at 42 • , but is quite broad, with the width depending on the electron energy.The 1σ directional resolution is defined as the angle that includes 68.3% of events in the distribution of the angular difference between their reconstructed direction and their true direction.It is estimated from simulated electron events that are generated uniformly in the whole ID volume with random direction, and reconstructed in the standard fiducial volume.Figure 10  vertex and direction reconstruction methods and their performances can be found elsewhere [25,46].

B. Energy
Electron energy is reconstructed based on the number of hit PMTs within 50 ns of the expected photon arrival time from the reconstructed vertex.Then, an effective number of hits, N eff , is calculated by applying a correction for the relative difference of PMT performance, contributions of dark rate and late arrival hits, the effective surface area of PMTs, and light attenuation in the water.Finally, N eff is converted to a recoil electron energy using a function derived from simulations of mono-energetic electrons in the standard fiducial volume in the SK-IV detector.
The method of calculating N eff was improved since the previous analysis [25] and is expressed as where N 50 is the number of hit PMTs included in the 50 ns window, i is an index running over these hit PMTs and N all (N alive ) is the number of total (live) ID PMTs at the time the event was recorded.Descriptions of the other variables in the above N eff definition are following.
X i : Hit occupancy correction calculated for the i-th hit-PMT and its eight surrounding PMTs.Assuming the Poisson distribution, the correction is calculated as where x i is the fraction of hit-PMTs in the 3 × 3 PMT block surrounding the i-th hit-PMT.The details are explained in Appendix B in Ref. [47].
ε tail : Correction for late hits due to reflected and scattered light.It is calculated as where N 100 is the number of hit-PMTs in a 100 ns window, R dark is the average dark rate overall live PMTs, and t 50 is the interval of 50 ns.Since the latter 50 ns is dominated by scattered light and background hits, it is not appropriate to apply the corrections in Eq. ( 4).So, we only correct statistically for the extra light outside the first 50 ns.
ε i dark : Correction for dark hit contribution.This parameter is calculated using the measured dark rates for the corresponding PMT as, where r i dark is measured dark rate for i-th hit-PMTs, and t 50 is the interval of 50 ns.
C i : Correction for PMT gain increase.This correction for the observed data was empirically determined and expressed as where G i is relative gain of i-th hit-PMT with respect to the gain in October 2008.Those relative gains are evaluated for five PMT groups separated by their production time, as shown in Fig. 3.A similar correction is also applied for the simulated data, since the relative gain increase is also included in the SK-IV detector simulation in this analysis.
S(θ i , φ i ): Correction for PMT coverage due to the fact that each PMT has a spherical surface, as a function for polar angle (θ i ) and azimuthal angle (φ i ) of incident photons in the local coordinate of each PMT.This correction was evaluated with a Geant4 [48][49][50] based simulation of the detector geometry.This correction is improved from the previous analyses by introducing more realistic geometry to the simulation.Details of improvements to this and the attenuation (P i ) correction are given in Appendix A.
P i : Correction for light attenuation between the vertex and the i-th hit-PMT.In the previous analysis, a uniform light attenuation length was assumed, but this analysis accounts for variation of the optical properties of the water.
QE i : Correction for relative detection efficiency for the i-th hit-PMT.This factor includes all effects leading to a light sensitivity variation between individual PMTs.This correction was re-evaluated using special data taken with a Ni-Cf source in February 2018.During this data taking, water convection in the SK-IV detector was intentionally evoked to realize uniform water quality across the detector volume.
Several changes are made to this analysis to improve energy reconstruction.The correction for PMT gain increase (C i ) is newly introduced for this analysis and improves the stability of the energy scale to be within ±0.5% over more than nine years of operation.The corrections for PMT geometrical acceptance (S(θ i , φ i )) and light attenuation (P i ) significantly improve the uniformity of the energy scale.The correction for the dark rate contribution, ε i dark , was also improved by incorporating PMT-by-PMT variation (in the previous analyses, an averaged value for all the PMTs was used.) Figure 11 shows non-uniformity of N eff for simulated electron events.The standard deviation of relative position dependence of N eff was improved from 1.7% to 0.5%.More details of the improved energy reconstruction method are given in Appendix A. The LINAC calibration system consists of an electron gun, a linear accelerator, beam pipes, collimators, magnets, and a beam trigger.The details of the system are described elsewhere [45].The last part of the beam pipe is inserted vertically into the SK detector and a single electron is injected at a time into the detector.The LINAC can inject electrons directly into the detector, which mimics the electrons produced by elastic scattering interactions of the 8 B or hep solar neutrinos.
Mono-energetic single electrons are injected into the SK-IV detector at several positions.The energy of injected electron ranges from 4.4 MeV to 18.9 MeV in total energy.
Since the previous report [25], LINAC calibrations were conducted in in 2016 and 2017.We analyzed calibration data taken in SK-IV with the improved event reconstruction algorithm as described in Sec III.For the LINAC data analysis, a LINAC trigger signal from the trigger counter around the end of the beam pipe is required.To reject multiple electron events, separation by hit timing distribution after subtracting TOF (defined in Sec.III A) is applied.Since the timing of electrons occasionally overlaps, we rejected events whose energy is beyond 3σ from the mean of the energy distribution.

Determination of the absolute energy scale in SK-IV
In this analysis, we use good-quality LINAC calibration data taken in 2010, 2012, 2016, and 2017 while the calibration data taken in 2009 is not used because of the lower quality of the water transparency at that time.The calibration data was taken at nine different positions in the SK detector.By comparing the peak of effective hits distribution (N eff in Eq. ( 4)) from these calibration data and simulated events, the absolute energy scale of the MC simulation code for SK-IV was tuned to match the LINAC data.In particular, the absolute correction factor for PMT quantum efficiency is determined by this analysis.Figure 12 shows the difference of the N eff between the calibration data and the MC simulation after the tuning.
To evaluate the position dependence of the energy scale during the SK-IV period, the weighted mean is obtained from the difference of N eff between the data and MC simulation, where the injecting beam energies are 6.989, 8.861, 13.644, and 18.938 MeV. Figure 13 shows the difference of data and MC simulation after taking the weighted mean for each position.The remaining position-dependent energy scale uncertainty in SK-IV is calculated by taking the volume average for nine positions and the result is a difference of 0.40%.Therefore, ±0.40% is used as the systematic uncertainty on the energy scale.
After the MC simulation is tuned, it can be used to derive the conversion from N eff to electron energy.For this purpose, mono-energetic electrons are generated with random vertex position and direction in the whole ID volume, and then the reconstructed energy distribution of events reconstructed in the standard fiducial volume is examined.Figure 14 shows the reconstructed energy distribution of the mono-energetic electron simulations.Each distribution of reconstructed energy is fitted with a Gaussian function and the peak energy and the deviation are obtained.The fit results are used to create the conversion function, and determine the energy resolution function of the detector, which is found to be: where E is the reconstructed electron total energy in MeV.This is comparable to the energy resolution evaluated in the past publications: σ(E) = −0.123+ 0.376 √ E +0.0349E in SK-III [46], and σ(E) = −0.0839+0.349 √ E + 0.0397E in SK-IV 1664 days [25].From here on, the energy scale tuned with this method will be used to estimate various performance metrics of the SK-IV detector.

Energy resolution measurement
After determining the energy-scale function, the energy resolution of the LINAC calibration data and MC are compared.Figure 15 shows the typical energy distributions of the LINAC calibration.The energy distributions are fitted with a Gaussian function and its peak value and deviation are obtained.Then, the energy resolution is calculated by dividing the deviation by the peak value.Figure 16 shows the systematic uncertainties of the energy resolution that comes from the difference between the data and MC.As energy is higher, the uncertainty also increases.We used a polynomial function to describe the energy resolution as a function of the energy.We estimated energy resolution systematic uncertainty in SK-IV is at most 3% level in the range of the solar neutrinos.This difference is taken into account in the analysis as an energy resolution systematic uncertainty.

Angular resolution measurement
Although the LINAC calibration is mainly used to determine the absolute energy scale of the detector, the electron beam makes it possible to evaluate the angular resolution for electrons in the water, since the direction of the LINAC beam is known.Figure 17 shows typical opening angle distributions between the direction of the beam injection and the reconstructed direction.Because of the multiple scattering of an electron in water, the distribution of the opening angle spreads depending on the electron's energy.
Using the same definition of directional resolution as Sec.III A, the angular resolution of an electron is estimated.Figure 18 compares the angular resolution obtained from LINAC data and from MC simulation.This difference between the data and MC simulation in angu-

B. DT calibration
The LINAC uses a beam pipe inserted vertically through the calibration ports on the top of the SK de-tector.As such, electrons can only be injected in the downward vertical direction and so the LINAC system cannot be used to test the directional dependence of the energy scale.Instead, the DT (deuterium-tritium neutron generator) calibration source [51] is used to monitor the directional stability of the energy scale of the SK-IV detector.In addition, the DT device is designed to be portable and easy to operate and this allows us to calibrate the energy scale at many more positions than the LINAC calibration.While the LINAC is used at 9 positions, the DT calibrations are taken at 35 different positions.
The DT device generates neutrons via the reaction of 2 H+ 3 H → 4 He+n.The neutrons are captured on 16 O in the water and 16 N is created via (n, α) reaction [52].The decay of 16 N, whose Q-value is 10.4 MeV, mainly emits an electron with an energy of 4.3 MeV and a single γray with energy of 6.1 MeV.These particles are emitted isotropically.This allows probing directional systematic uncertainty in the SK energy reconstruction.
In the actual calibration procedure, the DT device is lifted about 100 cm from the target position by the crane of the SK detector soon after its neutron generation.Although the 16 N emits an electron and γ-ray isotropically, the upward events are affected by the shadow due to the DT device above the target position.In order to reduce such shadowing effect when calculating the  average N eff , we removed the event whose direction is 0.8 < cos θ DZ < 1.0, where the θ DZ is defined as the zenith angle respect to the vertical (Z) axis of the detector (hereafter, we call this angle as the detector zenith angle).Since the half life of 16 N decay is 7.13 s, we collected the data for 40 s after the lifting of the DT device was completed.The former 20 s is used for the signal region and the latter 20 s is used for the background subtraction.

Reconstructed recoil electron kinetic energy [MeV]
Figure 19 shows the typical distribution of the reconstructed energy and the directional dependence of the effective hits obtained from the DT calibration data by setting the calibration device at the central position of the SK detector together with the MC simulation.In this analysis, the detector azimuthal angle is defined as the angle in the X-Y plane.The differences are sufficiently small for all azimuthal angles and cosine of the detector zenith angles.
In order to check the stability of the energy scale throughout the SK-IV data set, we evaluate the time dependence of the energy scale by comparing the DT calibration data and the MC simulation.In this analysis, a position-weighted average is performed to obtain the one value from the deviation measured at various positions in the detector.For this purpose, we assigned a geometrical weight to each position where DT calibration data are taken.That weight is based on the fraction of the volume nearest each calibration position contributes to the total 22.5 kt fiducial volume.There is about −0.5% difference between the DT calibration data and the MC simulation, hereafter referred to as the offset.Therefore, the DT calibration is used to evaluate relative variations from the offset, not absolute energy scale calibrations.
At the moment, we do not know the origin of the offset.
Figure 20 shows the time variation of N eff determined by the DT calibration.The bottom panel of Fig. 20 also shows the difference between the data and the MC simulation.There may be a slight increase over the datataking period, which probably originates from the modeling of the water transparency and the top-bottom asymmetry in the MC simulation described in Sec.II A. Although the difference between data and simulation fluctuates at the ±1% level, the difference reduced to zero during the period of the convection study on February 2018 when the water quality throughout the detector is expected to be uniform, supporting the idea that the residual variation is related to non-uniformity of the water transparency.
Because the DT calibration data were taken at various positions in the SK detector, we can also evaluate the position dependence of the energy scale by comparing each calibration position.Figure 21 shows the variation of the energy scale with height and radius-squared.Both of them are consistent within the ±0.4% fluctuation from the offset, which is the systematic uncertainty of the position dependence estimated from the LINAC calibrations in Sec.IV A.
Figure 22 shows the directional dependence of the energy scale on the detector azimuthal angle and detector zenith angle.For the azimuthal angle, the directional systematic uncertainty is estimated by the variation from the offset.As seen in the top panel of Fig. 22, the fluctuation of the azimuthal angle dependence is less than ±0.1% level throughout the detector.
For the detector zenith angle, the systematic uncertainty on the directional dependence of the energy scale is estimated using the same way as the azimuthal angle.However, the N eff in the most upward bin (0.8 < cos θ DZ < 1.0) is noticeably reduced, an effect that was also seen in the previous analysis [25].This is attributed to the increased shadowing of PMTs from the DT generator.This is not implemented in the current DT simulation, so this bin is excluded when setting the anisotropy systematic.We therefore assign a systematic uncertainty of ±0.10% on the energy scale from the directional dependence.

C. Summary of the systematic uncertainty of the energy scale in SK-IV
Table II shows a summary of the systematic uncertainties in the energy scale determination in SK-IV.The position dependence of the energy scale is estimated to be ±0.40% by taking the volume average of the difference of the effective hits in data and MC at different LINAC injection positions, as shown in Fig. 13.The accuracy of the LINAC calibration itself is set by the beam energy determination, which uses a Ge detector and is estimated  to be accurate to ±0.21% [45].Extrapolating this calibration to other time periods also has an uncertainty because of variations in the water transparency over time as shown in Fig. 1.These are continuously monitored with the decay electrons from stopping cosmic-ray muons, but the precision is limited by the number of stopping muons recorded during the reference period of LINAC operation.The energy scale uncertainty from this is estimated as ±0.11%.Finally, the directional dependence of the energy scale is evaluated by using the DT calibrations and the resulting uncertainty is ±0.10% as described in Sec.IV B. In total, we estimate an uncertainty of ±0.48% on the absolute energy scale in SK-IV.

V. DATA ANALYSIS A. Event selection
Events used in the analysis must pass a series of selection criteria.A basic explanation of each selection step is given here, with more details to be found in the previous publication [25].Using the full data sample of SK-IV, the selection criteria are optimized to maximize the significance S/ √ S + B, where S and B are defined as the number of signal and background events after the selection cut, respectively.

Run selection for solar neutrino analysis
The typical run length in SK-IV is one day.For this solar neutrino analysis, a series of run selections are applied to select good-quality data.As in previous analyses, short runs (of less than five minutes), calibration runs, runs with hardware and/or software troubles, and runs with strange OD rates or strange numbers of bad PMT channels are rejected.After the good run selection, the total live time for the solar analysis in SK-IV is 2970 days, from September 2008 to May 2018 as listed in Table I.This corresponds to an extra 1306 days more than the previous publication (1664 days).Since May 2015, the experiment has run with a lower software trigger threshold as explained in the next section.In February and May 2018, convection studies were performed to obtain better estimates of systematic uncertainties; during these periods events below 5.49 MeV are rejected since a large number of background events originating from the Rn in the detector are observed [37].Table III summarises the data sets of solar neutrino analysis in SK-IV.

Trigger scheme
The online software trigger in SK-IV is based on the number of coincident PMT signals.A lower threshold is desirable for the solar neutrino analysis, but causes the data rate to rise sharply, so events must pass a software trigger based on the number of hit-PMT within a 200 ns window.In the case of the solar neutrino events, the energy of recoil electrons is approximately proportional to the total number of the hit-PMTs since most of the hit-PMTs receive one photon.
In May 2015, the software trigger threshold was changed from 34 hit-PMTs to 31 hit-PMTs [31,53].Figure 23 shows the time variation of the estimated trigger efficiencies near the trigger threshold using the MC simulation.For events above 4.99 MeV the trigger efficiency is basically 100%.For events between 3.49 and 3.99 MeV, there was a noticeable inefficiency due to the 34-hit threshold, this is significantly reduced by moving to the lower 31-hit threshold.The trigger efficiency is highly correlated with the gain of PMTs (Fig. 3) and water transparency (Fig. 1).After changing the threshold of data acquisition (vertical light-blue dashed line), the efficiency of 3.49-3.99MeV is improved significantly.

1st reduction
The first step of the data reduction is the elimination of obvious background events and noise originating from hardware in the detector.The method of selection is the same as in the previous analysis.
• Since the PMTs are operated at high voltage, an arc discharge can occasionally occur between the dynodes during data taking.When this happens, other PMTs also register optical signals.If the maximum charge of a single PMT is larger than 50 p.e. and there are more than three hits on the surrounding 24 PMTs, the event is identified as a 'flasher' event, and removed.
• Events that occur less than 50 µs after a previous event are rejected to remove decay electrons from stopping cosmic-ray muons, and instrumental noise from PMT after pulses.
• Several calibration sources remain inside the detector during data taking.Events that are triggered by external calibration triggers and scheduled calibration events are obvious and rejected.In addition, the calibration sources inevitably have higher radioactivity than the water they displace.Therefore events below 4.99 MeV are excluded if their vertex position is less than 200 cm from the calibration sources or 100 cm from the cables supporting the sources (all cables run parallel to the Z-axis from the top of the detector to the source position).

Spallation cut
Cosmic-ray muons pass through SK at ∼2 Hz, and will occasionally shower; some of these showers produce hadrons such as pions or neutrons (hadronic showers) [54][55][56].The showers can interact with or break apart 16 O nuclei, generating radioactive isotopes (spallation).Those isotopes that undergo βs and/or γ decays have similar reconstructed energies as solar 8 B and hep neutrino interactions.Their half-lives range from milliseconds to tens of seconds, with the most abundantly produced isotope, 16 N, having a half-life of 7.3 s.This "spallation background" dominates all others above 5.5 MeV and motivates a number of selection cuts to reduce its impact.
Previous solar neutrino analyses of SK data reduced this background using three variables calculated and checked against any cosmic-ray muons in the previous 100 s: the closest distance between the muon track and the candidate; the time difference between the muon and the candidate and the excess light from the muon compared to that expected from minimum ionization.A likelihood was formed from these three variables, and a cut was placed on this likelihood.The cut was tuned to have 20% loss of solar neutrinos (effectively the same as a dead time), and was 90% effective at removing spallation.The details of the previous method can be found in Ref. [25,34].In this analysis, the "spallation cut" is complemented and improved.The improvements are summarised in this section, and full details will be provided in the recent publication [57].Firstly, a completely new method directly tags hadronic showers responsible for isotope production.Such showers typically produce a large number of neutrons, and in pure water, these neutrons are captured on hydrogen via n + 1 H → 2 H + γ (2.2 MeV), and this can be used to tag neutrons.While the total light produced by a 2.2 MeV γ is small (typically only seven detected photo-electrons) and therefore difficult to see, the large neutron multiplicity allows this method to remove 54% of spallation-induced events with only a 1.4% loss of signal efficiency.The updated cut depends on the location of the solar neutrino candidate relative to a cluster of tagged neutrons, the multiplicity of the neutron cluster, and the time between the candidate and the muon responsible for the shower.Due to the difficulty of triggering on such a small signal, the data acquired through a special Wideband Intelligent Trigger (WIT) [58,59] must be used; its trigger efficiency for a 2.2 MeV γ is about 13%.The WIT was only available for 388 days, towards the end of the SK-IV period.
Secondly, a new method uses the spallation decays to veto themselves.Events within the standard fiducial volume that pass most event quality cuts, and above 4.99 MeV are used to create a sample of "veto events".Solar neutrino candidates within 4 m and 60 s of the veto events are removed.This removes 47% of spallation events with an effective dead time of 1.3%.
Thirdly, the existing spallation cut was updated and re-tuned to more effectively tag spallation remaining after the two new cuts were applied.The updates included: a better muon fitter; the removal of artificial saturation of PMTs; updated probability density functions (PDFs) for the three input variables and the addition of a fourth variable.This new input variable utilizes the distance between the solar neutrino candidate and the peak of the muon energy loss along the track.This peak typically corresponds to the position of electromagnetic showers.The resulting likelihood cut value was tuned separately for cases with and without neutron cluster data, to achieve the same (90%) reduction in spallation-induced events as used previously.Figure 24 shows a comparison of the likelihood of the signal and background of the non-neutron data period.
The improvement is quantified in the reduction of The multiple spallation cut is already applied, so only 82% of remaining spallation needs to be removed to achieve 90% overall spallation removal effectiveness.
the effective dead time: 8.6% (10.5%) for the dataset with (without) neutron clusters.The total SK-IV effective dead time from these cuts is 10.2%.This leads to an increase of 12% in the number of solar neutrino signal events (∼7000 events) in the 2970-day sample.Figure 25 shows the difference in the number of events accepted using the updated spallation cuts, which includes a large contribution to the solar neutrino signal sample.

Summary of the event selection
After the spallation cut, we applied the same event selections as the previous analysis [25], with cut points re-tuned, to obtain the final data sample.
Figure 26 shows the event rate of the SK-IV 2970-day final data sample.The event rate in 3. 49 energy ranges are most affected; event rates in higher energy ranges are relatively stable throughout the data set.
Figure 27 shows the energy distribution after each reduction step of the SK-IV 2970-day data set.Background events above 5.5 MeV are mostly removed by the spallation cut.On the other hand, background events below 5.5 MeV are mostly removed by the ambient cut, external cut, and tight fiducial volume cut [25].
Figure 28 shows the signal efficiency after each reduction step as a function of reconstructed recoil electron kinetic energy.The improvement of the signal efficiency above about 6 MeV is due to the new spallation cuts, as described in Sec.V A 4. In the lower energy regions, tighter cuts are applied to remove more background events, which was optimized by the cut point tuning in this analysis.

B. Signal extraction
The observed number of solar neutrino signal events are extracted from the final data sample using an extended maximum likelihood function fit.The fit is made to the angle, θ Sun , between the reconstructed event direction and the Sun at the time of detection.This distribution is used since neutrino-electron elastic scattering is peaked in the forward direction, so the solar neutrino signal corresponds to a peak on top of a uniform background.This extraction method is similar to previous analysis [25] but in the SK-IV 2970-day analysis the binning with multiple scattering goodness (MSG) parameter, defined in [25], is newly considered.The likelihood function is defined as where N energy is the total number of the energy bins, FIG.27.Reconstructed electron kinetic energy distribution of the observed data after each reduction step.The square, triangle, rhombus, cross, and circle correspond to after 1st reduction, spallation cut, ambient cut, external cut, and tight fiducial volume cut, respectively.The dotted lines correspond to the SK-IV 1664-day final data sample in the previous analysis [25] and can be compared to the circle markers.N MSGi is the total number of the MSG bins, n ij is the total number of events in the {ij} bin, b ij (s ij ) are the background (signal) probability density function, Y ij is the fraction of the total signal events (S) in the {ij} bin, B ij are free parameters corresponding to the number of background events in the {ij} bin, S is a free parameter corresponding to the total number of solar neutrino events in all energy and MSG bins.The fitting parameters, S and B ij , represent the number of signal and background events, respectively.These parameters are obtained by maximizing the likelihood function.The uncertainty of this extraction method is evaluated with dummy solar angle distributions produced with solar MC events.The difference between generated and extracted solar neutrino events is taken into account as the signal extraction systematic uncertainty.Details are explained in Appendix B.

C. Systematic uncertainty
The uncertainties on solar neutrino measurements are estimated using simulation programs and calibration data.For the efficiency of the reduction steps, LINAC data and MC events are compared.For the vertex shift, Ni-Cf calibration data and MC are used in the same way as previous analyses.
Table IV shows a summary of the systematic uncertainty on the observed solar neutrino event rate in SK across the whole energy range.The total systematic uncertainty on the solar neutrino flux in SK-IV is estimated as ±1.4% in the 3.49-19.49MeV region.Table V shows TABLE IV.Summary of the systematic uncertainty on the observed solar neutrino event rate for SK-IV and compared to SK-I [34], SK-II [30], and SK-III [46]. the uncertainty in each energy range in SK-IV.These are treated as energy-uncorrelated uncertainties in the oscillation analysis.Energy-correlated uncertainties, like energy scale, resolution, and expected solar neutrino spectrum shape, are estimated with the simulation by artificially shifting the relevant parameters.To evaluate the systematic uncertainties related to the energy scale determination, the energy scale in the simulation program is artificially changed according to the estimations in Sec.IV C. Using this modified simulation program, the solar neutrino events are generated and then applied same reduction process.Finally, we investigate the difference of the extracted solar neutrino signal events.For the energy resolution, the same treatment is done with the uncertainty estimated in Sec.IV A 3. In addition, the theoretical uncertainty of the 8 B spectrum shape is also considered [41].These uncertainties are correlated in the solar neutrino energy spectrum.Figure 29 IV.

VI. FLUX AND ENERGY SPECTRUM MEASUREMENTS
In this section, the results using the SK-IV 2970-day data set as well as the combined results with other phases are described.The basic conditions for the solar neutrino simulation events used in this analysis are given in Sec.II D. The analysis in this section does not take into account the effects of neutrino oscillations.For this reason, it is explicitly indicated as "MC(Unoscillated)" in figures.

A. Total 8 B neutrino flux
The direction of a recoil electron from elastic scattering is strongly correlated to that of the incident solar neutrino.Figure 30 shows the solar angle distribution of the final data sample from 3.49 to 19.49 MeV obtained in this analysis.The number of extracted solar neutrino events is 65, 443 +390 −388 (stat.)± 925 (syst.).This number corresponds to a 8 B solar neutrino flux of Φ8 B,SK-IV = (2.314± 0.014 ± 0.040) ×10 6 cm −2 s −1 (11) assuming a pure electron neutrino flavor component, without neutrino oscillation.The ratio of the extracted events to that expected by assuming SNO's neutral current flux [42] is 0.441 ± 0.003 (stat.)± 0.007 (syst.).The 8 B solar neutrino flux measured in SK-IV is therefore consistent with those from previous phases within their total uncertainties.Figure 31 shows a summary of the measured flux ratios in SK.
The combined flux ratio from all the SK phases is At lower neutrino energies, the MSW effect in the Sun is weaker and neutrino flavor change reverts to a vacuumlike mechanism.For solar neutrinos this transition is expected to occur around 3 MeV, so SK's lowest-energy events are of particular interest to confirm our under-  [60][61][62][63], Borexino [64], KamLAND [65], and SNO+ (pure water) [66].The left (right) vertical axis shows the measured flux without neutrino oscillation (the ratio of the measured flux to the SNO's neutral current flux [42]).
standing of the MSW effect.Figure 33  signal-to-noise ratio is small because of the contamination from the background events, the peak of the solar neutrino signal is clearly observed over the background rate.The number of extracted events in the energy range of 3.49-3.99MeV is The statistical significance is about 10 sigma.This is achieved by reducing the background events, mainly from Rn daughters, in the central region of the SK-IV detector [33].

B. Yearly flux
The solar activity cycle is a roughly 11-year periodic change of sunspot numbers and reconfiguration of the magnetic field at the surface of the Sun.Although the standard solar model predicts the production rate of solar neutrinos is constant on this timescale, it does not consider the periodical activities of the Sun, such as the rotation inside the Sun or the variation of the sunspot numbers.The combined SK data sample contains the period from 1996 to 2018 and this long-term observation covers nearly two solar activity cycles, cycle 23 and 24, so it can be used to check for any correlation.
Figure 34 shows the yearly averaged fluxes observed in the different phases of the SK detector together with the corresponding sunspot numbers.To test the correla- tion between the observed yearly fluxes and the sunspot numbers, the χ 2 between the average flux value and the yearly flux values is defined as: where the summations are over p, the SK phase (I to IV); n p is the number of operating years in each phase; and t is year within a SK phase.The combined SK average flux (with no annual variation) is r, while d p,t is the observed yearly flux in year t of SK phase p.The statistical uncertainty on each d p,t is σ p,t .Systematic variations of SK phases are described by the nuisance parameters α p , where τ p is the systematic uncertainty for each SK phase p.Using the observed data in the energy range from 4.49 MeV (exceptions: 6.49 MeV for SK-II and 5.49 MeV for the 2018 data point) to 19.49 MeV, the minimum χ 2 for a steady flux (parameter to be varied) is calculated with the total experimental uncertainties as: where N d.o.f (= 22) is the degree of freedom for χ 21 .This corresponds to a probability of 58.9%.The solar neutrino rate measurements in SK are fully consistent with a constant solar neutrino flux emitted by the Sun.

C. Flux Measurement for Day and Night
Due to flavor-specific interactions of solar neutrinos with the Earth's matter, the solar neutrino interaction rate measured by SK depends on the solar zenith angle (θ z,solar ), which is the angle between the vertical (Z) direction of the SK detector and the neutrino direction from the Sun when an event occurs (i.e. the time of day).In most cases, the Earth matter effects lead to a "regeneration" of electron flavor, i.e. the electron flavor survival probability P ee is larger during the night compared to the day.The apparent day-time flux and the night-time flux of solar neutrinos in SK are measured separately in the SK-IV 2970-day data set to test this Earth matter effect: Φ day 8 B,SK-IV = (2.284± 0.020 ± 0.032) ×10 6 cm −2 s −1 , Φ night 8 B,SK-IV = (2.341± 0.019 ± 0.033) The live time of the day (night) is 1434 days (1536 days).The day/night asymmetry parameter is then calculated as = −0.025± 0.012 (stat.)± 0.014 (syst.).( 17) In addition, we fit the amplitude of the expected zenith angle variation to the observed data.Here, we use the convention that during the day-time cos θ z,solar ≤ 0 and during night-time cos θ z,solar > 0. The day/night asymmetry parameter extracted from the best-fit amplitude to the solar zenith angle variation (see Sec. VIII) is A SK-IV, fit D/N = −0.0262± 0.0107 (stat.)± 0.0030 (syst.).( 18) , this parameter depends on the oscillation parameters.The expected asymmetry depends also on the oscillation parameters as well as the energy range in which the flux is measured.With current oscillation parameters, the day/night asymmetry is expected to be at the few-percent level above the MeV region while no asymmetry is expected in the keV region.In a previous publication [68], an indication for a non-zero day/night flux asymmetry was found in SK; the best fit was at the few-percent level.Borexino also measured the day/night flux asymmetry of solar 7 Be neutrinos [69] and the difference is consistent with zero.The SNO collaboration [42] measured a day/night asymmetry consistent with zero, and consistent with the SK result.All three results favor the MSW-LMA solutions.Figure 35 shows the solar zenith angle distribution of the observed flux divided into five-day and six-night bins.

D. Spectrum results
Figure 36 shows the observed signal spectrum as well as the expected spectrum calculated assuming the SNO NC flux and the hep flux of the standard solar model [43,63].In the calculation of the expected spectrum, no neutrino oscillation is assumed and the signal efficiency in Fig. 28 is corrected.To test the contribution of hep neutrinos, the expected spectrum without hep neutrinos is also shown.Although the total hep neutrinos flux is 3 orders of magnitude smaller than the 8 B neutrinos, its contribution can be seen in the the highest energy bins of the recoil electron spectrum since the end-point of the hep neutrino energy spectrum (18.8 MeV) is slightly higher than that of 8 B neutrinos (16 MeV) [70].
Figure 37 shows the energy spectrum of solar neutrinos measured in the SK-IV 2970-day data set.The top panel  panel of Fig. 38.

VII. OSCILLATION ANALYSIS A. SK-IV Analysis
The solar neutrino oscillation analysis is based on the same methods as [25] which were first described in Ref. [34,71].The survival probability is calculated in a two-neutrino framework: P 2ν ee (A mat , θ 12 , ∆m 2 21 ), where θ 12 and ∆m 2  21 are the vacuum oscillation parameters, and A mat is the potential caused by non-zero electron density when the neutrino travels in matter.
To account for three-neutrino effects-specifically a non-zero value of θ 13 -this is modified, following [72]: In fitting neutrino oscillations, we use a constraint on θ 13 derived from reactor neutrino experiments [73]: Since solar neutrino measurements are not sensitive to changes in sin 2 θ 13 less than 0.005, the closest calculation point to this value (0.020) is effectively used.The SK analysis constrains neutrino flavor oscillation by measuring the rate of elastic scattering of 8 B and hep neutrinos with electrons, the spectrum of the recoiling electrons, and the time variation of the interaction rate.Such variations of rate occur due to matter effects on neutrino oscillations in the Earth, and therefore depend on the time of day.The most stringent constraints are obtained when combining all SK operating phases, however, SK-IV dominates in the combined fit: SK-IV observed the largest number of solar neutrinos, it has the lowest energy threshold, and the smallest systematic uncertainties.An unbinned likelihood L is used to fit the number of solar neutrino interactions to the angular distribution (see Eq. ( 10)).The likelihood is modified to account for oscillation-induced spectral distortions and rate time variations and depends then on oscillation parameters, neutrino fluxes, and nuisance parameters describing energy-correlated systematic uncertainties due to the detector energy scale, energy resolution as well as the neutrino energy spectrum: where r i (cos θ z,solar ) is the expected solar neutrino event rate as a function of cos θ z,solar in the i-th energy bin, and r ave i shows the average of the event rate over all the solar zenith angle bins.The spectral distortions enter this likelihood via the factors Y ij describing the expected fraction of signal events in the bin labeled with i and j.
A time-independent likelihood L ave is defined by replacing the neutrino interaction rates r i (cos θ z,solar ) with the time averages r ave i , so the zenith angle dependence factor r i (cos θ z,solar )/r ave i is set to 1.The large number of events renders the maximization of this likelihood (which depends on several nuisance parameters) computationally difficult.Therefore, the log-likelihood is formally split into a time-dependent log-likelihood ratio and a time-independent log-likelihood: The first two terms are identified with the time-variation log-likelihood ratio log L time = log L − log L ave .(log of the ratio of time variation over no time variation).Of course, by definition, L ave depends only on the interaction rate and recoil electron spectrum, so to save computation time, it is maximized by minimizing an equivalent χ 2 .This spectral χ 2 is defined as χ 2 spec = −2 log L ave and is then approximated by a binned sum of squared differences of expected and observed event rates divided by the uncertainties (see below in Eq. 24) augmented with several nuisance parameters describing neutrino fluxes and energy-correlation uncertainties.Those nuisance parameters don't affect L time as much as χ 2 spec , so the minimization can be factorized: Minimizing χ 2 spec determines the nuisance parameters, and once best-fit values are obtained L time can evaluated using just the best-fit values from χ 2 spec .This calculation method for the time variation derived from the matter effects in the earth is the same as described in [34,71].The total equivalent χ 2 of SK is then where χ 2 time = −2 log L time .The spectrum part is where β and η are, respectively, the dimensionless scaling parameter of the 8 B and hep neutrino fluxes, and are the subject of the minimization.The expectations b i (h i ) are the ratios of the oscillated 8 B (hep) neutrino interaction rates over the unoscillated combined neutrino interaction rate in MC in energy bin i. d i are the observed ratios with uncertainties σ i .The energy-correlated spectral distortion factor f i (τ, ǫ, ρ) is controlled by three constrained, dimensionless nuisance parameters: τ is the systematic deviation of the 8 B neutrino spectrum (common to all phases); ǫ is the systematic energy scale deviation and ρ is the systematic energy resolution deviation.An extra term can be added to represent external constraints on the neutrino fluxes.The constraint on β comes from the measurement of the NC reaction of the SNO experiment, (5.25 ± 0.20) × 10 6 cm −2 s −1 [42], and the hep constraint η of (7.88 ± 15.76) × 10 3 cm −2 s −1 is based on the standard solar model [43] while the uncertainty roughly corresponds to the limit reported by SNO [74].

B. SK I/II/III/IV Combined Analysis
Figure 40 compares the SK-IV spectrum to the spectra from previous phases.The energy scale and resolution of each phase are individually calibrated and evaluated.Within the uncertainties of this calibration, all four phases agree with each other.Not only these energycorrelated but also the energy-uncorrelated systematic uncertainties are evaluated for each phase separately.To combine all phases of SK, the spectral χ 2 is assembled from the spectral χ 2 from each of the four phases of the experiment.The χ 2 is defined as in [76].The time variation likelihood ratios log L time,p of all four phases p are just added.When combining SK data from all phases, strong oscillation parameter constraints are obtained from the recoil electron spectrum shape and day/night variation alone without external constraint on the elastic scattering rate (but requiring consistency within uncertainties of the different phases of SK).This is done by removing the external constraint on β (based on the NC interaction rate with deuterium measured by SNO) in Φ.
Figure 41 compares the allowed regions at 1, 2, and  3σ confidence level contours obtained from spectrum and measured day/night rate variation with those using the absolute elastic scattering rate in addition.Even without flux constraint, the no oscillations scenario is excluded by the SK combined analysis at greater than 3σ.Both sets of contours are very consistent, but the addition of the absolute rate of course greatly improves the oscillation constraints: allowed areas are smaller and closed intervals are obtained even at 4σ and 5σ.Those allowed areas are exclusively in the "Large Mixing Angle" region (sin 2 (θ 12 ) ≈ 0.3 and ∆m 2 21 ≈ 6 • 10 −5 eV 2 ) Using 5805 days of SK data containing more than 100,000 solar neutrino interactions, SK measures the oscillation parameters to be sin 2 θ 12,SK = 0.324 +0.027 −0.023 and ∆m 2  21,SK = (6.10+1.26 −0.86 × 10 −5 ) eV 2 (Without the total rate: sin 2 θ 12 = 0.33 +0.19  −0.13 and ∆m 2 21 = (6.1 +1.6 −1.8 × 10 −5 ) eV 2 ). Figure 42 shows, on a linear scale, the allowed regions by combining all data taken during the four SK experimental phases.When compared to the Kam-LAND measurements, there is a slight tension (about 90% C.L. for a ∆χ 2 ≈ 2.3) in ∆m 2  21 as shown in Fig. 42.A significant discrepancy in the oscillation parameters obtained from neutrinos (SK solar neutrino data) versus anti-neutrinos (KamLAND reactor anti-neutrino data) would imply CPT violation.The SK and KamLAND data are consistent with CPT conservation within their estimated uncertainties. Figure 43 compares the measured with the best-fit MSW spectra.Both the solar best-fit value (see section VII D) for ∆m 2  21 as well as the KamLAND measurement lead to good agreement with the SK spectral data.Note that the different SK phases apply different distortions due to uncertainty in energy scale and resolution.Therefore, a statistical average spectrum can only used for illustration purposes since the predictions of each phase differ, in particular near the 8 B endpoint.Figure 44 shows this statistical average together with the MSW predictions.

C. Super-Kamiokande and SNO combined analysis
The neutrino oscillation analysis method to combine data with other experiments is similar to previous SK publications [25,34,71].Since SNO data also contains solar 8 B neutrino interactions in a very similar energy range, the combined analysis of SK and SNO plays the most important role in the global analysis.The SNO collaboration has performed a fit of all SNO data to six parameters (Φ B , c i , and a i [42] ) describing neutrino oscillations in the SNO experiment: three of these parameters are quadratic coefficients of the day-time electronflavor survival probability (developed around 10 MeV: c 0 , c 1 , and c 2 ), two more parameters describe a linear approximation of the day/night asymmetry of this survival probability as a function of energy (also developed around 10 MeV: a 0 and a 1 ).The sixth parameter Φ B is the total 8 B neutrino flux of all (active) flavors.SNO published the best-fit values of the six parameters, the fit errors, and the correlation matrix.In our anal- ysis, these six parameters are mapped to s i , and their errors and correlation matrix are assembled in the covariance matrix V SNO .For a given set of oscillation parameter θ 12 , θ 13 , ∆m 2  21 , the day-time survival probability and the day/night asymmetry of the survival probability is calculated as a function of energy.The day-time probability is fitted with a quadratic function, and the day/night asymmetry with a linear function.The fit uses the published relative sensitivity of SNO as a function of energy.From these fit parameters, an expected set of is added to the SK χ 2 , then the SK+SNO χ 2 is defined as In this equation, the hep flux is constrained in Φ to (7.9 ± 1.2) × 10 3 cm −2 s −1 to be consistent with SNO's assumption in obtaining the SNO parameters while the β constraint of Φ is removed and β is identified with the SNO parameter Φ B .Then β, η, ǫ p , ρ p , and τ are minimized.The definition of parameters is the same as in [76].Figure 45 compares the oscillation constraints of SK data with those of SNO data: SNO data constrains the mixing angle θ 12 better than SK data, while SK data constrains ∆m 2  21 better than SNO data.Figure 45 also shows the SK+SNO combined fit, which significantly improves the θ 12 constraint compared to the one from just SNO data, while the improvement of ∆m 2  21 compared to just using SK data is small.The combined fit also eliminates alternate small mixing angle allowed regions (appearing in the SNO only analysis at 4σ) and small ∆m 2 21 allowed regions (appearing in the SNO only analysis at 1σ and in the SK only analysis at 5σ).This SK+SNO analysis gives the most stringent constraints on θ 12 and ∆m 2   21   with neutrinos independent of solar model predictions of the solar neutrino fluxes: sin 2 θ 12,SK−SNO = 0.305±0.014and ∆m 2  21,SK−SNO = (6.10+1.04 −0.75 × 10 −5 ) eV 2 .

D. Global oscillation analysis
In addition to SNO, other solar neutrino data from Borexino [77], the Homestake experiment [78], Gallex/GNO [79], and SAGE [80] are considered.Unlike in previous publications, Borexino data and radio-chemical solar neutrino data (Homestake, GALLEX/GNO and SAGE) are now treated separately from each other.To analyze Borexino data we extract from Ref. [77] pp, pep, and 7 Be neutrino fluxes, errors and correlations.For the correlations we assume a general Gaussian covariance matrix and match it to the published MC pair correlations of fit parameters in Figure 5 of Ref. [77]; Figure 46 demonstrates this match.These fit parameters describe neutrino interactions and radioactive background rates: Borexino calculates sin 2 (Θ 13 )=0.0218±0.0007sin 2 (Θ 12 )=0.305±0.014∆m 2  21 =(6.10 spectra for each neutrino and background species, and (scales by the corresponding parameter) they are fit to the observed Borexino data spectrum.After matching the pair correlations, we then reduce the Gaussian covariance matrix to a smaller matrix that describes just pp, pep, and 7 Be neutrino interaction rates by marginalization: the neutrino interaction rate vector is then r Borexino = r pp r pep r7 Be = 134.0 2. 43 48.3 counts (100 tonne) −1 day −1 and the inverse covariance matrix is The corresponding inferred neutrino fluxes from Borexino data are φ Borexino = φ pp φ pep φ7 Be = 6.10 0.0127 0.499 These fluxes already take into account the neutrino oscillations using sin 2 θ 12 = 0.306, sin 2 θ 13 = 0.02166, and ∆m 2 21 = 7.50 × 10 −5 eV 2 [77].To calculate neutrino flux predictions for Borexino data for different oscillation parameters, we integrate the neutrino-electron elastic scattering cross section in the region where that particular neutrino species gives the largest contribution to the total Borexino spectrum.We form ratios of the predicted rate with a particular oscillation parameter set over the predicted rate without oscillations.In the case of the "standard parameters" (meaning the ones assumed by Borexino in [77]), these ratios are 0.689, 0.614, and 0.638 for pp, pep, and 7 Be, respectively.We then use solar model predictions with varying oscillation parameters to predict the Borexino event rates and fit those predictions to the observed rates using the Borexino covariance matrix V rate as well as the solar model neutrino flux vector r SSM and covariance matrix V SSM (converted to Borexino interaction rates).The components of the flux vector r are the solar neutrino species pp, pep, 7 Be, 8 B, 13 N, 15 O, 17 F, and hep.The relative solar model covariance matrix (i.e.constraining the fraction of the neutrino flux deviations) in units of 10 −3 (so that the first fraction is 1.0 For V rate and r Borexino , 0 values are padded for neutrino species that are not pp, pep or 7 Be.
We then fit the Borexino spectrum to the neutrino fluxes constrained by the solar model with The minimum χ 2 is then which is Borexino's contribution to the global χ 2 .We convert the rate vector and covariance matrix back to the neutrino flux vector and covariance matrix.The SSM+Borexino neutrino flux vector and covariance are modified by the SK+SNO determination of the 8 B and hep neutrino fluxes (and covariance), so that there is no impact of either the SSM 8 B or hep neutrino flux uncertainty on the analysis.The radio-chemical covariance matrix V RC is obtained from the flux covariance matrix V via the cross sections (and errors) of the target isotopes of the radio-chemical experiments (for some details about the covariance method see Ref. [81]).The radio-chemical rate measurements are then fit via where R obs and R exp are the observed and expected signal rate, respectively.The indices n and m run over Chlorine and Gallium (Gallex/GNO and SAGE are combined into one R obs ).The χ 2 of the solar global fit is defined as To extract the best values of θ 12 and ∆m 2 21 we combine solar experimental neutrino data and KamLAND reactor anti-neutrino data [24] (solar+KamLAND) by simple addition of the χ 2 functions, assuming no correlation of the solar and KamLAND results: ) Both Eq. ( 28) and Eq. ( 29) are evaluated with an external θ 13 constraint of sin 2 θ 13 = 0.0218 ± 0.0007.Figure 47 shows the oscillation parameters allowed by SK and SNO data, all solar data, KamLAND data, and all solar+KamLAND data.
The best-fit oscillation parameters from all solar experiments and KamLAND are sin 2 θ 12,global = 0.307 ± 0.012, ∆m 2  21,global = (7.50+0.19 −0.18 ) × 10 −5 eV 2 .The oscillation results of this analysis show a tension of ∆m 2  21 between the neutrino and anti-neutrino of about 1.5σ as shown in the right panel in Fig. 47, and it is slightly stronger for the global solar analysis compared to the SK+SNO analysis.
The global solar neutrino analysis has some sensitivity to θ 13 independently from reactor anti-neutrino measurements since the high-energy solar neutrino branches ( 8 B and hep) undergo MSW flavor conversion while the lowenergy solar neutrinos (pp, 7 Be, and pep) change flavor by averaged vacuum oscillations.Figure 48 shows the resulting contours of the mixing angles θ 12 and θ 13 .The best-fit value from all solar data of sin 2 θ 13 = 0.032 +0.021 −0.022 , is statistically consistent with zero as well as the reactor anti-neutrino measurements.When combining all solar data with KamLAND data (whose anti-neutrinos are subject to non-averaged vacuum oscillations), the bestfit value is almost unchanged, but the preference for a non-zero value gets somewhat stronger: sin 2 θ 13 = 0.030

VIII. DAY/NIGHT ASYMMETRY AMPLITUDE FIT
Studying the earth matter effect by making separate solar neutrino rate measurements during the day and the night results in comparatively large systematic uncertainties.Those arise from detector asymmetries (the directionality of solar neutrino recoil electrons means that different parts of the detector are illuminated during the day compared to the night) and directional energy scale effects as well as systematics due to the angular distribution of background events.Such systematic errors can be reduced by fitting the amplitude of the rate variation and extracting the corresponding day/night asymmetry from the fit to the amplitude.This fit is done using the likelihood of Eq. ( 21).The rate variation (calculated from standard oscillation parameters) r(cos θ z,solar ) in each bin is scaled by an amplitude factor α such that the corresponding day/night asymmetry scales with α, but the average rate is unchanged [71].Each energy bin rate variation is scaled by the same constant α.The corresponding day/night asymmetry for the SK-IV data is A SK-IV, fit D/N = −0.0262± 0.0107 (stat.)± 0.0030 (syst.).(30) This fit value is based on the expected zenith variation shape from ∆m 2 21 = 6.1 × 10 −5 eV 2 with an expected asymmetry of −0.0238.Figure 49 shows the dependence of the fit on ∆m 2 21 .In the region of interest, the fitted amplitude agrees well with the expected amplitude.The significance of a non-zero day/night asymmetry is 2.4σ.The systematic uncertainty includes energy scale, energy resolution, event selection, the density of electrons from the Earth model [82], and background angular distributions.Table VI summarizes the systematic uncertainties for the day/night amplitude fit.The first four contributions are based on Ref. [68], where the dominant last contribution was reevaluated.The amplitude method was improved in three important ways: Firstly, below 7.49 MeV, the data was split into three MSG regions which improves the stability to fluctuations of intrinsic radioactive background.Secondly, data between 3.49 and 4.49 MeV was added, resulting in a combined fit to 39 data samples (8 × 3 samples below 7.49 MeV, and 12 samples between 7.49 and 13.49MeV, 3 samples between 13.49 and 19.49MeV).
Lastly, the shapes of the angular distribution for the background PDFs vary as a function of solar zenith angle in five regions for day-time events, five regions for event candidate times where solar neutrinos pass only through the Earth's mantle, and another region for times where solar neutrinos pass also through the Earth's core.The amplitude fit allows free variation of the number of background events in each of the eleven regions as well as each of the 39 energy bins (429 different background numbers are fit).To estimate the systematic error due to the uncertainty in the angular background PDFs, we use data to measure the background distribution in detector coordinates (zenith and azimuth of the reconstructed direction).The uncertainty of that measurement is propagated to the calculation of the background PDF (as a function of cos θ Sun ), and, subsequently, the amplitude fit result.Due to the larger event statistics as well as these improvements, this important systematic uncertainty was reduced from 0.006 to 0.0027.Figure 50 shows the energy dependence of the asymmetry from amplitude fit using the SK-IV 2970-day data set.The energy region between 3.49 and 4.49 MeV differs by 2σ from the average.However, exclusion of this energy region does not significantly change the best-fit result (A SK-IV, fit D/N = −0.0241± 0.0109 for the higher energy threshold).D/N = −0.0286± 0.0085 (stat.)± 0.0032 (syst.).(31) Here, the asymmetry parameter is expressed based on the SK-I energy range, so the expected asymmetry is a bit stronger: −0.0242.Since this result differs from zero by 3.2σ, we find evidence for the existence of earth matter effects on solar neutrino oscillation.The asymmetry parameter depends on the expected zenith angle variation shapes and therefore on the oscillation parameters.The dependence on the mixing angle is negligible.The dependence on ∆m 2  21 is somewhat stronger.Figure 52 shows that dependence by analyzing the data taken throughout four SK phases.In particular, at the solar+KamLAND best-fit value of ∆m 2 21 = 7.5 × 10 −5 eV 2 , the combined SK day/night amplitude fit corresponds to a slightly smaller asymmetry of A SK, fit D/N = −0.0274± 0.0083 (stat) ± 0.0032 (syst.)where an asymmetry of −0.0172 is expected.Zero asymmetry differs from this measurement by 3.1σ.For reference, at the previously favoured [68] ∆m 2 = 4.8 × 10 −5 eV 2 we obtain A SK, fit D/N = −0.0273± 0.0086 (stat.)± 0.0032 (syst.)where an asymmetry of −0.0384 is expected.Zero asymmetry differs from this measurement by 3.0σ.

IX. ENERGY SPECTRUM ANALYSIS
The MSW effect [16,17] of the higher energy 8 B neutrinos is a unique feature of solar neutrino flavor physics as well as an important test of standard weak interaction theory.It predicts almost complete adiabatic conversion of the electron flavor state at neutrino production in the core of the Sun to the second mass eigenstate when neutrinos leave the Sun.Therefore, the electron flavor survival probability becomes just sin 2 θ 12 (in the two neutrino approximation), more or less independent of neutrino energy provided it is far above the "resonance energy" at the solar core.Lower energy 8 B neutrinos (as well as pp, pep, and 7 Be solar neutrinos) on the other hand undergo regular vacuum neutrino oscillations which average out to an energy-independent electron flavor survival probability P ee of 1 − 1 2 sin 2 (2θ 12 ).This results in an energy-dependence of that survival probability from a lower value of ≈ 0.3 at high energy (> 10 MeV) to a higher value of ≈ 0.6 at low energy (< 1 MeV) called the "upturn".Each phase of Super-Kamiokande measures the ratio of observed neutrino-electron elastic scattering over the no-flavor change expectation as a function of recoil electron energy (see Figures 40 and 43).The upturn of P ee leads to a distortion of the measured ratio.In order to test the presence of such an upturn, P ee (E ν ) is parameterized [25,42] in several different ways: Using the P ee,par (E ν ) (par is either quad, cubic, or exp), the modified expected energy spectrum of the recoil electron b i,par and h i,par (corresponding b i and h i in Sec.VII A) are made for each experimental phase, considering the energy resolution function (Eq.( 9)) and the average day/night asymmetry.Then, the spectral χ 2 is calculated by using the b i,par and h i,par .Finally, the fitted function giving the smallest chi-square (defined as χ 2 min ) becomes the best fit.Among ∆χ 2 = χ 2 − χ 2 min less than 1, the regions between the maximum and the minimum P ee,par (E ν ) at each E ν are defined as 1σ region.
Figure 53 shows the electron neutrino survival probabilities as a function of neutrino energy, obtained using data from all the four SK run periods.The thick blue(medium gray) line is the Pee(Eν) distribution expected from neutrino oscillation at all solar + KamLAND best-fit parameter set, and the thick green(light gray) line is the Pee(Eν) distribution expected from neutrino oscillation at all solar best-fit parameter set.
As shown in Table VII, both χ 2 values are well-approximated by cubic function χ 2 values of 66.34 ("sol.+KL")and 65.76 ("solar").The best-fit energy-independent ("indep") P ee = 0.336 has χ 2 = 67.20,so the solar (solar+KamLAND) best-fit equivalent cubic function is favored by 1.2σ (0.9σ) over an energy-independent P ee .The quadratic approximation is also reasonable: χ 2 = 66.34 ("sol.+KL")and χ 2 = 65.63 ("solar").Of course, the best-fit energyindependent P ee = 0.336 is the same, so the solar (so-lar+KamLAND) best-fit equivalent quadratic function is favored by 1.3σ (0.9σ).The exponential approximation is similar: χ 2 = 66.30("sol.+KL")and χ 2 = 65.71("solar") with similar conclusions.In summary, the SK spectrum measurement favors the existance of an "upturn" by 1.2σ.Table VIII shows a summary of the fitting results of the coefficients.The SNO constraint on the "upturn" is obtained and combined with the SK results as follows: from the SNO parameters s i , the first three are calculated from the quadratic fit parameters c 0 , c 1 , and c 2 of P ee,quad while with day/night asymmetry fitting parameters a 0 and a 1 are set to the oscillation best fit.The SK+SNO χ 2 is then the same as in Eq. (26). Figure 54 shows the electron neutrino survival probability distributions as a function of neutrino energy obtained from Eq. ( 32) with SK and SNO data.The SK+SNO combined result on quadratic P ee coefficients favor a distorted spectrum by 2.1σ. Figure 55 shows the SK+SNO combined result in the context of other solar neutrino survival probability measurements (assuming the standard solar model predictions of the unoscillated neutrino fluxes).The SK+SNO result fits in well with the other data as well as the MSW+neutrino oscillation prediction.

X. CONCLUSION
The fourth phase of the Super-Kamiokande (SK) has measured the solar neutrinos from September 2008 to May 2018.The total operation period of SK covers almost two solar activity cycles of 23 and 24.
In order to obtain precise solar neutrino measurements, several improvements to the analysis are applied in SK-IV.The data acquisition threshold was lowered  The solar neutrino rate measurements in SK are fully consistent with a constant solar neutrino flux emitted by the Sun.
The SK-IV time variation data fit results in a day/night asymmetry of A SK−IV,fit D/N = −0.0262± 0.0107(stat.)± 0.0030(syst), while a fit to all SK data gives A SK,fit D/N = −0.0286± 0.0085(stat.)± 0.0032(syst).This is a 3.2σ direct evidence for the existence of earth matter effects on solar neutrino oscillation.The fit assumes ∆m 2 21 = 6.1 × 10 −5 eV 2 , but similar conclusions hold over the region of interest.
SK-IV data measures the solar oscillation parameters to be sin 2 θ 12,SK−IV = 0.308 +0.030 −0.029 ∆m 2  21,SK−IV = 6.9  2  21 between all solar experiments and KamLAND persists at about 1.5σ.Further precise measurements may shed light on this tension in the future.
As described in Sec.II C, a constant upward drift of the PMT gain was observed.This effectively decreased the hit detection threshold at the front-end electronics, QBEEs [31], and affected the global energy scale as well as the position dependence.This effect was evaluated using decay electrons from stopping muons inside the detector.An empirical correction factor, 1/(1 + 0.226G i ) was applied, where G i is the relative gain of the i-th PMT.This minimizes the time variation of the light yield from decay-electron events after correcting for light attenuation.Figure 56 shows the relative size of the effective number of hits for PMTs, grouped by their production time before and after correction.This demonstrates that this method effectively corrects the observed time dependence of the number of hits for each PMT group.This correction of the time variation was cross-checked with an independent set of calibration data taken with the Ni-Cf source.We confirmed that the new correction of the gain drift successfully removes the time dependence of the energy scale.

Correction for effective PMT coverage
Because of the curved structure of the PMT sphere, the effective surface area of the PMTs varies significantly depending on the photon incident angle.The PMT coverage is approximately 40% for normal incident photons and increases as the incident angle gets shallower.The correction of this effect is implemented into the calculation of N eff as the (S(0, 0)/S(θ, φ)) term.In the previous analyses, this effect was estimated by a simple simulation of PMT geometrical structure shape arranged on a flat surface.In this analysis, this was improved with a Geant4 based simulation that includes the light propagation in the acrylic cover and the PMT glass structures.In addition, the PMT arrangement on the curved structure was also considered for the correction for the barrel PMTs. Figure 57 shows the S(θ, φ) values used for this analysis.

Correction for light attenuation in water
Correction for light attenuation in water was also improved with a more accurate representation of the water properties including its position dependence.As described in Ref. [28], the interaction cross section of optical photon in water as a function of the wavelength is characterized with the following three components: absorption (σ abs ), symmetric scattering (σ sym ) and asymmetric forward scattering (σ asym ).The absorption cross section, σ abs , is further modeled as a function of time (t) and vertical position in the detector (z) as, Here, the time-dependence of the position-averaged absorption cross section, σ 0 abs (λ, t) is inferred from the attenuation length measured with the decay-electron sample [25].The coefficient for the z dependence of the absorption, β(t), is computed using an empirical function of TBA in percent measured with the auto-Xe calibration system and water attenuation length in centimeters (L) as, where, The coefficients in the above formula were derived using simulated hits from a Ni-Cf source placed at the detector center.The survival probability of photons which originate from a vertex position (v x , v y , v z ) with a given wavelength λ is estimated by solving the following differential equation.N (r) is number of photons surviving at r, and C scat is an empirical correction factor to effectively characterize photon loss due to scattering.Based on a simulation study of mono-energetic electrons uniformly distributed in the detector, C scat = 0.6 is chosen as it minimizes the position dependence of the overall energy scale.
When a photon path is contained either within z ≥ −11 m or z ≤ −11 m, Eq. (A5) can be analytically solved.The survival probability after integrating over the path length of r i for the i-th PMT, p(r i , λ, t), is described as, (A6) where, Finally, the photon survival probability for the i-th PMT, P i (t), is evaluated by integrating p(r i , λ, t) over the wavelength, as, P i (t) = λmax λmin w 0 (λ)p(r i , λ, t)dλ, (A8) where w 0 (λ) is a weight function that represents the product of the Cherenkov light emission spectrum and the PMT's photon detection efficiency.In the previous study, this P i (t) was simply P prev i (t) = exp[−r i /L(t)] where L(t) is water attenuation length evaluated from the decay-electron sample.The new definition used for this analysis significantly reduces the position dependence of the energy scale.Combined with other improvements described in this section the variation of the energy scale across the detector was improved, characterized by a re-duction in the standard deviation from about 1.7% in the previous analyses to about 0.5% level in this analysis.

Appendix B: Signal extraction
The observed number of solar neutrino signal events are extracted from the final data sample using an extended maximum likelihood function fit.The distribution of θ Sun ( which is the angle between the reconstructed direction and the direction from the Sun at that time) is used from 3.49 MeV to 19.49 MeV in kinetic energy.The likelihood function is defined as Eq.(10), which is formed from the following quantities: b ij : The background probability density function for kth event in the i, j-th bin.This is obtained from the φ and θ distribution of data.For the day/night amplitude fit, this function is evaluated separately in five different solar zenith angle regions for the day and six different regions for the night.
s ij : The signal probability density function evaluated for the k-th event in the i, j-th bin.This function is extracted by the solar neutrino MC simulation.
: Y ij : The fraction of signal events in the j-th MSG bin in the i-th energy bin based on the MC simulation.This is obtained from the signal energy spectrum, which is the energy distribution of the solar neutrino MC final sample.The signal fraction is calculated from the number of signal events of each energy and MSG bin divided by the total number of signal events.
: B ij : Free parameter corresponding to the number of background events in the i-th energy bin and the j-th MSG bin.
: S: Free parameter corresponding to the total number of solar neutrino events in all energy and MSG bins.
The fitting parameters, S and B ij , represent the number of signal and background events, respectively.These parameters are obtained by maximizing the likelihood function.
The recoil electrons from solar neutrinos undergo multiple Coulomb scattering in water.This multiple scattering distorts the pattern of Cherenkov rings.For the higher energy region, the true energy range of the spallation background and solar neutrino events is similar, so multiple scattering has a similar impact on signal and background events.But when comparing events in a lower reconstructed energy region, the impact of multiple scattering is large for beta-decay backgrounds, such as bismuth-214 [33], because their true energy is often lower than the solar neutrino signal.Therefore we have introduced the MSG binning in Eq. (10).
Figure 58 and Figure 59 show the solar angle distributions of MSG sub-groups.As expected, the peak corresponding to the solar neutrino signal is more prominent in the high MSG-value bins but much smaller in the low MSG-value bins.This tendency is very clear in the lowenergy samples.
Figure 60 and Figure 61 also show the solar angle distribution of the observed events above 7.49 MeV, where the MSG sub-group is not used.
Appendix C: Supplementary explanation of the solar neutrino results

Spectrum results from SK
The expected rate of the solar neutrino interactions throughout the detector is estimated by considering the solar neutrino flux and the detector performance.The observed rate and the expected rate are summarized in Table IX.The ratios of the observed rate to the expected rate are summarized in Table X.

FIG. 1 .
FIG.1.The time variation of water transparency measured by decay electrons.This parameter is described as L in Appendix A.

FIG. 2 .
FIG. 2. (Color online) Top:The time variation of the topbottom asymmetry (TBA) throughout SK-IV full period.The red dots (blue upward-pointing triangle) show the TBA value measured by the auto-Xenon calibration system (Ni-Cf calibration source)[28].The black vertical line shows when the temperature of the supply water was changed, as described in the main text.Bottom: The difference between two TBA measurements, which is normally within ±0.5% during the SK-IV data set.

FIG. 3 .
FIG. 3. (Color online) The time variation of the relative PMT gain.Each panel shows PMTs from different production years.

FIG. 4 .FIG. 5 .
FIG. 4. (Color online) The time variation of the average dark rate of PMTs in each run.Each panel shows PMTs from different production years.

FIG. 6 .
FIG. 6. (Color online) The probability density function of the timing residual (ti − t tof − t0 = ∆t) for the single photoelectron signal.The second and the third peaks around 40 ns and 110 ns are caused by the PMT's after pulses.

FIG. 7 .
FIG.7.The vertex resolution as a function of the true electron kinetic energy.The events whose reconstructed vertex is in the standard fiducial volume (more than 200 cm from the ID wall) are sampled.

FIG. 8 .
FIG. 8. (Color online)The PMT acceptance as a function of the incident angle of the photon.The dashed-black line shows the acceptance in the case of the PMT itself while the solidred line shows that of the PMTs equipped with the acrylic cover and FRP case used after SK-II.

15 FIG. 11 .
FIG. 11. (Color online) Relative size of reconstructed N eff for simulated electrons with 10 MeV/c momentum generated uniformly in the entire ID volume.The top (bottom) panel shows results with the analysis method used for the previous (this) analysis.The vertical and horizontal axis represents reconstructed Z and R 2 = X 2 + Y 2 positions, respectively.

5 FIG. 12 .FIG. 13 .
FIG. 12. (Color online)The difference of the N eff between the calibration data and the MC simulation.The markers show the position of the calibration.The red (dashed) horizontal lines show ±0.5%.

5 FIG. 14 .FIG. 15 .FIG. 16
FIG. 14. (Color online) (Top) The energy distributions of the mono-energy electrons, reconstructed in the standard fiducial volume.(Bottom) The energy resolution [%] as a function of the electron's true total energy.

FIG. 17 .
FIG. 17.The typical opening angle distributions of the LINAC calibration data and the MC simulation in SK-IV.The black solid (gray dashed) line shows the angular resolution of data (MC), which contains 68.3% of events.The data samples are the same as Fig.15.

FIG. 18 .
FIG. 18. (Color online) The systematic uncertainties of the angular resolution estimated by the LINAC calibration.The definitions of the markers are the same as Fig. 16.

FIG. 19 .
FIG.19.The typical distribution of the DT calibration at position (X, Y, Z) = (+0.353,−0.707, 0.0) m taken on December 2016.The top panel shows the reconstructed kinetic energy of electrons, the middle panel shows the azimuthal angle dependence of the effective hits, and the bottom panel shows the zenith angle dependence of the effective hits.The effective hit is obtained by the peak position from the fitting with a Gaussian function.
FIG. 20. (Color online)The time variation of the effective hit (N eff ) of the DT calibration data and MC simulation (Top) together with their difference (Bottom).The vertical dashed line shows the period of the convection study.The two horizontal red-dotted lines show ±0.5% from the offset (greendot-dashed line).

FIG. 22 .
FIG. 22. (Color online) The directional dependence of the difference of the effective hit between the calibration data and the MC simulation.The value is obtained by considering the volume weight of the calibration position.The two horizontal red-dashed lines show ±0.1% from the offset (green-dotted line).
FIG. 23. (Color online)The time variation of the trigger efficiencies estimated by the MC simulation of 8 B solar neutrinos.The trigger efficiency is highly correlated with the gain of PMTs (Fig.3) and water transparency (Fig.1).After changing the threshold of data acquisition (vertical light-blue dashed line), the efficiency of 3.49-3.99MeV is improved significantly.

FIG. 24
FIG. 24.(Color online)Comparison of the log-likelihood for the signal and background of the non-neutron data period, without the WIT system.The spallation signal is shown by the red solid line, accidentals by the blue dotted line, and the tuned cut value is shown in the vertical green dashed line.The multiple spallation cut is already applied, so only 82% of remaining spallation needs to be removed to achieve 90% overall spallation removal effectiveness.

FIG. 25
FIG. 25. (Color online)Difference in solar neutrino signal events between 5.99-19.49MeV using the updated spallation cuts compared to the previous method.The red solid (blue dashed) line is the difference in total (background), respectively.The increase in background is from the change of cut tuning unrelated to spallation.

FIG. 26 .
FIG. 26.The time variation of event rates in the 2970-day final data sample.The vertical dashed line shows the threshold change from 34 hits to 31 hits.The events below 5.49 MeV during convection periods in 2018 are removed.

FIG. 28 .
FIG.28.Signal efficiency after each reduction step as a function of the reconstructed electron kinetic energy.The filledsquare corresponds to the number of events after the trigger cut in the 22.5 kt standard fiducial volume.Other markers are the same as Fig.27.

FIG. 29 .
FIG.29.The energy-correlated systematic uncertainties in SK-IV 2970-day data set.The gray, light-gray, and crosspattern show the systematic uncertainty of the 8 B spectrum shape, energy scale, and energy resolution uncertainties, respectively.

FIG. 30 .
FIG. 30.Solar angle distribution for the energy range of 3.49-19.49MeV using SK-IV 2970-day data set.The black points show the observed data, the dark and light gray shades show the background and the best-fit signal, respectively.

FIG. 33 .
FIG.33.Solar angle distribution for the energy range of 3.49-3.99MeV using SK-IV 2970 days data sample.The definitions are the same as in Fig.30.
FIG. 34.(Color online)Yearly solar neutrino flux measured by SK.The red-filled circle points show the SK data with statistical uncertainty and the gray striped area shows the systematic uncertainty for each phase.The horizontal black solid line (red shaded area) shows the combined value of measured flux (its combined uncertainty).The black-blank circle points show the sunspot numbers from 1996 to 2018.The sunspot numbers are taken from Ref.[67].

FIG. 35 .
FIG. 35.(Color online) Solar zenith angle dependence of the measured 8 B solar neutrino flux in SK-IV.The vertical bars show the statistical uncertainty.
FIG. 36.(Color online) Energy spectrum of the solar neutrino signal using the SK-IV 2970-day data set.The blue circle markers show the observed event rate as a function of the recoil electron energy.The red triangle markers show the expected rate from the 8 B + hep MC simulation without the oscillation effect.The signal efficiency shown in Fig. 28 is corrected.

FIG. 37 .FIG. 38 .
FIG. 37. (Color online) The measured energy spectrum in SK-IV.The blue points and bars show the observed rate divided by the expected event rate assuming no neutrino oscillation, with statistical and energy-uncorrelated uncertainties.The red bars and gray bands show the energy-uncorrelated and energy-correlated systematic uncertainties.ofFig.38shows the recoil electron spectrum taken in day or night time in SK-IV.Based on the spectrum results between day-time and night-time, the spectral straight day/night asymmetry is obtained as shown in the bottom

2 FIG. 41
FIG. 41. (Color online) SK constraints on solar oscillation parameters with (blue (medium gray)) and without (green (light gray)) using the absolute elastic scattering rate.The shaded regions are for 2σ confidence level.Also shown are 1σ and 3σ allowed areas.The dotted (dashed) contours are for 4σ (5σ) C.L.

FIG. 43 .
FIG. 43.(Color online) Comparison of the Recoil Electron Spectra of SK-I (top left), SK-II (top right), SK-III (bottom left), and SK-IV (bottom right) with MSW predictions.The green(light gray) (blue(medium gray)) histograms are for the solar (solar+KamLAND) best-fit ∆m 2 21 value.The thick (thin) lines are distorted (not distorted) by adjusting the energy scale, resolution, and neutrino spectrum within uncertainties.

1 FIG. 49
FIG. 49. (Color online) SK-IV day/night amplitude fit dependence on ∆m 2 21 .The black line (grey band) shows the best-fit value of day/night asymmetry (its uncertainty).The red solid curve shows the expected day/night asymmetry.The green solid (blue dashed) box shows the 1σ range allowed by solar experiments (solar experiments and KamLAND).

2 FIG. 50 .
FIG. 50.SK-IV day/night amplitude fit dependence on energy.The grey band shows the combined value.

Figure 51 compares
Figure 51 compares the results of the day/night flux asymmetry measured at each SK phase.Combining all SK data of day/night flux amplitude fit results in A SK, fitD/N = −0.0286± 0.0085 (stat.)± 0.0032 (syst.).(31) Here, the asymmetry parameter is expressed based on the SK-I energy range, so the expected asymmetry is a bit stronger: −0.0242.Since this result differs from zero by 3.2σ, we find evidence for the existence of earth matter effects on solar neutrino oscillation.The asymmetry parameter depends on the expected zenith angle

FIG. 51
FIG. 51. (Color online) Day/night asymmetries from amplitude fits by the four different SK phases.The black-thick (redthin) bar shows the statistical (total) uncertainty and the grey band shows the combined value.

FIG. 53
FIG. 53.(Color online)The electron neutrino survival probability as a function of neutrino energy obtained from all SK data.The blue(medium gray) (red(dark gray)) region obtained form the Pee,quad(Eν) (Pee,exp(Eν)) function.The green(light gray) region is from the Pee,cubic(Eν) function.The thick blue(medium gray) line is the Pee(Eν) distribution expected from neutrino oscillation at all solar + KamLAND best-fit parameter set, and the thick green(light gray) line is the Pee(Eν) distribution expected from neutrino oscillation at all solar best-fit parameter set.

FIG. 54
FIG. 54. (Color online)The electron neutrino survival probability as a function of neutrino energy.The green(light gray) (blue(medium gray)) region is obtained from the Pee,quad(Eν) function with SK (SNO) spectrum data.The red(dark gray) region is obtained from the same function, but with SK+SNO data.The thick blue(medium gray) and green(light gray) lines are the same as Fig.53.

)
When the path crosses the z = −11 m boundary, this survival probability is calculated separately for z ≥ −11 m and z ≤ −11 m regions.Then they are multiplied together to evaluate the attenuation over the entire path.

N
energy : Total number of energy bins The kinetic energy from 3.49 MeV to 19.49 MeV is divided into twenty-three bins (N energy = 23).The width of the energy bin is defined as follows: / bin ( 3.49 ≤ E < 13.49MeV) 1.0 MeV/ bin (13.49≤ E < 15.49MeV) 4.0 MeV/ bin (15.49≤ E < 19.49MeV) N MSG i : Total number of MSG bins Below an electron kinetic energy of 7.49 MeV, the MSG parameter (g MS ), which runs from 0.0 to 1.0, is divided into three bins in each energy bin (N MSG i = 3).In higher energy bins, no separation by the MSG parameter is applied (N MSG i = 1).The width of the MSG bins are defined as follows: < g MS < 0.35 (E < 7.49 MeV) 0.35 < g MS < 0.45 (E < 7.49 MeV) 0.45 < g MS < 1.00 (E < 7.49 MeV) 0.00 < g MS < 1.00 (E ≥ 7.49 MeV)n ij : Total number of events in the i-th energy bin and in the j-th MSG bin.
Currently at Department of Physics and Wisconsin IceCube Particle Astrophysics Center, University of Wisconsin-Madison, Madison, WI 53706, USA.i Currently at Division of Science, National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka-shi, Tokyo, Japan.j Currently at Department of Physics, Boston University, Boston, MA 02215, USA.k Currently at Department of Physics, Illinois Institute of Technology, Chicago, 60616, IL, USA.l Currently at Department of Physics and Astronomy, York University, Toronto, Ontario, Canada.m Currently at EP Department, CERN, 1211 Geneva 24, Switzerland.
a Currently at Research Center for Neutrino Science, Tohoku University, Sendai 980-8578, Japan.b Currently at High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan.c Currently at School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai, China.d Currently at Fermi National Accelerator Laboratory, Batavia, IL 60510, USA.e deceased f Currently at School of Physics and Astronomy, University of Glasgow, Glasgow, G12 8QQ, United Kingdom.g Currently at Department de Fisica Teorica, Universitat de Valencia, and Instituto de Fisica Corpuscular, CSIC, Universitat de Valencia, 46980 Paterna, Spain.h n Currently at SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025-7090, USA o Also at BMCC/CUNY, Science Department, New York New York, 1007, USA.p Also at University of Victoria, Department of Physics and Astronomy, PO Box 1700 STN CSC, Victoria, BC V8W 2Y2, Canada.

TABLE I .
Experimental phases of SK.The live times are the total duration of the good observation periods for solar neutrino analysis.The energy threshold is based on the electron kinetic energy.

TABLE II .
Summary of the systematic uncertainties of the energy scale in SK-IV

TABLE III .
Summary of the data set in SK-IV.During the two convection study periods, a higher reconstructed energy threshold of 5.49 MeV is used.

TABLE VI .
Summary of systematic uncertainties for the day/night amplitude fit.
× 10 −5 eV 2 .Color online) The electron neutrino survival probability as a function of neutrino energy and solar neutrino measurements.The red region is obtained from the Pee,quad(Eν) function with SK+SNO spectrum data.The thick blue and green lines are the same as Fig. 53.The light blue (filled triangle) and gold (open triangle) data points are the average pp and CNO neutrino survival probabilities inferred from the radio-chemical solar neutrino data as well as SK+SNO and Borexino 7 Be measurements, respectively.The dark blue (filled square) and dark red (open circle) data points represent the average pp and 8 B neutrino survival probabilities, respectively.The green (open square) and turquoise (filled circle) data points represent the 7 Be and pep neutrino survival probabilities from Borexino data, respectively.