Direct observation of coherent elastic antineutrino–nucleus scattering

Data taking

Data from the CONUS+ run 1 analysis presented in this study correspond to the period from November 2023 to July 2024. Three of the four detectors, named C2, C3 and C5, are considered for analysis with a total active mass of 2.83 ± 0.02 kg (ref. ²). After selection cuts to remove time periods with unstable noise conditions, deficient radon flushing and a few days with a strong contribution of microphonic events, the exposure considering the active mass is 327 kg days with reactor on and 60 kg days with the reactor off.

The evolution of the main environmental parameters during reactor on and off is shown in Extended Data Fig. 1 for the three detectors. The shape of the noise peak reconstructed in the lowest channels of the data acquisition system (DAQ) was found to follow a Gaussian distribution¹ and its full-width at half maximum (FWHM) was monitored over time with variations below 2 eV_ee (electron equivalent energy). The dataset was selected to ensure that the noise rate variations are below 20%.

In a dedicated study, we confirmed the correlation between noise rate and cryocooler power. Stability was improved compared with the CONUS setup in Brokdorf (KBR) by replacing the two-fan ventilation system with a water-cooled chiller system for the pulse tube cryocoolers¹. Even at cryocooler power variations of up to 30 W, no correlation was observed with the count rate in the region of interest above 160 eV_ee.

Detector response and energy scale

The trigger efficiency was determined by injecting artificial signals produced by a pulse generator with the same rise time as the physical signals. The pulses are injected through a specific circuit implemented in the HPGe preamplifier. A detailed scan allowed us to measure the detector response as a function of the energy¹. The trigger efficiency remains more than 90% down to 140 eV_ee for all detectors. The evolution of the trigger efficiency curve parameters during run 1 was studied with different measurements, remaining stable with differences of less than 2% throughout the run.

The energy was calibrated using the binding energies of the K-shell (10.37 keV) and L-shell (1.30 keV) from the decays of ⁶⁸Ge/⁷¹Ge inside the HPGe diodes, considering a linear behaviour in this energy range. At shallow depth, the ⁶⁸Ge/⁷¹Ge radioisotopes are continuously produced by cosmic and muon-induced neutrons¹. Thus, it is possible to monitor the stability of the energy scale during the whole measurement of this in situ activation, observing variations below 2%. A specific ²⁵²Cf irradiation was performed at the end of run 1, collecting in 45 days of the measurement more than 5,000 events from the K-shell and 700 events from the L-shell in each detector. In this way, an energy calibration uncertainty below 5 eV_ee is achieved. The HPGe diode also produces 158 eV X-rays from the binding energy of the M-shell of the ⁶⁸Ge/⁷¹Ge decays. Using the ratio of the K- and L-shells, a total of 100 events per detector was expected after irradiation with ²⁵²Cf. Although an indication of such events is seen, no conclusive signal was observed because of the proximity to the noise edge and the lack of statistics.

The linearity of the DAQ chain in the sub-keV_ee region was investigated using the pulse generator signals. The results are shown in Extended Data Fig. 2. Deviations from a pure linear behaviour are observed below 250 eV_ee, at which a few eV_ee variations have a strong impact on the CEνNS signal. They were attributed to two DAQ-related effects³⁴. The energy non-linearity was corrected in the CEνNS analysis and was taken into account during the event energy reconstruction.

The energy resolution of the detectors at low energies was also evaluated with the pulse generator signals. For the C2 detector, a resolution of (48 ± 1) eV_ee (FWHM) was found, whereas the C5 and C3 detectors have a resolution of (47 ± 1) eV_ee (FWHM). Further contribution to the total energy resolution comes from the statistical fluctuation of the number of electron–hole pairs produced in the Ge crystals in the case of an event, which is

with ϵ = 2.96 eV, the average energy needed to create a single electron–hole pair in Ge, and the Fano factor F = 0.11.

Selection cuts

