Ionization efficiency at sub-keV energies for crystals and noble liquids

We study the ionization and light yields produced by nuclear recoils at low energies in pure crystals and noble liquids in the context of Lindhard’s integral equation, incorporating the effects of binding energy, improved modeling of the electronic stopping, and electronic straggling. We consider three different models for the electronic stopping power that incorporate Coulomb repulsion effects at low energies, and Bohr electronic stripping for high energies. Finally, we discuss possible new effects near threshold


Introduction
Dark Matter, Neutrinos and other rare searches [1] required to be detected by direct sub-keV recoil energy spectrum from terrestrial and astrophysical sources. For pure element ionization detectors like Si [2], Ge [3], liquid Argon (LAr) [4] and liquid Xenon (LXe) [5], the typical recoil energy analyzed is less than 1 keV, where at this regime only a fraction of this energy goes to ionization. The convertion of the total recoil energy to visible or electronic energy is given by the ionization efficiency or quenching factor. Lindhard integral equation with binding energy [6] already have proven to success in describing low energy measurments for ionization efficiency in silicon [7], this rediscovered approach can be used to obtained the ionization deposit energy that an non ionizable particle gives to a pure media when interacts with a single nuclei, by separating electronic from nuclear processes.
Although Lindhard integral equation 1 with constant binding energy describes accurately the data of silicon below 4 keV, it required an outranged value from the expected one, Frenkel pair energy [8]. Furthermore, other studies [9,10] evidenciate that the electronic stopping power computed by Lindhard [11] might be overestimated at low energies.
For DM or CEνNS searches with Si and Ge ionization detectors the ionization efficiency (quenching factor) plays an important role for calibration. Different quenching factor leads to different shape of energy spectrum specially near threshold (100 eV ee ). Quenching factor has the effect to move the events of the spectrum to low energy, where some of them can pass the threshold. This lead to a systematic error fluctuation in low energy detection experiments, where usually spectrum rate for different quenching factor models are reported, [12,13].

Ionization and Lindhard Integral Equation
When a non ionizing particle, like DM or a neutrino, interacts with a crystal detector it deposit a total recoil energy E R . When this happens part of the recoil energy is used to disrupt the atomic-binding U, so then the ion now moves with an kinetic energy E (small than E R ). This process continue until the energy of the ion is not enough to disrupt the atomic binding. In this work we are going to use Lindhard units ϵ = C E, where C = 11.5/Z 7/3 (1/keV). The energy is divided into atomic motion (ν) due to all subsequent ions collisions and ionization energy due to excite electrons (η). So we define the ionization efficiency or quenching factor by the quotient of the energy that goes to ionization over the total recoil energy deposited in the material, We are going to assume as work hypothesis that the ionization occurs only when the ion have the necessary energy to disrupt the atomic bonding of the crystal and moves freely. We claim that this argument is reasonable for energies up to 60 eV (recoil energies). The basic integral equation deduced by Lindhard [14] takes into account, in this case for atomic motionν; that the atomic motion energy given by the ion with energy E is equal to the atomic motion energy given by the ions after scattering, where one of them takes in to account ionization and the other target ion have to spend some energy in disrupting the binding,

008.2
Where T n is the nuclear recoil energy given by the incident ion to the ion at rest, i T ei is the energy transfer to electron during the ion colision, U is the atomic binding energy, σ n,e is the inelastic nuclear ion-ion corss section and is a fucntion of T n and i T ei andν e is the atomic motion due to electrons (negligible). As we can see at low energies the part of atomic motion given by electron can be neglected. Lindhard gives an asymptotic approximate solution to Eq.(2) that is just valid at high energies. 2 Lindhard approximate solution works fine in many cases, e.g. Si, Ge, etc, for energies above 10 keV, but already for Si this parametrization fails below 4 keV [7].

