Direct measurement of the ionization quenching factor of nuclear recoils in germanium in the keV energy range

This article reports the measurement of the ionization quenching factor in germanium for nuclear recoil energies between 0.4 and 6.3 keV$_{nr}$. Precise knowledge of this factor in this energy range is relevant for coherent elastic neutrino-nucleus scattering and low mass dark matter searches with germanium-based detectors. Nuclear recoils were produced in a thin high-purity germanium target with a very low energy threshold via irradiation with monoenergetic neutron beams. The energy dependence of the ionization quenching factor was directly measured via kinematically constrained coincidences with surrounding liquid scintillator based neutron detectors. The systematic uncertainties of the measurements are discussed in detail. With measured quenching factors between 0.16 and 0.23 in the [0.4, 6.3] keV$_{nr}$ energy range, the data are compatible with the Lindhard theory with a parameter $k$ of 0.162 $\pm$ 0.004 (stat+sys).


Introduction
Recent advances in the detection of very low energy signals opened up new possibilities to study coherent elastic neutrino-nucleus scattering (CEνNS) or to detect low mass dark matter candidates. Their observation relies on the detection of nuclear recoils in the keV energy range, where energy losses result from the combination of ionization and atomic collisions. Depending on the chosen detector technology, the associated signal -collected in the form of ionization energy, scintillation light or heat -may differ from the one obtained from gamma calibration sources of the same energy. It is therefore of primary importance to precisely know the detector response to nuclear recoils in order to accurately reconstruct their energy. For detectors relying only on the detection of the ionization signal, such as High-Purity Germanium Detectors (HPGe) operated at liquid nitrogen temperature, one quantity of specific interest is the dimensionless ionization quenching factor defined as the ratio of the ionization energy generated by nuclear recoils over the one generated by electron recoils of the same energy. This quantity has been extensively measured for nuclear recoils in the 10-100 keV range [1][2][3][4][5][6] and it follows the energy dependence predicted by Lindhard et al. [7]. For experiments aiming at detecting reactor antineutrinos via CEνNS in HPGe detectors [8][9][10][11], recoils of a Corresponding author: aurelie.bonhomme@mpihd.mpg.de, conus.eb@mpi-hd.mpg.de a few keV nr (nuclear recoil energy) are expected, producing ionization signals in the sub-keV ee (electron equivalent energy) range. At these energies, precise experimental measurements of the ionization quenching factor are still lacking and the overall validity of the Lindhard theory has been questioned [12][13][14][15]. Moreover, recent measurements have shown a deviation with respect to the prediction in silicon [16,17] and germanium [18,19], although previous measurements [20][21][22][23][24] below 10 keV nr were in agreement within uncertainties.
One of the experiments looking for the yet not observed CEνNS signal from reactor antineutrinos is CONUS. It consists of four 1kg-sized HPGe detectors [8] deployed at the nuclear power plant of Brokdorf (Germany). For CONUS, the quenching factor is not only crucial for the interpretation of the Standard Model neutrino data [25], but it is also the current main systematic uncertainty in the search for new physics [26]. Because of their threshold of ∼ 250 eV ee , experiments like CONUS will be able to detect the ∼ 1.5 keV nr recoils induced via CEνNS by the highest energy reactor antineutrinos if the quenching factor is larger than 0.15. The experimental determination of this quenching factor in this range is therefore tremendously important to support the on-going and future experimental program since it determines directly the sensitivity of these experiments and their ability to perform precise measurements of the CEνNS signal. This article reports a direct measurement of the ionization quenching factor of nuclear recoils in germanium for keV recoil energies. A low mass HPGe detector was positioned as an active target in monoenergetic neutron beams at the PTB Ion Accelerator Facility (PIAF) [27] of the Physikalisch-Technische Bundesanstalt (PTB) in Braunschweig (Germany). Liquid scintillator (LS) based neutron detectors allowed to select nuclear recoil energies between 0.4 keV nr and 6.3 keV nr in the germanium target via the coincidence detection of the scattered neutrons. The article is structured as follows: after the experimental setup description in Sect. 2, the analysis is detailed in Sect. 3 and the quenching results are presented in Sect. 4. The systematic uncertainties are thoroughly discussed and as a further cross-check a comparison of the integrated recoil spectrum with a Monte-Carlo (MC) simulation is included.
All uncertainties reported in this article correspond to one standard deviation, i.e. with a coverage factor of 1.

Experimental setup 2.1 Layout of the experiment
The experiment was set up at beam line 2 in the experimental hall of the PTB ion accelerator facility PIAF. Fig.  1 shows a sketch of this setup. The HPGe detector acts as an active target deployed in the neutron beam, produced by proton-induced reactions in a lithium target. Neutrons scattered from the Ge nuclei are detected in LS detectors in coincidence with the ionization signal produced in the HPGe detector. Thus, the energy deposited in the HPGe target can be traced back using the kinematical relations for elastic neutron scattering. The LS detectors had to be shielded against the neutron source. Therefore, the lithium target was inserted into a massive shield made of borated polyethylene with a conical opening aligned with the nominal axis of the proton beam. A lithium-6 glass scintillation detector ( 6 LiGl) and a barium fluoride scintillation detector (BaF 2 ) were employed to monitor the neutron beam by time-of-flight (TOF) measurements. The different components of this setup, their characterization and calibration are detailed in the following subsections.
The analysis presented in this article is based on data collected in autumn 2020. Details of the runs are given in Tab. 1. Prior to the actual measurement of the quenching factor (runs 1 to 12 corresponding to different neutron beam energies), the LS detectors were characterized in a separate experimental campaign at another PIAF beam line (runs a, b, c).