Selection cuts are applied to reduce background events while keeping the CEνNS signal. Four different selection cuts are applied to the data. First, the muon-veto system allows for efficient suppression of the impact of cosmic radiation, using a 450-μs anti-coincidence window between the veto and HPGe signals. The average rate detected in the muon veto during reactor on is (274 ± 1) Hz and decreases to (214 ± 1) Hz in reactor off periods. The corresponding average dead times in both periods are 12.3% (reactor on) and 9.6% (reactor off). Second, inhibit signals are generated when the increasing baseline has reached saturation of the dynamic range of the transistor reset preamplifier for each HPGe detector¹. An anti-coincidence window of 1–2.5 ms (depending on the detector) is applied to veto unwanted spurious HPGe detector signals generated shortly after the resets. This cut suppressed 30% events after the muon veto anti-coincidence at low energy, becoming negligible above 5 keV_ee. The dead time induced by this cut is calculated, combined with the previously mentioned muon-veto dead time, to avoid the overlapping of both veto windows. An additional dead time of 0.5–2.1% (depending on the detector) is estimated. Third, the time difference distribution of events is studied in each channel as proposed in ref. ². Finally, an anti-coincidence cut is applied between different HPGe detectors with a 5-ms time window. The probability of a neutrino interacting with different detectors is negligible, whereas for other backgrounds, such as muon-induced neutrons created in the shield, simultaneous hits in several detectors at once can be expected.

The rejection efficiencies of these selection cuts are summarized in Extended Data Table 1 for the three detectors in different energy regions. The total dead time induced by the selection cuts is between 12.8% and 14.4% in the reactor on periods.

Background model

The background model used in the analysis of the CONUS+ data is based on Monte Carlo simulations using the Geant4-based framework MaGe³⁹, following the approach in ref. ²⁹. A complete decomposition of the background in both reactor on and off data was done for all detectors used in the analysis. In the following, the main sources of background are described.

Cosmogenic neutrons

As described in ref. ³⁰, the impact of cosmogenic neutrons with energies up to 100 MeV was investigated by first propagating the expected neutron spectrum at the KKL location (based on refs. ^40,41) with a flux of (1.4 ± 0.2) × 10⁻² neutrons s⁻¹ cm⁻² through a reactor building model. Then the resulting flux is tracked inside the CONUS+ room. This flux was used as the basis for the next simulation, in which neutrons were started isotropically from a half-sphere around the CONUS+ shield, and their contribution to the CONUS+ background was measured. The simulations show a contribution of (21.6 ± 3.1) counts day⁻¹ kg⁻¹ in each detector in the energy range between 0.4 keV_ee and 1 keV_ee, which corresponds to approximately half of the background counts in this region.

Cosmogenic muons

For the muon simulations²⁹, the expected flux and muon spectrum at an overburden of 7.4 m water equivalent (w.e.) were calculated from those at the surface of Earth^42,43 and propagated through the shield. The resulting spectrum was validated by comparing it with the CONUS+ data without the applied muon-veto cut, which showed good agreement. The muon-veto cut was then applied by multiplying the simulation output by a factor of 0.01 for all energies greater than 2 keV_ee, corresponding to a muon-veto efficiency of 99%.

For energies below 2 keV_ee, a different approach was taken. Here, simulations show an inefficiency in the tagging ability of the muon-veto system because of the setup of the shield. The outer muon veto is located under a layer of lead in the CONUS+ shield. As such, it is possible for muons to pass through this outermost lead layer without hitting one of the muon-veto layers. These muons can induce electromagnetic showers in the outer lead layer, which propagate through the shield and are registered in the Ge detectors. However, because no muon passes through any muon-veto layer in such an event, the energy deposition in the plastic scintillator plates will be much lower, resulting in a greatly reduced tagging efficiency of these events. Simulations show that at energies below 0.4 keV_ee, up to 80% of all muon-induced background comes from these events. Based on this simulation output, this inefficiency was modelled with a polynomial and accounted for in the final muon-veto efficiency. The resulting efficiency drops towards lower energies, with its minimum being 97% below 0.4 keV_ee. Using this approach, the overall background contribution of the cosmic ray muons was found to be (17.4 ± 0.3) counts day⁻¹ kg⁻¹ in each detector in the energy range between 0.4 keV_ee and 1 keV_ee, which corresponds to approximately a third of the background counts in this region.