Simplify Integro Differential Equation
In order to handle the integral equation Eq.(2), we make the following basic assumptions (most of them made by Lindhard); I we consider that nuclear recoil energy and the energy given to electrons are small compared to the initial kinetic energy of the ion E, II effects of electronic and atomic collisions can be treated separately. In contrast to Lindhard we are not going to neglect the effects of binding energy.
Furthermore, in order to get a solvable integro-differential equation, we need to relax approximation I up to second order. Also we consider nuclear stopping power using the universal nuclear cross section, that can have some variations depending on the inter atomic potentials used. With this and other details we can deduce a second order integro differential equation, that can be solve by using the shooting method [6], where, t = ϵ 2 sin 2 (θ /2) (center of mass frame), σ n is the nuclear scattering cross section and, S e = kϵ 1/2 is electronic stopping power of the medium, if the electronic stopping is zero at all energies, then the quenching factor also might be zero. This equation predicts an energy threshold at 2u. By firt considering a constant binding energy model and the electronic stopping power given by Lindhard, we succes in describe Si and Ge measuremnts at low energies [6].

High energies effects
One of the limitations of the aforementioned approach in Si, relies in having a cut off too high compared to the expected threshold given by the energy to create a Frenkel-pair (≈ 30 eV) [8].
To affront this limitations, in this study we considered a varying binding energy model and Coulomb repulsion effects for low energies and electron stripping at higher energies. Hence, Lindhard electronic stopping is not valid at low energies (≈ 5 <keV). Also we are going to considered electronic straggling effects in to the cascade process, this is done by expanding at second order in i T ei the termν E − T n − i T ei in Eq. (2). The main effect of including electronic straggling is to low the quenching factor near threshold and slightly increase it at high energies. These effects can be added directly into Eq.(3).
For high energy the electronic stopping has to be corrected by velocity effects that reduce the effective number of electron in the ions, this can be modeled by using the Bohr stripping criteria [15]. The oscillatory effects for electronic stopping as function of the number of electrons of the incoming ion are very well known. Since due to Bohr stripping the effective number of electrons reduce, implies that for energies much lower than the Bragg peak (<20 Mev in Si) this effect may be active. This oscillations are caused by an appearance of a phase shift,Firedel sum rule [16], to maintain neutrality of electron Fermi gas, it can be tested that this effect produce a better agreement with the model and the available data for electronic stopping in silicon.

Low energies effects
At low energies, it is very well known that the electronic stopping tents to damp in a non proportional velocity dependence [17], due to at very low relative velocities (compared to v 0 ) colliding nuclei will not penetrate the electron clouds of each other strongly, this defines the Coulomb repulsion effects for electronic stopping power.
The general formula for electronic stopping according kinetic theory [18] S e = (Ξ)N mv where R is the distance of closest approach, that is computed by solving E = 1 2 V (R) for a given inter-atomic potential, σ t r is the transport cross section, N e the electron density, v F the Fermi velocity of the electron gas and Ξ a geometrical factor previously mentioned by Tilinin [9]. This factor is only relevant for high atomic number elements, like Ge, where for Si the effect can be neglected and approximate Ξ ≈ 1. Here we consider three different models for the analysis, Tilinin [9] (Transport cross section and kinetic theory), Kishinevsky [19] (based on Firsov model with inter-atomic interaction included ) and Arista [20] (based on the dielectric function formalism, based in Lindhard ideas). Furthermore we consider four inter-atomic potentials;Thomas-Fermi, Moliere, AVG and Ziegler [21,22]. Both the scale factor ξ e and the inter-atomic potential used affects the electronic stopping and the binding energy, specially at low energies.
For the binding energy model we include the Frenkel energy, i.e. the energy to create a free ion in the crystal lattice, mainly important at low energies, and we include inner excitation of the electron clouds, that are mainly important at high energies. In general high binding energies, tent to reduce the quenching factor.
Usually Density Functional Theory is used to model electronic stopping and binding energies. The kinetic energy for electrons usually is taken to be the average free electron energy of a Fermi gas (3/5)E F where this assumption may work at high energies. But for low energies interaction, Tilinin [9] makes the observation that only electrons near the Fermi energy can be excited, due to Pauli exclusion principle, since inner electron occupy the energy levels. So if we take in to account this observation, we can change in the model that the average electron kinetic energy to be just E F . By doing this it can be shown that the atomic scale change by a factor of 5/3. This implies for the electronic stopping the appearance of an scaling factor of ξ = 2.15. Before this Lindhard observed the need to ad a corrective factor to electronic stopping in the range among one or two. This argument may give a reasonable explanation about the physical origin of this scale factor in the context of density functional theory. In this study we are going to consider this atomic scaling factor ξ, as a free parameter among the range 1 to 2.15, furthermore this also affects the binding energy model and nuclear recoils interactions.