Neutron beams
The neutron beams used for this experiment were produced by bombarding metallic lithium targets with proton beams from the 2 MV Tandetron accelerator of PIAF. The advantage of metallic lithium target compared to the commonly used lithium fluoride targets is the increased The proton beam is coming from the right and hits the Li target, producing a neutron beam in the forward direction. Neutrons scatter off germanium nuclei in the HPGe target detector and are detected in coincidence with the LS detectors placed behind in a half-circle. Monitoring detectors (Mon.) placed at 0 degrees allow to measure the neutron energy distribution of the beam. neutron yield per unit target thickness and the reduced yield of parasitic high-energy photons from 19 F(p,αγ) 16 O reactions which could cause an unwanted background increase. The proton beams were produced in pulsed mode with a repetition frequency of 1.25 MHz using the chopper/buncher system in the injection beam line. The resulting proton beam had pulse durations of about 2.5 ns and beam currents around 1 µA. A reference signal for TOF measurements was derived from an inductive beam pick-up located close to the neutron production target. For proton energies between the threshold of 1.88 MeV and 2.37 MeV, the 7 Li(p,n) 7 Be reaction produces only monoenergetic neutrons. For higher proton energies, transitions to the first excited state in the residual 7 Be nucleus result in a second neutron contribution of lower energy. The energy of the proton beam was measured using a 90 o analyzing magnet, which is regularly calibrated using a set of resonances and reaction thresholds.
A direct measurement of the neutron energy distribution using the TOF method was chosen as reference for the present experiment. Two monitoring detectors were employed for this purpose: a 6 LiGl detector and a fast BaF 2 detector. The 6 LiGl detector has a diameter of 38 mm and a thickness of only 3 mm. The small thickness and diameter make this detector suitable for TOF measurements in the 100 keV to 1 MeV energy range because the crossing  Table 1: Parameters of the neutron beams used for the characterisation of the LS detectors and for the measurement of the ionization quenching factors. The nominal proton and mean neutron energies are denoted by E p,0 and E n,0 , respectively. The experimentally measured neutron energy E n,exp from the two employed calibration detectors is reported and the adopted value for the analysis in this article is named E n . The difference between the proton energy E p calculated from E n,exp and the nominal proton energy E p,0 is denoted by δE p . The nominal lithium mass per unit area used for this calculation is m Li .
time is less than 0.7 ns. In this detector, neutrons are detected via the 6 Li(n, t) 4 He reaction. The large Q-value of 2468 keV facilitates a separation of neutron-and photoninduced events by a simple pulse-height threshold. Fig.  2 shows a two-dimensional pulse-charge versus TOF distribution for a 2240 keV proton beam corresponding to a nominal neutron energy of 500 keV. The random background of neutron-induced events above the pulse-height threshold is mainly due to epithermal room-return neutrons. Due to the 1/v-dependence of the cross section, the detection efficiency for these neutrons is much larger than for the fast neutrons. Without the large experimental hall at PIAF and the low-mass grid floor this background would have been much more intense. The time resolution of the 6 LiGl detector is about 4 ns, i.e. larger than the duration of the proton beam pulses. Therefore, the fast BaF 2 detector (51 mm in diameter and length, time resolution of about 0.5 ns) was used to optimize the time structure of the beams. This detector can also detect neutrons via inelastic scattering and the detection of the deexcitation photons. Thus, it can be used to cross-check the 6 LiGl measurements, but the length of the detector made it necessary to simulate the time distribution of inelastic scattering events in the detector using the MCNPX code [28] and include it in the analysis. The distance between the center of the monitoring detectors and the neutron production target was measured with a laser distance meter and varied between 650 mm and 1540 mm with an uncertainty of 2 mm.
The characteristics of all neutron beams used for the present experiment are summarized in Tab. 1. The expected mean neutron energy E n,0 calculated from the nominal proton energy E p is significantly higher than the experimental mean neutron energy E n,exp . E n,0 was calculated from the proton energy by taking into account the energy loss in the target, ranging between 9 keV and 16 keV for the proton energies used for the present experiment. The uncertainty of E n,0 mainly results from the uncertainty of the mass per unit area m Li , determined via an In this way, they are discriminated from the ambient γ-ray background associated to a low charge.
indirect measurement using co-evaporation of lithium on a quartz balance during the preparation of the targets by physical vapour deposition (PVD) with typical uncertainties of about 10 %. The uncertainties of the experimental neutron energy E n,exp reflect the uncertainties of the peak positions of the gamma and neutron peaks in the TOF distributions, the uncertainty of the flight distance and the uncertainty resulting from the differential non-linearity of the TOF measurement. For all proton beams for which measurements with the two monitor detectors were available, the measured mean neutron energies agreed within their uncertainties. The difference δE p between the proton energy calculated from the experimental mean neutron energy and the nominal proton energy is about -10 keV, except for the runs 3 to 6, where the nominal and the experimental neutron energy agreed within their uncer-  tainties. Similar discrepancies between the nominal and the experimental neutron energy were already observed in earlier experiments. Most likely, they are due to a not completely understood technical issue in the beam transport system at PIAF. For the analysis of the present experiment, the experimental mean neutron energies measured with the 6 LiGl detector were adopted whenever available.
For the other cases, the experimental energy measured with the BaF 2 detector was used. The only exception is the 800 keV neutron beam for which no experimental measurement of the mean neutron energy was available. Here, it was assumed that the effective proton energy was 10 keV less than the nominal proton energy of 2519 keV, resulting in a mean neutron energy of 790 keV for a lithium width of 100 µg cm -2 . A conservative estimate of (790 ± 11) keV was adopted for this dataset. At the position of the HPGe detector 20 cm behind the exit of the collimator, the neutron beam had a diameter of about 4 cm. The latter was estimated via a scan of the neutron counting rate along the beam axis performed with the 6 LiGl detector.