Leakage test component

During the final run of the CONUS experiment at KBR, an additional background component had to be introduced in the background model²⁸. This component was present after ventilation of the cryostats with argon gas to avoid mechanical deformation during a regular leakage test at KBR in July 2019 (ref. ²⁹). Of the four detectors used in the CONUS experiment, two (C2 and C3) are used for the analysis presented in this work and are affected by this additional background. The simulations and the background model show that an additional component with the same shape is still present in the background of these two detectors, but is absent in C5, which was not at KBR at the time. This additional background is constant during reactor on and off periods. Therefore, the leakage test component was again included in the background model of C2 and C3 by modelling it using a function with two parameters, as in ref. ²⁸. The resulting impact is below 10% in the energy region between 0.4 keV_ee and 1 keV_ee.

Other background components

The remaining background in each detector is made up of many different components, similar to the situation in KBR²⁹. There are no hints that the detector upgrades or the movement of the setup were introducing any contamination. The background model includes cosmogenically induced isotopes in the copper parts (⁵⁷Co, ⁶⁰Co and ⁵⁴Mn) of the cryostat and the Ge crystals (⁵⁷Co, ⁶⁸Ge, ⁶⁸Ga, ⁶⁵Zn and ³H), Radon inside the detector chamber, ²¹⁰Pb inside the cryostat and shield, metastable Ge states (^71mGe, ^73mGe and ^75mGe) and inert gases coming from the reactor (⁸⁵Kr, ¹³⁵Xe and ³H). The results of the simulation of these components were scaled to be in accordance either with the rates of gamma lines produced by the decay of these isotopes in the spectrum (for example, for radon) or with screening measurements performed before the installation. The listed contributions are typically very subdominant in the region of interest, with the decay of radon inside the detector chamber being the only exception. This contribution results in a background rate of (1.9 ± 0.1) counts day⁻¹ kg⁻¹ between 0.4 keV_ee and 1 keV_ee for the C5 detector (C2: (2.8 ± 0.1) counts day⁻¹ kg⁻¹; C3: (2.6 ± 0.1) counts day⁻¹ kg⁻¹). These values correspond to approximately 5% of the background in this energy region. Radon decays have a high impact on energies above 100 keV_ee, at which they can contribute up to 60% of the measured background. Slow pulses arising from decays on the surface of the diode and in the transition layer are included in the model²⁹. Energy depositions from these events can be stopped within the transition layer, and the released charge diffuses slowly into the active volume, resulting in long rise times and incomplete charge collection. Their impact is accounted for by registering the exact coordinates of an interaction in the Ge crystal. If the coordinates place it within the transition layer of the crystal, the energy of the event is shifted towards lower energies using a sigmoid-like function. Details on this procedure can be found in ref. ².

Model differences in reactor off data

