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

Abstract

This article reports the measurement of the ionization quenching factor in germanium for nuclear recoil energies in the keV range. Precise knowledge of this factor in this energy range is highly 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 using 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 keV$_{nr}$ to 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).

1 Introduction

Recent advances in the detection of very low energy signals opened up new possibilities to study coherent elastic neutrino-nucleus scattering (CE$\nu$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$\nu$NS in HPGe detectors [8,9,10,11], recoils of 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$\nu$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 $\sim$ 250 eV$_{ee}$, experiments like CONUS will be able to detect the $\sim$ 1.5 keV$_{nr}$ recoils induced via CE$\nu$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$\nu$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.

2 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. Figure 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 (⁶LiGl) and a barium fluoride scintillation detector (BaF₂) 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 Table 1. Prior to the actual measurement of the quenching factor (runs 1–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).

Fig. 1

Photo and scheme (not at scale) of the experiment installed at the PTB facility. 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

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 $\delta E_p$. The nominal lithium mass per unit area used for this calculation is $m_{Li}$

2.2 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 neutron yield per unit target thickness and the reduced yield of parasitic high-energy photons from $^{19}\mathrm {F(p,}\alpha \gamma )^{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 $\upmu$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 ⁷Li(p,n)⁷Be reaction produces only monoenergetic neutrons. For higher proton energies, transitions to the first excited state in the residual ⁷Be nucleus result in a second neutron contribution of lower energy. The energy of the proton beam was measured using a 90$^{\circ }$ 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 ⁶LiGl detector and a fast BaF₂ detector. The ⁶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 time is less than 0.7 ns. In this detector, neutrons are detected via the $^6\mathrm {Li(n,t)}^4\mathrm {He}$ reaction. The large Q-value of 2468 keV facilitates a separation of neutron- and photon-induced events by a simple pulse-height threshold. Figure 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 ⁶LiGl detector is about 4 ns, i.e. larger than the duration of the proton beam pulses. Therefore, the fast BaF₂ 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 ⁶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.

Fig. 2

Pulse charge versus time-of-flight for the ⁶LiGl detector placed in the beam line for monitoring of the beam energy. Neutron signals can be distinguished by their high charge. In this way, they are discriminated from the ambient $\gamma$-ray background associated to a low charge

Fig. 3

a Energy distribution in ADC bins for Pixel-1 from a calibration measurement with a Fe-55 source. The red line is the fit of the spectrum. The two peaks correspond to X-ray energies of 5.90 and 6.49 keV. b Relative variation of the mean peak position of the Fe-55 calibration line at 5.90 keV over the duration of the measurement campaign, relative to one of the measurements done on October 1, 2020

The characteristics of all neutron beams used for the present experiment are summarized in Table 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 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 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 $\delta 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–6, where the nominal and the experimental neutron energy agreed within their uncertainties. 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₂ 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 $\upmu$g cm$^{-2}$. A conservative estimate of $(790\,\pm \,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 ⁶LiGl detector.

2.3 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 $\times$ 9 mm$^2$. For this geometry and neutron energies between 150 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 $\upmu$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$_{\alpha }$ 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:

$$\begin{aligned} \sigma (E) = \sqrt{\sigma _0^2 + {\mathcal {F}}\cdot \varepsilon _{e-h}\cdot \text {E}} \end{aligned}$$

(1)

where $\sigma _0$ coincides with the pulser resolution, ${\mathcal {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 $\varepsilon _{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 $\sim$ 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 $\times$ 10$^7$ cm$^{-2}$ over the whole experimental campaign. The energy resolution of the 5.90 keV reference line remained stable within $\pm \,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.

Fig. 4

Calibrated background spectrum collected in the HPGe detector after beam irradiation: lines from activation are visible at 1.30, 9.87 and 10.37 keV$_{ee}$. The 1.30 keV$_{ee}$ line (red, in the inset) is used to validate the energy scale at low energy. The line position was tested with several background descriptions (blue), fit ranges and binnings

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:

$$\begin{aligned} \varepsilon _{trig} = 0.5\cdot \left( 1+\text {erf}\left( \frac{E_{ee}-t_1}{t_2}\right) \right) \end{aligned}$$

(2)

where erf is the error function with typical parameters $t_1$ = 170 eV$_{ee}$ and $t_2$ = 65 eV$_{ee}$.

Fig. 5

Top: measured trigger efficiencies for the four pixels (points) and best-fit model (dashed lines) using the analytical description of Eq. (2). Bottom: deviations from a purely linear energy scale measured with two different pulse generators (points) described by the black line. The ± 12 eV$_{ee}$ uncertainty band derived for the energy scale is shown in gray and it largely covers the small variations observed between the measurements and pixels. Both effects are included in the detector response matrix used to compute the expected energy distributions (cf. Sect. 3.2)

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 $\sim$ 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 $\mu$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 $\upmu$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 $\lesssim$ 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.

2.4 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³. 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 light-tight aluminium housing. Neutrons are detected via the scintillation light produced from the prompt recoil signal, providing a fast response ($\sim$ 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 ⁷Li(p,n)⁷Be reaction and a 70 $\upmu$g/cm² metallic lithium target and ns-pulsed proton beams (cf. Table 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]. Non-monoenergetic 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 $\gamma$-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 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 $\gamma$-rays) is obtained and allows to get rid of about 80% of the ambient $\gamma$-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.

Fig. 6

PSD ratio versus charge (in ADC units) for a LS detector placed in the neutron beam line during commissioning. The black line indicates the threshold used for the analysis. The neutrons, inducing proton recoils in the liquid, are identifiable by their high PSD ratio ($> rsim$ 0.1). The lower band consists of interactions from ambient $\gamma$-rays

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 $\theta$ is:

$$\begin{aligned} E_{nr} = \frac{2E_0}{(A+1)^2}\left( A+\sin ^2(\theta )-\cos (\theta )\cdot \sqrt{A^2-\sin ^2(\theta )}\right) \nonumber \\ \end{aligned}$$

(3)

where $E_0$ is the initial neutron energy and $\theta$ is the scattering angle of the neutron as indicated in Fig. 1. Scattering angles between 17$^{\text {o}}$ and 45$^{\text {o}}$ were selected, which allowed to probe nuclear recoils between 0.4 and 6.3 keV$_{nr}$, as illustrated in Fig. 8.

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

Fig. 8

Nuclear recoil energy obtained in germanium for incoming neutrons with energies 250–800 keV as a function of their scattering angle, according to Eq. (3). The expected angular spreads for the chosen setup are represented in the inset. The dots indicate the angles chosen for the measurement

2.5 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 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 $\upmu$s and 2.4 $\upmu$s, respectively. The energy of each event was reconstructed by a trapezoidal shaping filter with a shaping time of 4 $\upmu$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.

3 Analysis

3.1 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 neutron-induced 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 $\sigma$ 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 $\gamma$-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\sigma _n$ were selected, where q$_{n}$ and $\sigma _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 $\upmu$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.

Fig. 9

Data selection for the coincidence analysis based on the time difference between a LS and a HPGe signal for two different expected nuclear recoil distributions (red at 4.9 keV$_{nr}$ and gray at 0.6 keV$_{nr}$). The energy reconstructed in the HPGe detector is represented on the y-axis. Correlated signals are selected for $\varDelta$T$<1.6\,\upmu$s (blue line) enclosing the smearing of the distributions at low energy, whereas accidental coincidences are estimated in 10 windows starting from $\varDelta$T$>3.2\,\upmu$s (yellow line)

3.2 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 $\sigma _{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$^{{\circ }}$ to 1.5$^{{\circ }}$ 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 ${\mathcal {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.

Fig. 10

Illustration of trigger efficiency effects for an initial recoil energy of 0.8 keV$_{nr}$ under the assumption of an ionization quenching factor q of 0.16 (red) or 0.20 (blue). 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 non-linearities 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 $\sigma$ 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 Table 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

In order to normalize the models ${\mathcal {M}}_{ee}$ to the data, a common signal coincidence rate $n^d_0$ per unit of beam 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 $\sigma$ window around the signal peak. It is higher for the 250 keV dataset with a contribution at the level of $\sim$ 10%. It has therefore been modeled for each dataset d by a simple flat contribution of common amplitude $b_0^d$ 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 $\chi ^2$ function:

$$\begin{aligned} \chi ^2_i =&\sum _{j=0}^{N_{bins}} \frac{[d_j - n^d_i\,c_i\,f_{\theta _i}\,{\mathcal {M}}_{ee,j}(q_i, \sigma _{nr,i})- c_i\,b_{0}^d - a_j^i]^2}{\sigma _j^2} \nonumber \\&+ \frac{[n^d_i - n^d_0]^2}{\varPhi _d^2} \end{aligned}$$

(4)

Table 2 Constraints in terms of coincidence rate used in the fit. The constraint on the scattering cross-section reflects the discrepancies between the neutron databases. Constraints concerning the LS detectors are derived from the characterization detailed in Sect. 2.4

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 $\sigma _j$. Apart from the ionization quenching factor $q_i$ and the width $\sigma _{nr,\,i}$ (expressed in units of  keV$_{nr}$), which are left free in the fit, the expected model ${\mathcal {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_{\theta _i}$, a dimensionless factor accounting for the differential angular elastic neutron scattering cross-section. Also included are the contributions from the accidental coincidences ($a_j^i$) and the correlated background ($b_0^d$). 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 $\varPhi _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 summarized in Table 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 anti-correlation 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.

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

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 angular dependency of the scattering cross-section, whereas the drop for $\lesssim$ 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}$. Figure 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 $\sigma _{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:

$$\begin{aligned}&\sigma _{nr,\,missing}^2 \equiv \sigma _{nr,\,exp}^2~-~\sigma _{nr,\,MC}^2 \end{aligned}$$

(5)

$$\begin{aligned}&{\mathcal {B}} \equiv \sigma _{nr,\,missing}^2/E_{nr}^2 \end{aligned}$$

(6)

where $\sigma _{nr,\,MC}$ is the width predicted from the MC simulation and $\sigma _{nr,\,exp}$ are the ones derived from the data, i.e. the values of $\sigma _{nr,i}$ extracted from the fit of the experimental energy distributions (cf. Eq. 4). The quantity $\sigma _{nr,\,missing}$ was found to scale roughly with $\sqrt{E_{nr}}$ and amounts to $\sim$ 0.3 keV$_{nr}$ at 2 keV$_{nr}$ and 0.55 keV$_{nr}$ at 5 keV$_{nr}$. The corresponding values of the broadening factor ${\mathcal {B}}$ agree with the ones measured in [45], another recent Ge-based setup.

Fig. 12

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

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.

Fig. 13

Coincidence energy distribution obtained in the HPGe detector after coincidence selection E$_{nr}$ = 0.4 keV$_{nr}$ (blue) along with the estimated accidental contribution (gray). The upper limit at 90% C.L. corresponds to a quenching factor of 0.28 and the corresponding expectation is shown as dashed red line

4 Results and discussion

4.1 Combination of data

The ionization quenching factor obtained from the fits for each coincidence distribution are provided in Tables 5a, b and 6a, b 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.

Fig. 14

Ionization quenching factor as a function of the nuclear recoil energy. The data points are obtained for the four different data sets, each one corresponding to a different nominal neutron beam energy. Due to the low statistics for the beam energy at 250 keV only upper limits were extracted for the lowest energy data points. They are represented by the orange arrows. Indicated error bars are uncorrelated uncertainties (statistics and spatial coordinates measurements for each angle), whereas the correlated uncertainties (beam energy, energy scale and trigger efficiency) are represented by the blue band. The best-fit of these data points within a Lindhard theory description is obtained for $k\,=\,0.162\,\pm \,0.004$ (stat + sys) and is illustrated by the black curve

Table 3 Best-fit of quenching data separated into different data sets (divided into beam energies and pixels)

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

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:

$$\begin{aligned} E_{ee} = q(E_{nr})\cdot E_{nr} = \frac{k\,g(\varepsilon )}{1+k\,g(\varepsilon )}\cdot E_{nr} \end{aligned}$$

(7)

where $E_{nr}$ is the recoil energy in keV and, for germanium, $\varepsilon \,=\,11.5\,Z^{-7/3}E_{nr}$ and $g(\varepsilon )\,=\,3\,\varepsilon ^{0.15} + 0.7\,\varepsilon ^{0.6} + \varepsilon$. 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 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 Table 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 Table 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 $\varDelta E_{ee} = \pm \,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} \lesssim 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 accidental contribution varies by about one order of magnitude between the pixels. The full analysis was then repeated for each pixel independently. The results (reported in Table 3, III) are compatible with the results under I and II and do not indicate any underestimated source of uncertainty.

Table 4 summarizes all the above mentioned sources of systematic uncertainties and their contributions to the final results. For $E_{nr} > rsim$ 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 non-linearities 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 $\sim$ 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.

4.2 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 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 ${\mathcal {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.

Fig. 15

Measured integrated nuclear recoil distribution for E$_{n}$ = 249 keV (orange points) and E$_{n}$ = 490 keV (blue points) and associated Geant4 simulations (solid lines) with JEFF 3.3 (gray) and ENDF-BVII.1 (purple). Superimposed are the simulations including the extension of the quenching factor model found in [18] (dashed lines)

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 $\theta$ is given with its geometrical uncertainty $\delta _{\theta }$. The corresponding selected mean nuclear recoil energy $E_{nr}$, its uncertainty $\delta E_{nr}$ and its width $\sigma _{nr}$ are obtained from a MC simulation of the setup. The $\chi ^2$/ndf and p-values of the fit are reported along with the deviation $\varPhi _{norm}$ in sigma units of the normalization nuisance parameter, i.e. $\varPhi _{norm} = (n_{i, BF}^d-n_0^d)/\varPhi _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 $\delta q_{stat}$, whereas $\delta q$ includes the geometrical uncertainty. The dashed lines indicate how the points were combined and shown on Fig. 14

An overall reasonable agreement was found for the four data sets, especially concerning the position of the end-point 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 collimator 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.

Table 6 Analysis results for the runs 7, 8, 9 (a) and 10, 11 (b). The description of the columns is identical to the one in Tab. 5

4.3 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 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\,\pm \,0.004$ (stat + syst), as presented in Fig. 14.

The uncertainties of the measurement were discussed in detail. Their contributions are summarized in Table 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$\nu$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 $\sim$ 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.