HPGe target detector
A dedicated n-type HPGe detector was designed in cooperation with the company Mirion Technologies (Canberra) SAS in Lingolsheim (France) to fulfill three main requirements. A small detector mass (10 g) was chosen in order to minimize multiple neutron scattering inside the crystal which would smear out the fundamental reconstruction of the neutron scattering angle. Moreover, such a small capacitance detector garantees low noise and thus a low energy threshold. The crystal is 6 mm thick and segmented into four quadratic pixels, each with dimensions 9×9 mm 2 . For this geometry and neutron energies between 150 keV and 800 keV, the probability for multiple scattering in one pixel ranges between 10 % and 20 % only. Then, to mitigate neutron scattering from structural materials around the detector as much as possible, the active crystal was embedded in a low-mass end-cap, arranged in such a way that the collimated neutron beam only hits the germanium crystal and two end-cap walls. These are made of very thin Be windows (300 µm) on two opposite sides. These windows allow for low-energy calibrations with e.g. a Fe-55 source. A typical energy spectrum of a Fe-55 source, emitting two low energy X-rays at 5.90 keV and 6.49 keV is shown in Fig. 3a. It demonstrates the excellent energy resolution achieved (FWHM of 135 eV ee at 5.90 keV), with a pulser resolution of about 80 eV ee . The Ge crystal is cooled with an electrical cryocooling system at the nominal temperature of 88 K.
Over the whole measurement campaign, the energy scale stability of the HPGe target detector was monitored via regular calibrations with a Fe-55 source. The peak position was stable below the 0.1 % level during the main measurement campaign as shown in Fig. 3b. For the physics runs used in the analysis (indicated within dashed lines), the RMS of the 9 calibration points was used for the uncertainty on the peak positions. After the measurement, the setup was dismantled and the detector was moved to another room to continue acquiring data to study the background produced by activation during neutron exposure. This change in experimental conditions is responsible for the small shift observed in the Fe-55 peak position in the last date cf. Fig. 3b. Therefore, this post-irradiation measurement was considered as an independent data set and only served as a validation of the calibration procedure. The energy scale was determined for each pixel independently using the information of four lines: the 5.90 keV and 6.49 keV lines from the Fe-55, the K α X-ray line at 9.87 keV which can originate from fluorescent X-rays in Ge escaping from the adjacent pixels and the 10.37 keV line from the deexcitation following electron capture (EC) in 68,71 Ge, produced in the crystal by cosmic ray activation and neutron beam irradiation. The two latter intrinsic lines were observed with sufficient statistics at a rate of about a few hundred counts per day in the 10 g crystal mass. They are visible in the background spectrum shown in Fig. 4. A linear calibration allowing a non-zero offset was employed, typically varying between -80 and 25 ADC bins, and a slope of about 2600 ADC bins/ keV ee . As expected the linear fit model allows a good description of the data for all four pixels, with p-values between 0.18 and 0.95. Depending on the pixel, the statistical uncertainties of the free offsets range in (9-12) eV, with a mean value over all pixels of 10.6 eV. The L-shell line at 1.30 keV, which originates mainly from atomic deexcitation after EC in 68,71 Ge, was used only as cross-check of the calibration procedure. For this purpose, the whole background data were merged since the line was barely visible in the independent data sets and pixels because of lacking statistics. The corresponding spectrum is shown in Fig. 4. With a step-like shape describing the background, the peak position is consistent with the expected value from literature within a statistical uncertainty of ± 10 eV ee . Changes of the binning or of the fit range do not alter this agreement, with deviations of the mean peak position staying within ± 10 eV ee . When using a simpler and less probable description of the background consisting in a linear model, a maximal deviation from the expected value of 12 eV ee was observed. Based on the uncertainty derived from the linear calibration fits and on this cross-check, an absolute deviation of ± 12 eV ee was therefore adopted as a conservative systematic uncertainty on the energy scale.
From calibration data, the energy resolution of the detector was modeled as: where σ 0 coincides with the pulser resolution, F = 0.13 is the material specific Fano factor and its value was chosen such to reproduce the difference between the pulser and the Fe-55 line resolutions. For ε e−h = 2.96 eV, the energy required to create an electron-hole pair at 90 K [29], this parametrization gives a good description of the experimental data in the region of interest up to ∼ 10 keV ee . The temporal stability of the energy resolution was carefully monitored throughout the experiment. HPGe detectors are indeed known to be sensitive to neutron damages. Because of incomplete charge collection, a deterioration of the resolution and the appearance of tails on the low-energy side of peaks might be observed. A critical threshold neutron fluence of about 10 9 cm −2 is often reported in the literature for n-type HPGe detectors, which are expected to be less sensitive to damages than p-type detectors [30][31][32]. However, this has only been taken as a rough indication since the effect is often considered for the irradiation of large kg-sized coaxial detectors with fast neutrons (see e.g. [33]) and a dependence on the crystal temperature also seems to exist [34]. To be on the safe side, the integrated neutron fluence was therefore limited to 6×10 7 cm −2 over the whole experimental campaign. The energy resolution of the 5.90 keV reference line remained stable within ± 3 eV ee during the entire measurement campaign. Moreover, the 59.5 keV line from regular measurements with an Am-241 source was used to monitor the energy resolution at higher energy, where effects due to incomplete charge collection become more visible due to the dominance of the second term in Eq. (1). The width of the line remained stable with a resolution (FWHM) of 416 eV ee at 59.5 keV and did not show any indication of a degradation of the charge collection during and after the neutron irradiation.
Since no photon sources are available for calibration purpose between the threshold and 1.30 keV ee , detailed investigations of the detector response were performed with artificial signals produced by pulse generators with the same risetime as for physical signals. Fine-grained scans allowed to precisely measure the trigger efficiency (cf. Sect. 2.5) for the four pixels as a function of the expected energy, derived from a careful calibration of the pulse generator amplitudes. As shown on the upper figure in Fig. 5, they exhibit a similar behaviour. A very good description of the experimental trigger efficiency curves was obtained via: where erf is the error function with typical parameters t 1 = 170 eV ee and t 2 = 65 eV ee . With the same pulser scans, the linearity of the DAQ chain in the sub-keV ee region was investigated. Deviations might appear at low energy in dependence of the reconstruction algorithms [35]. Such deviations from a purely linear behavior were observed below ∼ 400 eV ee and are shown in the bottom panel of Fig. 5. They were attributed to two nearly independent DAQ-related effects. First, as the amplitude of the signals decreases, artificial delays of the trigger time stamp up to a few µs were observed. This effect is due to the use of a specific trigger logic (cf. Sect. 2.5) and is observed in quenching data as well, see Fig. 9. Because of the relatively short shaping time of 4 µs, these delays affect the reconstructed energy value obtained by trapezoidal shaping. The corresponding negative non-linearity effect was less than 5 eV ee at 400 eV ee and reaches its maximum of about 25 eV ee at 250 eV ee . Second, the drop of the trigger efficiency for 250 eV ee is responsible for an artificial positive deviation from linearity, overwhelming the above-mentioned effect. This will be further discussed in the analysis section (Sect. 3.2) and is illustrated in Fig. 10.
Being precisely measured, all above mentioned effects are corrected for in the quenching analysis (Sect. 3.2) by taking them into account in the form of a detector response matrix. The energy resolution is implemented following Eq. (1). For the trigger efficiency, the analytical description of Eq. (2) is used and the non-linearies are described by the black model shown in Fig. 5. This matrix is then used to compute the expected energy distributions.