The background model accounts for the experimental differences during data collection with the reactor off. The first of these differences is induced by the fact that during a reactor outage, the drywell head from the containment structure surrounding the reactor core is placed directly above the CONUS+ room. This drywell head is made of 3.8 cm steel and, therefore, increases the overburden of the experiment by approximately 0.3 m w.e., which results in a reduction of 19% in the flux of cosmogenic neutrons and a reduction of 3% in the flux of cosmic ray muons. As a result, these two background contributions are reduced accordingly. The second difference in the background model of the reactor off period comes from a more effective removal of radon in the detector chamber. During the course of run 1, shield tightness and radon-free air flushing were improved. As a result, the radon contribution in reactor off time is reduced by approximately a factor of 4–6 compared with that for reactor on. The radon contributions are scaled to match the count rates in the gamma lines induced by the radon decay. Moreover, reactor-correlated background components, such as reactor neutrons and high-energy gammas, for example, from ¹⁶N, were investigated. Their impact was found to be negligible in all energy regions, including the region of interest. The resulting background model for all three detectors can be seen in Fig. 2 and Extended Data Fig. 3. The time-dependent contributions of the M-shell line of the ⁶⁸Ge/⁷¹Ge decays are included in the background model. The difference between the M-shell contribution in reactor on and off periods was estimated based on the measured K-shell count rates. This difference was found to be below 1% of the observed neutrino signal.

In Extended Data Fig. 4, the reactor on and off datasets are compared with the background models. If the calculated difference between on and off phases is added to the reactor off data, good agreement with the measured on data is found above the signal region. For all detectors, background models and data are fully consistent above 0.4 keV_ee. As the data in the reactor off phase are statistically limited, the background model still plays an important part in the analysis. Extended Data Fig. 5 shows the fractional contribution of the main components to the total background rate as predicted in the model.

Quenching

The ratio of the ionization energy released by nuclear recoil in a CEνNS event and the ionization energy of electrons of the same energy is given by the quenching factor. In the CONUS+ analysis, the energy-dependent signal quenching is described by the Lindhard model³² with a quenching parameter k = (0.162 ± 0.004) as determined in ref. ³⁴. Alternative quenching descriptions are also tested as described in the supplementary material of ref. ²⁸, including a linear and cubic functional form to describe the increased quenching factor compared with the Lindhard theory found in ref. ³³. It was shown that the Migdal effect is subdominant with respect to the CEνNS signal in our region of interest³⁶.

The signal predictions for the different quenching descriptions are shown in Extended Data Fig. 6, together with the difference between the data in the reactor on phase and the background model. A signal rate of (2,600 ± 300) events is expected for the linear function, whereas for the cubic function, (550 ± 50) events are predicted. Both numbers are significantly higher than the neutrino rate extracted from the CONUS+ data. The standard Lindhard model provides the best description of the reactor-correlated excess at low energy. In Extended Data Fig. 6, we have included an additional data point at 200 eV_ee, which is obtained using the information shown in fig. 4 of ref. ¹⁹. There, a CEνNS signal was approximated with an exponential using two parameters, the amplitude at 200 eV_ee (A_0.2) and a decay constant. The favoured value for A_0.2 is shown within a 2σ contour in a two-dimensional plot of the two parameters. We scaled this A_0.2 value to our exposure and corrected for the difference in neutrino flux. As expected, it matches the description of Lindhard with a linear increase added at low energies, but it is in clear conflict with the CONUS+ data. The lower error bar corresponds to the smallest value of A_0.2 in the 2σ contour. This value is also ruled out by the CONUS+ data points.

Likelihood fit and systematic treatment

A likelihood function is used to determine the CEνNS signal in the CONUS+ run 1 data.

where \({{\mathcal{L}}}_{{\rm{O}}{\rm{N}}}\) and \({{\mathcal{L}}}_{{\rm{O}}{\rm{F}}{\rm{F}}}\) are the binned likelihood functions for reactor on and off periods. Gaussian pull terms for the systematic uncertainties are also included. Here, the first term represents the pull terms for correlated parameters, namely, the trigger efficiency parameters, and the second term represents pull terms for uncorrelated parameters, such as the active mass of the detectors, the reactor neutrino flux and the uncertainty on the energy scale calibration of the spectra. Parameters that were experimentally determined are pulled to their measured values. The background scaling factor b, an additional fit parameter for the overall normalization of the background model, is also included and pulled to 1. The binned likelihood functions assume a Poisson distribution and have the form