5 Ionization Efficiency for Noble Gases
For noble gases we can apply our model for the quenching factor for LXe and LAr, using recent measurements of the ionization yield Q ER y and the scintillation efficiency L y [23]. Reconstruction is done by exploiting the full anticorrelation between the S1 (scintillation photons n γ ) and S2 (ionized electrons n e ), where it is usually to assume that each excited or ionized atom leads to one scintillation photon or electron, hence for the total quanta defined by the number of ions N i and excitons N e x produced, N i + N ex = n γ + n e independent of recombination.
The fraction of ionizations due to recombination is predicted by the Thomas-Imel box model [24].This model has been shown to work well for spatially small tracks. The charge yield is proportional to Eq.(1) by the electron recoil energy E er = f n E R and the energy W i to create an electron-hole pair in the liquid, where γ is a free parameter of the model typically of the order of 10 −2 . In analogy for the light yield, where N ex /N i will be considered as a constant, the typical values are of the order of ≈ 1/2. We compute the charge and light yields for LAr and LXe, using the constant binding energy model and S e = kϵ 1/2 , as a first attempt to explain recent low energies measurements.

Results and Applications
We show in Fig.(1), the preliminary result for Si, compared with the constant binding energy model and Lindhard's. Now we can describe five orders of magnitude of data. We fit the interatomic scale parameter to electronic stopping power (S e ) data for Si-Si ions [25][26][27], that scales the electronic and nuclear stopping and has important effects at low energies. In this study we obtained Ξ = 1.46 to be a best fit for S e data. Even though the model is independent of Lindhard's, we see a very good match with our model at high energies. Excess signals produced by flat background can be expected to appear at low energies due to new effects, binding energy, considered by Eq.(3) than Lindhard's model fail to reproduce. For flat nuclear recoil spectrum or background signals, e.g. thermal neutrons, we can explain an EXCESS signal considering that typically Monte Carlo (MC) simulations used Lindhard's theory to reconstruct from recoil to visible energy in experiments. If in substitution we use our model for quenching factor, we reconstruct a signal that have an EXCESS at low energy compared from this MC simulations. The physical origin of this excess relies in the interatomic interactions due to binding energy (energy stored in a defect of the crystal), electronic and nuclear stopping.
Defect production in crystals is a recent are of study [28], that has potential to be applied for low energy detection. With our model and the Kinchin and Pease theory [29] , is very straight forward to compute the number of Frenkel pairs in a crystal in Si. We note that the expected number of events for our model is great than the prediction by using Lindhard's model.
For germanium we can use Eq.(4) to compute the electronic stopping power and introduce it to Eq.(3) (with straggling) for f n computation. But as mentioned before Tilinin approach fails for large atomic numbers like Ge. In this case it is necessary to introduced and model

Conclusions
We present a general model based on integral equations for ionization in pure crystals and noble liquids. The model predicts the turnover of f n at low energies, already observed in Xe for E R < 1 keV. We incorporate corrections due to electronic straggling and atomic scaling in the Int. Diff. Eq. For silicon Coulomb effects allow us to fit the data up to 3 MeV and have a threshold near Frenkel-pair creation energy. For germanium our model shows potential to explain recent measurements [31]. We show charge and light yields for LXe and LAr consistent with actual data. Much work can be done from here, e.g directional quenching factor, straggling forν, etc.