Neutron detectors
For the detection of the scattered neutrons, eleven LS detectors were designed and assembled at MPIK. They consist of cylindrical PTFE cells with very thin front walls.
Diameter and height of the cells are 5 and 6 cm, respectively, resulting in an active volume of about 120 cm 3 . The cells were filled with EJ-301(Eljen Technology) LS, which is identical to the well-known NE213 scintillator [36] and was chosen for its fast response, high light yield and good pulse shape discrimination (PSD) properties. The cells are optically coupled to 2" photo-multiplier tubes (PMTs). Detector cells and PMTs are embedded inside a lighttight aluminium housing. Neutrons are detected via the scintillation light produced from the prompt recoil signal, providing a fast response (∼ ns) adapted for the timing requirements of the coincidence measurement.
The neutron response of these detectors was characterized in detail in a dedicated commissioning campaign to ensure that all the detectors were equivalent and their data can be combined. The measurements were carried out in several monoenergetic neutron fields produced in open geometry using the 7 Li(p,n) 7 Be reaction and a 70 µg/cm 2 metallic lithium target and ns-pulsed proton beams (cf. Tab. 1). For neutron emission angles of 0°, the neutron fluence at the position of the detector cells was measured using a de Pangher [37] long counter which is the reference instrument for routine fluence measurements at the PTB. For non-zero emission angles, the neutron fluence was calculated from the measured 0°value using the known angular dependence of the neutron yield [38,39]. Nonmonoenergetic neutrons of lower energy resulting from neutron scattering in the production target were discriminated using the time-of-flight (TOF) technique. Residual background resulting from forward scattering neutrons from air and structural materials was determined in a separate measurement with a shadow cone placed between the production target and the detectors. The trigger threshold was calibrated in electron-equivalent energy using Cs-137 and Ba-133 photon sources emitting γ-rays in the range from 80 keV to 667 keV. For a trigger threshold of 12 keV ee , the neutron detection efficiency of the detectors was about 75 % over the energy range of interest from 250 keV to 800 keV with a slight variation of about 5 %. Neutron-and photon-induced events were separated by measuring the signal decay times using the ratio of the signals integrated over a short and a long gate period. Representative PSD distributions for two different neutron beam energies are shown in Fig. 6. A good discrimination between proton recoils (induced by neutrons) and electron recoils (ambient γ-rays) is obtained and allows to get rid of about 80 % of the ambient γ-ray background while keeping a neutron efficiency of 85-95 % depending on the beam energy.
Based on these commissioning data, the small discrepancies observed between the eleven detectors and variations over time were quantified in terms of relative variations in the neutron detection efficiency. The impact of small shape discrepancies in the neutron response distribution between the detectors -shown in Fig. 7 -was mitigated to 3 % by choosing an energy threshold of 250 ADC bins corresponding to 12 keV ee , indicated by a dashed line. During the measurement campaign, the stability of the LS detectors, including scintillator liquid, high- voltage, PMT and trigger thresholds, was regularly monitored using the Cs-137 and Ba-133 photon sources. The small gain variations between the detectors amount to less than 1.5 %. Finally, time variations in the response are responsible for an additional 1.5 % effect.
The detectors were positioned in a circular array at about 45 cm away from the HPGe target with a better coverage for the smallest angles -corresponding to the lowest recoil energies. The nuclear recoil energy E nr deposited in the germanium target of atomic mass A for a neutron scattering angle θ is: where E 0 is the initial neutron energy and θ is the scattering angle of the neutron as indicated in Fig. 1. Scattering angles between 17 o and 45 o were selected, which allowed to probe nuclear recoils between 0.4 and 6.3 keV nr , as illustrated in Fig. 8.