where N is the number of bins in the region of interest, n_i is the bin content in the histograms of the measured data and μ_i is the bin content of the model. For a single detector, the model μ for reactor on data is calculated from the sum of the background model n^b (scaled with b) and the predicted CEνNS spectrum n^s (scaled with the signal parameter s) by taking into account the live time of the experiment (t_ON and t_OFF), the active volume of the detector (m_act), the dead-time correction (c_dt), the fission flux (θ₂) and a multiplication factor (θ₃), summarizing uncertainties of the detector response. In total,

The combined fit minimizes \(-\log {{\mathcal{L}}}_{{\rm{ON}}}\) and \(-\log {{\mathcal{L}}}_{{\rm{OFF}}}\) for all three detectors simultaneously, whereas the signal parameter s, which indicates the number of CEνNS counts, is shared among the detectors.

Fits using just single detectors independently give consistent results as shown in Extended Data Table 2. To cross-check the result, fits were performed by two independent likelihoods with different approaches concerning the uncertainty on the quenching model. Likelihood A applied a predicted CEνNS spectrum with a fixed k parameter (k = 0.162), while introducing a fourth-order polynomial to vary the shape of the signal spectrum, applying Gaussian pull terms on each parameter. In likelihood B, the k value of the Lindhard model is a fit parameter with a pull term. Moreover, there are some differences in the treatment of the non-linearity corrections and the minimization algorithms between the two likelihood fits. The results of both likelihood implementations agree within 2%.

For the signal prediction, we choose the Helm parameterization of the nuclear form factor^44,45 and a reactor antineutrino spectrum based on the data-driven approach in ref. ³¹. The given antineutrino spectra are adjusted to the fission fraction of the Leibstadt reactor and augmented with measured antineutrino spectra at energies above 8 MeV (ref. ⁴⁶) and simulation data below 1.8 MeV (the threshold of inverse beta decay)⁴⁷. We assume no uncertainty on the nuclear form factor, but account for a 2% uncertainty on the weak mixing angle⁴⁸, leading to an overall uncertainty of 3.2% on the CEνNS cross-section. The shape of the applied antineutrino spectra, the thermal power of the reactor, fission fractions and energy releases per fission of the relevant fission isotopes contribute a combined 4.6% error to the expected event rate. Further uncertainties on quenching, detector active mass, trigger efficiency and the energy threshold are taken into account and listed in Extended Data Table 3. All contributions lead to an overall uncertainty of 17% on the signal prediction. To be conservative, we assume prediction uncertainties are fully correlated between single detectors. Currently, the dominant uncertainty is from the energy calibration. Future extended calibrations will allow us to reduce this uncertainty from 5 eV_ee to 3 eV_ee.

The likelihood fit itself gives a result of (395 ± 86) CEνNS counts in the combined fit of the three detectors. Additional systematic uncertainties not implemented as pull terms are evaluated in a second step and added in quadrature, giving the final uncertainty of  ±106 signal counts, as shown in the Extended Data Table 3. The non-linearity term was obtained by varying the correction parameters and checking the impact on the likelihood result. The uncertainties of the calibration points in Extended Data Fig. 2 were used to generate the parameter variations. A Gaussian fit was performed on the distribution of the central values of the new likelihood results, where its 1σ value was taken as additional systematic uncertainty. The term related to the background model was studied in the same way by varying muon flux, neutron flux and the leakage component. The variations were 14% (taken from ref. ⁴⁰), 6% (based on ref. ⁴³) and 10%, respectively. Moreover, the systematic uncertainty of the fit method was estimated by the difference between likelihood A and B.

The fractional contribution from antineutrinos of different energies to the signal expectation is quantified in Extended Data Fig. 7. At the current detector threshold, we are sensitive to antineutrino energies above 5 MeV. With lower detection thresholds, the steep rise in antineutrino flux towards low energies will result in markedly higher signal expectations¹.