DAQ
A DAQ system based on commercially available CAEN electronic modules was used to acquire data from the LS and the HPGe detectors as well as the beam information. In order to fulfill both the requirements of a TOF  Fig. 7: Energy distribution obtained for four representative LS detectors exposed to a 500 keV neutron beam during commissioning. The bottom distributions are obtained by shadowing the LiF target with a polyethylene cone and allow to estimate the background contribution from air scattering. The threshold at 250 ADC used for the analysis (black dashed line) was chosen such to minimize the discrepancies in terms of neutron efficiency between the detectors. analysis with fast PMT signals and of the slower shaping times used for the HPGe detector signals, a modular configuration was developed from digitizers offering pulse processing implemented at the FPGA level.
The HPGe signals were sampled at a rate of 100 MHz rate by a V1782 multi-channel analyzer. A combination of a slow and a fast triangular discriminator was used for the trigger. This gave better results in terms of detection efficiency of small signals compared to the standard algorithm. The fast and slow triangular discriminators had shaping times of 0.8 µs and 2.4 µs, respectively. The energy of each event was reconstructed by a trapezoidal shaping filter with a shaping time of 4 µs, optimized in terms of energy resolution.
The signals from the PMTs of the LS detectors were acquired with a V1725 module at a sampling rate of 250 MHz. The charge integration windows were chosen to maximize the particle discrimination by PSD. A time resolution of 6 ns (FWHM) for physical LS signals was achieved.
The two modules shared a common synchronization clock, allowing to identify coincidences between LS and HPGe events. For practical reasons, the reference signal provided by the proton beam pick-up was sampled by both modules and was saved only in the presence of a physical trigger (LS or HPGe).
All triggering events were saved and the coincidence selection between the signals from the HPGe and the LS detectors was performed offline. Typical rates encountered during a quenching data run were in the range of (250-550) Hz for the LS detectors. For the HPGe target, typical rates were always dominated by the noise trigger rate, ranging from 60 Hz to 1.2 kHz depending on the pixels.

Coincidence data selection
The data selection for the analysis relies on a three-fold time coincidence between the proton beam pick-up signal, the signals from the HPGe target and the LS detectors. In practice, this was achieved by first selecting the neutroninduced events in the LS detectors using their TOF and PSD distributions. A selection window of 20 ns was chosen for the TOF distribution. The size of this window was determined from the LS data without coincidence requirements and corresponds roughly to a 3 σ wide time window around the neutron arrival time. This selection was extended to 40 ns for the 250 keV data set. Neutron signals were discriminated from the ambient γ-ray background with an additional PSD selection. In order to keep a constant neutron detection efficiency, the PSD cut was derived individually for each LS detector from a reference calibration measurement with an AmBe neutron source. Events with a PSD higher than q cut = q n − 3σ n were selected, where q n and σ n are the mean peak position and standard deviations of the nuclear recoil PSD distributions. The time window for a coincidence between a HPGe event and a LS event was fixed to 1.6 µs to take into account the increased smearing of the trigger time stamp at low energies detailed in Sect. 2.3. Accidental background essentially consists of random noise. Its contribution was evaluated by opening 10 random coincidence gates outside of the main coincidence window. This selection is illustrated in Fig. 9. In this way, a coincidence rate of the order of 100 counts per hour and per HPGe pixel was obtained for each angle. The coincidence energy distributions of interest for the quenching analysis were obtained by adding up the contributions of the four pixels.

Analysis procedure
The expected nuclear recoil energy distribution for each angle is obtained by running a MC simulation taking into account the geometrical extension of the detectors, the mean free paths of the neutrons inside the LS and the width of the neutron beam. The mean nuclear recoil E nr for each detector position is taken as the mean of the distributions. The widths of the expected distributions were found to be systematically smaller than the ones observed in the experimental quenching data. This effect is further discussed at the end of this section. Due to the lack of a quantitative description of this spread, a conservative approach was adopted and the nuclear recoil energy distributions were simply modeled by a Gaussian with mean E nr derived from the MC distribution and a free width σ nr (in keV nr ). An uncertainty on E nr is estimated for each angle from the propagation of the experimental uncertainties of the spatial coordinates measurement on site. An angular uncertainty of 0.5 o to 1.5 o was obtained, which translates into a 0.05 keV nr to 0.2 keV nr uncertainty in nuclear recoil energy, depending on the angle and on the beam energy. These nuclear recoil energy distributions are convoluted with the response matrix of the HPGe detector to obtain the model M ee -including the quenching process -used to describe the data. On top of resolution effects, it is of primary importance to take into account the detector response at low energies detailed in Sect. 2.3. In particular, detected distributions in the region below 1 keV nr are expected to be artificially shifted towards higher energies due to the drop of trigger efficiency. If the effect is not properly considered, this may lead to a biased value of the ionization quenching factor as illustrated in Fig. 10.
In order to normalize the models M ee to the data, a common signal coincidence rate n d 0 per unit of beam The "detected" distributions are obtained by convolving the "true" distributions with the HPGe detector response, including the energy resolution and the effects of trigger efficiency and energy scale non-linearities discussed in Sect. 2.3. In particular, if the trigger efficiency is not taken into account properly, this can lead to apparent nonlinearities in the energy scale and biased results in terms of ionization quenching factor. Due to the poor discrimination between the detected shapes, rate and ionization quenching factor strongly anti-correlate at these energies. This is illustrated by the blue 1 σ contour in the inset plot obtained when fitting the low energy coincidence data distribution without any norm constraint. The addition of the norm constraints summarized in Tab. 2 clearly provides a sensitivity improvement and allows to reduce the statistical uncertainties on q by about a factor 2 as shown by the corresponding orange contours. charge was determined for each data set. Indeed, within a given data set d corresponding to a given neutron beam energy, the coincidence rate obtained for each angle is expected to reflect only the differential angular elastic neutron scattering cross-section in Ge. It was cross-checked that the correlated counts corrected by the neutron elastic scattering cross-section in Ge and the acquisition time for each angle were constant within a data set for the energies above the region of trigger efficiency loss. The common rate n d 0 was taken as the mean of the obtained values.
For the modeling of the background, the accidental component is directly estimated from data as described in the previous paragraph. An additional continuous correlated component arising from multiple scattering in the crystal and the surrounding materials populates the correlated energy distributions up to a few keV. It amounts to about 20 % of the total correlated rate between 0 and 12-13 keV ee (end of the dynamic range). For the 500, 640 and 800 keV datasets, this background contributes to (1 − 4) % to the counting rate in the region of interest defined as a 2 σ window around the signal peak. It is higher for the 250 keV dataset with a contribution at the level of ∼ 10 %. source of uncertainty cross-section 5.0 % LS stability 1.5 % deviation between LS detectors 1.5 % shape variations between LS detectors 3.0 % total 6.2 % It has therefore been modeled for each dataset d by a simple flat contribution of common amplitude b d 0 obtained by a combined fit of the coincidence distributions above the peak region. Two refined modelisations of this background were also tested as extreme cases: the first one with a quadratic increase extracted from the integrated recoil distributions (see Fig. 15 in Sect. 4.2) and the second one with the same model convoluted by the trigger efficiency. The impact on the result was found to be negligible.
The coincidence energy distributions for each angle i were treated separately with a model-independent ionization quenching factor q i -assumed to be constant in energy -by minimizing the following χ 2 function: where the index j is running over the bins of the coincidence energy distribution i and d j are the corresponding counting rates with their associated statistical uncertainties σ j . Apart from the ionization quenching factor q i and the width σ nr, i (expressed in units of keV nr ), which are left free in the fit, the expected model M ee, j depends on all detector model ingredients described in Sect. 2.3 such as trigger efficiency and energy non-linearities. It is scaled using the integrated beam current c i and f θi , a dimensionless factor accounting for the differential angular elastic neutron scattering cross-section. Also included are the contributions from the accidental coincidences (a i j ) and the correlated background (b d 0 ). The nuisance parameter n d i , driving the normalization of the model, is constrained in the fit via the addition of a Gaussian pull-term incorporating the normalization knowledge: n d i can vary around its nominal value n d 0 , with a tolerance driven by Φ d , coming from both the characterized differences in terms of neutron detection efficiency of the LS and the knowledge on the cross-section. For the latter, the following evaluated nuclear data libraries were considered: ENDF/B-VII [40], JEFF-3.3 [41], JENDL-4.0 [42] and CENDL-3.1 [43]. The corresponding uncertainty on the expected counting rate was taken from the dispersion between the different evaluated databases. These normalization constraints are sum- Fig. 11: Best fit coincidence rates (points) for each data set and corresponding expectation (square) calculated from the global norm n d 0 , the neutron scattering cross section and the trigger efficiency. Rates are given per hour. Differences in current and neutron yield of the Li(p,n) reaction as a function of the energy explain the overall rate difference between the data sets. Coincidence energy distributions obtained in the HPGe detector after coincidence selection with three LS detectors at different angles for the data set with the neutron beam energy E n,0 = 500 keV. The accidental contribution, here in gray, is determined by opening random windows and consists of electronic noise. The best-fit models are superimposed in dashed lines. marized in Tab. 2. They introduce sensitivity to the particularly delicate low energy region where, due to the steep trigger efficiency, the ionization quenching factor degenerates with the rate of detected events. The resulting anticorrelation is shown in the inset in Fig. 10. For data below 1 keV nr , the addition of the constraint via the pull-term allows to reduce the statistical uncertainties by a factor of 1.5 to 2 with respect to an analysis without any rate constraint.
Best-fit values of the coincidence rate for each quenching data set along with their expected values are shown in Fig. 11. They are in a very good agreement, confirming a good understanding of all involved effects and systematics. The decreasing trend above 2 keV nr reflects the angu-lar dependency of the scattering cross-section, whereas the drop for 1.5 keV nr is due to the loss of trigger efficiency and shows therefore the same behavior for all datasets. It becomes also clear that with counting rates of a few hundred counts per hour, the statistics is sufficient for energies above 1-2 keV nr , whereas due to trigger efficiency effects, the available statistics notably drops below 1 keV nr . Fig. 12 shows coincidence distributions for three exemplary scattering angles at fixed beam energy along with their best-fit model.
Regarding the width of the nuclear recoil distributions, best-fit values of σ nr exceed the MC predictions by (30−60) %. Possible explanations are an additional smearing of the quenching factor due to the stochastic nature of ion collisions or a mismodeling of the energy spread of the neutron beam. However, the latter case seems to be disfavored: the estimated proton energy losses inside the Li target have been compared with SRIM calculations [44] and validated by previous measurements at the PIAF facility. Moreover, MC simulations show that the contribution of scattered neutrons at the exit of the collimator is at the percent level, which cannot explain the large discrepancy observed. For a better quantification of the additional broadening and for a better comparison with other experiments, the following setup-independent quantities have been defined: where σ nr, M C is the width predicted from the MC simulation and σ nr, exp are the ones derived from the data, i.e. the values of σ nr,i extracted from the fit of the experimental energy distributions (cf. Eq. 4). The quantity σ nr, missing was found to scale roughly with √ E nr and amounts to ∼ 0.3 keV nr at 2 keV nr and 0.55 keV nr at 5 keV nr . The corresponding values of the broadening factor B agree with the ones measured in [45], another recent Ge-based setup.
The coincidence energy distribution corresponding to the lowest nuclear recoil energy was obtained with the lowest beam energy E n,0 = 250 keV and corresponds to an expected nuclear recoil energy of 0.4 keV nr . Because of the large statistical uncertainty, a conservative approach was favored and an upper limit in terms of the ionization quenching factor was derived from a likelihood analysis, making use of the normalization constraints from the 250 keV data set. The upper limit at 90 % C.L. was determined from the expected likelihood distributions obtained via a toy MC. The energy distributions of the 0.4 keV nr coincidence data and the model for the corresponding estimated upper limit for the quenching factor are shown in Fig. 13. The same procedure was applied for the coincidence data at 0.6 keV nr from the same beam energy data set. Both limits are included in Fig. 14.

Combination of data
The ionization quenching factor obtained from the fits for each coincidence distribution are provided in Tab. 5a, 6a, 5b and 6b for the data sets with E n,0 = 250, 500, 640 and 800 keV, respectively. In the last column the total uncorrelated uncertainties are reported, including both the statistical and the systematic uncertainties from the spatial coordinates measurement. For most of the cases, the latter dominates. For each beam energy, the same angles were measured several times in each hemisphere and LS detectors were switched in order to cross-check for systematic effects related to detector positioning and neutron detection. No significant discrepancies were found. Therefore, the corresponding data were merged for better visibility and are illustrated in Fig. 14.
The data points were combined via a fit within the semi-empirical Lindhard theory [7] which describes the energy dependence of the ionization quenching factor q(E nr ) with a single free parameter k: where E nr is the recoil energy in keV and, for germanium, ε = 11.5 Z −7/3 E nr and g(ε) = 3 ε 0.15 + 0.7 ε 0.6 + ε. This model allowed for a good description of the data and did not require the inclusion of modified theories with additional parameters, such as the adiabatic correction proposed in [12].
The complementarity between the different beam energies was exploited to access the same nuclear recoils energies with a different combination of systematic uncertainties. From Fig. 8, it can for instance be seen that the expected width of the coincidence energy distributions  Table 3: Best-fit of quenching data separated into different data sets (divided into beam energies and pixels).
strongly varies between the different energies. Moreover, a bias in the spatial measurement of the scattering angles would impact the results obtained from the different beam energies in a different manner. The systematic uncertainty for a given neutron beam energy E n estimated via TOF measurement is fully correlated between the corresponding data points. As reported on Tab. 3, I and II, best-fit values of k are consistent for the different beam energies.
Since the fit is driven by the higher energy points with the smallest uncertainties, best-fit values were also compared separately in a restricted low energy range below 1.8 keV nr and in a high energy range above 1.8 keV nr . Best-fit values of k are consistent for the two energy ranges as reported in Tab. 3, I vs. II, for the different data sets.
The uncertainties related to the HPGe detector (energy scale, detector response) are treated as fully correlated between all data points. The systematic uncertainty on the energy scale was derived by allowing for a constant shift ∆E ee = ± 12 eV ee in the reconstructed ionization energy, implying a larger effect towards lower energies. The value comes from the constraint on the energy scale obtained with the 1.30 keV ee L-shell atomic deexcitation line following Ge EC decays, as discussed in Sect. 2.3. In addition, to quantify the impact of a potential mismodeling of the detector response close to the threshold, the mean position of the trigger efficiency curve was shifted by ± 10 eV ee which is a conservative upper limit for the four pixels. The impact is negligible above 1.2 keV nr . Due to the anti-correlation between detected rate and ionization quenching factor (discussed in Sect. 3.2), the impact grows rapidly below this energy. It becomes the dominant contribution for the lowest energy points at E nr 0.8 keV nr , as illustrated by the enlarged uncertainty band in Fig. 14.
Advantage was taken of the pixelized structure of the HPGe detector to validate these estimations. Indeed, as they were individually calibrated and characterized, the four pixels can be considered as almost independent detectors. Moreover, different acquisition conditions in terms of noise rate were used on purpose such that the acciden-  Tab. 4 summarizes all the above mentioned sources of systematic uncertainties and their contributions to the final results. For E nr 2 keV nr , the statistical uncertainty is negligible and the correlated uncertainty on the energy scale of the HPGe detector is the major contributor to the total quoted uncertainty, followed by the systematic uncertainty on the beam energy and of the geometry of the setup. In the low energy range, statistical uncertainties are not negligible anymore. The correlated uncertainty on the energy scale is still the dominant contribution. Furthermore, as emphasized by the enlarged uncertainty band in Fig. 14, the uncertainty due to the modeling of the trigger efficiency dominates the sub-keV region. The large contribution of the uncertainty of the energy scale is related to the lack of photon sources in the sub-keV ee region which could be used for calibration. A smaller statistical uncertainty from the 1.30 keV ee activation line would have allowed to put a more stringent constraint on the energy scale. The present approach is therefore conservative and does also cover the uncertainties related to the small nonlinearities observed in the few hundreds of eV ee region, as illustrated by the gray band in the bottom panel of Fig. 5. For the sub-keV nr region, the precision of the measurement is intrinsically limited by the trigger efficiency of the HPGe detector. A precise knowledge of the detector response at these energies allows partly to overcome this limitation but HPGe detectors with thresholds below 100 eV ee are needed to access the sub-keV nr region with a precision better than ∼ 0.01 on the ionization quenching factor.
To summarize, the combined fit of all data points between 0.6 keV nr and 6.3 keV nr using the Lindhard theory yields a best fit value of k = 0.162 with a total uncertainty (stat+syst) of 0.004. It is represented by the black line in Fig. 14.

Analysis of the integral energy distribution of recoiling nuclei
In this last section an additional cross-check is proposed, considered as almost independent of the main result of this article. A complementary approach to determine ionization quenching factors consists in comparing the energy distribution of recoil nuclei integrated over all scattering angles (i.e. without coincidence requirements) with a MC simulation of the experiment. This approach was for instance used in [16,18,24]. Such an analysis relies on an accurate modeling of the setup and a quenching model. Its parameters can be tuned in the simulation in order to reproduce the observed data. Although we strongly favor the direct and model-independent technique of a coincidence measurement, advantage was taken of the high-statistics neutron recoil data integrated over all angles and for different neutron beam energies in the HPGe detector collected as a by-product during data collection. Note that the range of recoil energies probed in this way is much broader than the one studied by selecting only small scattering angles: from Eq. (3), the maximum recoil energy obtained for back-scattered neutrons equals e.g. 26 keV nr for E n,0 = 500 keV. Integrated energy distributions of recoil  Table 4: Summary of the uncertainties sources taken into account in the analysis. Statistics and geometry uncertainties are uncorrelated between the points, whereas beam energy uncertainties are correlated for each data set. HPGe detector modeling effects are fully correlated. Their impact in terms of absolute uncertainty on the quenching parameter k of the Lindhard theory is given in the two last columns, for two different energy ranges. nuclei without coincidence requirements were compared to a Geant4 [46] simulation of the setup including the detailed geometry of the detector end-cap. A simplified description of the neutron beam was implemented: the beam profile was inferred from TOF measurements at different distances and the energy profile was taken from calculations accounting for the thickness of the lithium [39]. As for the main analysis, the expected nuclear recoil distribution was convoluted with the detector response M ee as detailed in Sect. 2.3. Since this analysis is intended to be a cross-check of the main coincidence analysis, only the ionization quenching factor description within the Lindhard theory with the best-fit k = 0.162 was implemented. An overall reasonable agreement was found for the four data sets, especially concerning the position of the endpoint in energy, as shown in Fig. 15. Detector resolution effects -described by Eq. (1) and the constraint of the Fe-55 calibration at these energies -cannot explain the steeper edge of the end-point found in the simulation. An additional energy broadening resulting from the collima-tor could be the origin of the discrepancy. The impact of surrounding materials and of the beam profile were investigated and found to play only a secondary role for the spectral shape of the energy distribution. However, it is worth noting that the shape of the energy distribution strongly depends on the evaluated nuclear database. In particular, a much better agreement was found for the region below 1 keV nr when making use of JEFF 3.3 [41] instead of ENDF-BVII.1 [40]. This dependence on the databases in the MC simulation illustrates the difficulty of extracting quenching factors out of integrated distributions and it reinforces the need of kinematically constrained measurements like the one presented in this work.
A steep increase of the ionization quenching factor below 1 keV nr as measured in [18] was also implemented. The resulting energy distributions of recoil nuclei for the two evaluated nuclear databases are shown by the dashed lines in Fig. 15. The addition of such an enhancement only affects the lowest part of the energy distribution and does not significantly improve the overall description of the data.

Conclusion
A direct and precise measurement of the ionization quenching factor for nuclear recoils in germanium was presented. A HPGe detector with a mass of 10 g and a thickness of 6 mm was operated at 88 K and exposed to monoenergetic neutron beams with energies ranging from 250 keV to 800 keV at the accelerator facility PIAF of PTB in Braunschweig (Germany). Nuclear recoil energies between 0.4 keV nr and 6.3 keV nr were selected by the detection in coincidence of the scattered neutrons in LS detectors. In this energy range, data are compatible with the Lindhard theory prediction with k = 0.162 ± 0.004 (stat+syst), as presented in Fig. 14.
The uncertainties of the measurement were discussed in detail. Their contributions are summarized in Tab. 4. Special attention was paid to constrain them by performing multiple cross-checks, e.g. by using different monoenergetic neutron beams, by interchanging the neutron coincidence detectors and by taking advantage of the segmented structure of the target detector. In the region of interest of the current experiments looking for CEνNS at nuclear power reactors with expected nuclear recoil energies above 1.5 keV nr , a measurement of the quenching factor with a relative uncertainty of a few percent is provided. This was achieved by determining the scattering angles with uncertainties of about one degree and by precisely monitoring the energy of the neutron beam. For recoil energies of less than ∼ 1 keV nr , the fully correlated systematic uncertainties related to HPGe detector response modeling dominate. In particular, energy threshold effects -and therefore decreasing statistics -illustrated the difficulty of extracting ionization quenching factors at these energies even with improved state-of-the-art HPGe detectors.
Thanks to a precise characterization of the detector response and a very good control on the coincidence rates, this limitation could partly be overcome and quenching factors were extracted with a precision of (5-10) % from data in this low energy region.
As an additional independent cross-check, a comparison of the integrated recoil energy spectrum (i.e. collected in the HPGe without coincidence requirement) with a MC simulation of the setup was performed and an overall good agreement was found when making use of the JEFF 3.3 evaluated nuclear database.
This work significantly contributes to the understanding of the ionization quenching factor of keV nuclear recoils in Ge. Above 1 keV nr , the quenching factor was measured with a precision of a few percent. The data follow the Lindhard theory and are in agreement with several previous measurements that were performed at liquid nitrogen temperature in the 0.5-10 keV nr range [20][21][22][23]. However, this result is found to fully disagree with a recent precision measurement at 50 mK [19] in which a significant drop of the quenching factor below 7 keV nr was found. In this case the measurement technique and experimental conditions (temperature, electric field) were nevertheless significantly different from the ones reported in this article. Below 1 keV nr , our data do not confirm the outcome of another recent measurement in Ge at 77 K that suggests an enhanced quenching factor with respect to the Lindhard prediction [18]. Our analysis does include a refined description of non-linearities and trigger efficiency effects affecting this low energy regime, refuting in particular the claims raised in [47] and the proposed correction to obtain a better agreement with the result found in [18]. However, resolving these discrepancies is essential for the present and future experiments based on HPGe detectors in the search for coherent elastic-neutrino nucleus scattering and light dark matter.   Table 5: Analysis results for the runs 3, 4 (a) and 5, 6 (b). For each point, the irradiation time t is indicated and the selected angle θ is given with its geometrical uncertainty δ θ . The corresponding selected mean nuclear recoil energy E nr , its uncertainty δE nr and its width σ nr are obtained from a MC simulation of the setup. The χ 2 /ndf and pvalues of the fit are reported along with the deviation Φ norm in sigma units of the normalization nuisance parameter, i.e. Φ norm = (n d i,BF − n d 0 )/Φ d following the notations of Eq. (4). For each distribution, the best-fit of the energy independent quenching factor is given by q. Its statistical uncertainty is denoted by δq stat , whereas δq includes the geometrical uncertainty. The dashed lines indicate how the points were combined and shown on Fig. 14