Light and charge yields from nuclear recoils in noble liquids: From sub-keV energies to 1 MeV

Introduction

Over the last decade, dual-phase time projection chambers (TPCs) based on liquid argon (LAr) and liquid xenon (LXe) have emerged as the leading technology for the direct detection of dark matter (DM) [1, 2, 3] and the observation of coherent elastic neutrino-nucleus scattering (CEνNS) [4, 5, 6]. The exceptional scalability, high target density, and superior electron-recoil discrimination of these media make them ideal for probing physics beyond the Standard Model (BSM) via nuclear recoil channels.

The observation of CEνNS [7] has opened a new window for precision phenomenology. As a low-threshold energy process, CEνNS serves as a powerful tool for measuring the electroweak mixing angle at low momentum transfer and providing stringent constraints on non-standard interactions (NSI), neutrino magnetic moments, and sterile neutrinos [8, 9, 10]. However, the success of these precision programs and the possibility of distinguishing a DM signal from the so-called “neutrino fog” depend critically on accurate calibration of the detector response [11, 12]. Specifically, the conversion of nuclear recoil energy into detectable quanta (photons and electrons) in the sub-keV regime remains the dominant source of systematic uncertainty.

It is essential to have a comprehensive theoretical model that describes the ionization and scintillation yields of low-energy nuclear recoils. This is a difficult task. Standard Lindhard theory [13], while foundational, fails to describe the quenching behavior at energies below 10 keV because it neglects the underlying atomic structure and non-linear electronic stopping effects [14, 15, 16]. For its part, semi-empirical frameworks, such as the Noble Element Simulation Technique (NEST) [17], rely on extensive fitting to experimental data. Although effective, these models often lack a fundamental predictive basis and may require different parameterizations for different experimental setups.

Nuclear collisions are governed by the competition between electronic and nuclear stopping processes. Atom-atom interactions involve more degrees of freedom than those between electrons and atoms, and the energy deposited by them is not fully converted into electron-equivalent energy. In the following, we outline the main stages of the physical processes for the production of quanta in noble liquids [15, 18]:

1. The process begins when a non-ionizing particle hits an atom of a noble element and transfers an initial recoil energy. If this energy is greater than the minimum required to create quanta, further collisions will occur until the ions or valence electrons scatter with an energy below this minimum value. In cascade collisions, the momentum transferred to the nuclei and the momentum transferred to the electrons are proportional to the nuclear and electronic stopping powers, respectively.

2. Since the recoil velocity of atoms is much lower than the typical Fermi velocity of valence electrons, it is appropriate to describe nuclear collisions as an adiabatic, elastic, or inelastic process, depending on whether the total kinetic energy is conserved or not during the process.

3. In a binary atomic collision, the system has many degrees of freedom, and hence, the average energy to create an ion during the cascade is generally not constant. We classify the dielectric process as follows:

• Passive channel: During the binary collision, the projectile electron remains in its ground state. Since it is strongly bound, it behaves as a “heavy” particle that transfers momentum to the nucleus. A quasi-molecule is formed during the process because of the partial shielding of the nuclear charges by the inner electron cloud. This reduces the target

ionization cross section compared to a bare nucleus [19]. Subsequently, due to the exchange of momentum, the target electrons can remain in their orbitals, become excited, or ionize. Therefore, depending on the magnitude of the momentum interchanged, either elastic or inelastic processes can occur in this channel. The passive inelastic channel is also known as the “single-inelastic channel” because only one of the atoms is excited or ionized.

• Active channel: This is essentially determined by direct Coulomb interactions between target and projectile valence electrons, which are emitted with an average energy proportional to the valence band width. For the energies involved, the closest approach distance is not short enough for the resulting ions to form a quasi-molecule. Consequently, only inelastic processes can happen in this channel, which is also known as the “double-inelastic channel”.

The passive and active channels modify the ionization cross section in opposite ways. One or the other can dominate, depending on the relative magnitudes of the projectile and target nuclear charges, and the projectile’s energy.

(iv) The probability of elastic scattering between atoms is proportional to the ratio of electronic to nuclear stopping power. At low energies, elastic collisions in the passive channel dominate.

(v) The elastic electron momentum exchange determines the electronic straggling in atom collisions [20]. This process is proportional to the ratio of the electronic straggling to the electronic stopping power multiplied by energy. As energy increases, this ratio decreases.

(vi) The total electronic energy, also called the equivalent electronic energy (eVee), corresponds to the sum of the energies of all nuclear recoils in the cascade with electron emission. The remaining energy, associated with nuclear motion, is converted into heat. Electrons, ions, and excitons form a plasma along a path determined by the first recoiling atom.

(vii) In the absence of an external electric field, the electron escape probability is negligible. Consequently, all electrons generated by a nuclear recoil recombine with ions, contributing to the scintillation light.

(viii) By contrast, in a dual-phase TPC, the applied electric field displaces electrons into the gas phase (S2), while photons are captured in the liquid phase by an array of photomultiplier tubes (S1). A region is created where partial electron recombination occurs, the size of which decreases as the intensity of the external field increases, with the consequent reduction of the recombination probability.

(ix) The probability that an ion captures two electrons during recombination is negligible. Thus, we can consider that there is only one electron per ion and one photon per exciton.

(x) In a cascade with high exciton density, there exists a considerable probability of exciton-exciton collisions, producing a single survival excitation.

(xi) If the noble liquid is pure enough (less than 10 ppb) to reduce light-quenching agents such as water and oxygen, then electrons and photons travel through the TPC as the primary detection signal.

Using the second-order integro-differential equation [21], derived from the Linhard integral equation [13] for atomic motion, we incorporate the effects of electronic and nuclear stopping powers, binding energy, and electronic straggling. From the solution to this equation, incorporating realistic inter-atomic potentials, we implement a unified description of the light and charge yields. A distinctive feature of our model is that it characterizes the medium’s response across a broad energy range from sub-keV energies to 1 MeV, in terms of only one adjustable parameter, which is associated with the recombination process. Our method accounts for the energy and electric field dependence of the recombination process through a modification of the Thomas-Imel box model which is valid for different noble liquids. The results show a significant improvement over traditional empirical models, providing the theoretical precision required by the next generation of CEvNS and DM experiments.

The paper is organized as follows. In § 2, we describe the production of quanta in noble liquids and introduce the total quanta model for noble liquid detectors, covering energy conservation, the passive elastic and inelastic channels derived from the second-order integro-differential equation, the electronic scaling parameter $\xi_{\mathrm{e}}$, and the active-inelastic channel. § 3 presents the active Bates-Griffing electron scattering model, describing the energy dependence of the exciton-to-ion ratio and incorporating spin exchange effects in binary atomic collisions. Section IV describes the recombination process within an extended Thomas-Imel box model, where the box size is derived from electrostatic equilibrium, along with a polarizability model for the electric field dependence of $N_{\mathrm{ex}}/N_{\mathrm{i}}$ and an analytical calculation of Birks’ constant. In § 5, we show the unified recombination predictions for light and charge yields across different noble liquids, comparing our results with experimental data and NEST predictions. Finally, Sec. VI summarizes our main conclusions and highlights the implications for future dark matter and CEvNS experiments. Detailed mathematical derivations and atomic calculations are provided in the Appendices.

Total quanta model for noble liquid detectors

Energy distribution

In dual-phase noble liquid detectors, the nuclear recoil energy $E_R$ that a particle (e.g., DM or neutrino) transfers to the target nucleus is partitioned into three channels: excitation of atoms $N_{\rm ex}^{R}$ (excitons), formation of electron-ion pairs $N_i^{R}$ (ions), and atomic motion $E_\nu$. Consequently, the conservation of energy is expressed by the Platzman equation [22]

where $E_{\rm ex}^{R}$ and $E_i^{R}$ are the average nuclear recoil energies used to create excited states and electron-ion pairs, respectively.

For electron recoils, the transferred energy $E_e$ is divided between the production of excitons and ions; therefore, the total quanta is $N_q^e = N_{\rm ex}^e + N_i^e$. The average energy to create an exciton or ion can be defined as $W \equiv E_e/N_q^e$ and from this the mean electron ionization energy is given by $W_i^e = W(1+N_{\rm ex}^e/N_i^e)$, where $N_{\rm ex}^e/N_i^e$ is called the (electronic) exciton-to-ion ratio. For LAr, $W = 19.5 \pm1.0$ eV and $(W_i^e)_{\rm LAr} = 23.6 \pm0.3$ eV and $(W_i^e)_{\rm LXe} = 15.6 \pm0.3$ [23]. For LXe, we use the average value from Dahl [24] ($W = 13.7 \pm0.2$ eV) and the more recent value from Zurich II [25] ($W = 11.5 \pm0.25$ eV), $W = 12.84 \pm0.2$ eV, this value is closer to that calculated with the electron-free gas model, see Appendix B. This yields electronic exciton-to-ion ratios of $(N_{\rm ex}^e/N_i^e)_{\rm LAr} = 0.21$, and $(N_{\rm ex}^e/N_i^e)_{\rm LXe} = 0.06$.

In this work, we study the light and charge yields of the photons and electrons originating from nuclear recoils. In general, the exciton to ion ratio is not equal to the value obtained from electron interactions. The binary collision includes a pair of nuclei with screened charges; consequently, more excitation levels are available and $N_{\rm ex}^{R} \simeq N_i^{R}$. The nuclear quanta efficiency $f_n(E_R)$ (also known as the Lindhard factor) is defined as

where $E_{ee} = N_q^R W$ is the electron energy equivalent. Here, $N_q^R = N_i^R + N_{\rm ex}^R$ is the total nuclear quanta and $E_R = E + U$, with $E$ the kinetic energy of the moving ion and $U$ the energy absorbed by atomic electrons. From this equation, one obtains the relation $f_n(E_R)E_R = W_i^R N_i^R$, where

is the mean ionization energy, with $\beta= N_{\rm ex}^R/N_i^R$.

The basic parameters needed to determine $f_n(E_R)$ in noble liquids are the electronic stopping power, $S_e$, the nuclear stopping power, $S_n$, the electronic straggling, $\mathcal{W}$, the atomic binding energy, the spin exchange energy, $E_{\rm ex}$, and the atomic number, $Z$. Two contributions to $f_n(E_R)$ can be identified: passive and active. The passive component $f_n^{\rm PI}(E_R)$ is determined by the competition between $S_e$, $\mathcal{W}$, and $S_n$, during electronic momentum transfer in the quasi-molecule. In contrast, the active component $f_n^{\rm AI}(E_R)$ is governed by the direct electron-electron interaction, which induces ionization without the formation of the quasi-molecule, provided the projectile energy is greater than $W_i^R$. In the rest of the article we will restrict ourselves to the nuclear recoil process and omit the index $R$, except for the recoil energy $E_R$ of the nucleus that initiates the cascade.

Passive elastic and inelastic channels

We determine $f_n^{\rm PI}(E_R)$ from the solution to the second-order integro-differential equation derived from the Lindhard integral equation [21]:

where $\varepsilon= QE$, $\bar{\nu}(\varepsilon) = QE_\nu$, and $u(\varepsilon) = QU$, with $Q = 16.26[Z_1Z_2(Z_1^{0.23}+Z_2^{0.23})\,\mathrm{keV}]^{-1}$ [26]. The integration variable is $t = \varepsilon^2\sin^2(\theta/2)$, with $\theta$ the scattering angle in the center of mass frame. It is pertinent to point out that the original Lindhard formula for $\bar{\nu}(\varepsilon)$ [27] is an approximate solution to the first-order integro-differential equation, with $u=0$ and an electronic stopping power valid for energies above 10 keV (bare ions). This precludes its applicability to lower energies.

The solution to the second order equation in Eq. (4) that we derive here incorporates the effects of atomic binding energy, electronic straggling in nuclear collisions, and electronic stopping. For the nuclear stopping power, we use LSS theory [27], with ion-ion differential cross-section parametrized according to Ziegler $f(t^{1/2})$ [16]. The effect of the binding energy is modeled as the sum of an $\varepsilon$-dependent term, $u_{TF}(\varepsilon)=Q U_{TF}(\varepsilon)$, plus a constant term, $u_c$, fixed by an energy threshold proportional to the average energy needed to create an exciton in a noble liquid (see Appendix A). $U_{TF}(\varepsilon)$ is given by the average energy to excite electrons in the atomic collision, calculated from the Thomas-Fermi (TF) formalism, and, as shown in Figure 1, it depends on the atomic number of the noble element.

Vector (PDF)

We adopt the model of Tilinin [28], in which, for impact parameters small compared with $a$, the electronic stopping power is governed by the mean free path of the electron, $\lambda_T$, during the formation of the quasi-molecule, and the energy transferred elastically by the electron per collision, $\Delta E_T=(m(v+v_F)^2-mv_F^2)/2\simeq mvv_F$, where $m$ is the electron mass, $v$ is the ion velocity and $v_F$ is the Fermi velocity of the electron, with $v\ll v_F$ in the quasi-molecular system. $S_e$ is obtained by integrating the volume of the quasi-molecule and multiplying the result by the atomic density, $N_b$, namely,

with $\lambda_T^{-1}=n_{eT}\sigma_{tr}$, where $n_{eT}$ is the electronic density in the quasi-molecule and $\sigma_{tr}$ is the momentum transfer cross section (MTCS) of the electron. The electronic density determines the Coulomb repulsion effects [15, 21] and can be approximated by the sum of the contributions of the projectile, $n_{eT}^{P}$, and the target, $n_{eT}^{T}$. For atoms with the same mass and atomic number, we have $n_{eT}^{P}=2n_{eT}^{T}$. This model better describes the low energy region than the one in [27] and works well for silicon [21, 29]. However, it uses a semiclassical momentum transfer cross section for electron-atom scattering, and its applicability is limited to energies higher than a threshold proportional to the atomic number (1.4 keV for Ne and 70 keV for LXe). To examine the low-energy regime, we implemented a solution for the MTCS using a partial-wave analysis with the ZBL potential for the atomic numbers of Ne, Ar, Kr, and Xe. In Figure 2, we compared our computed values of $\sigma_{tr}$ with available experimental data for these noble-gas elements [30, 31, 32, 33, 34]. The comparison is presented as a function of the reduced energy, $\varepsilon_e=32.52 E_e\,\mathrm{keV}\,Z^{-4/3}$. In addition, we report the mean trend obtained by averaging across all species, which highlights the near-universal behavior of the MTCS at electron energies exceeding 10 eV. We have used the standard Fermi scale, $a=0.885a_0 Z^{-1/3}$, where $a_0$ is the Bohr radius, and add the data for other elements reported in [35]. In this way, the validity of our model extends from $\varepsilon_e=0.1$ to $\varepsilon_e=10^{-3}$. For $E_R=2\sqrt{2}ZE_e$, this interval corresponds approximately to nuclear recoil energies from 10 keV to 0.1 keV, for LXe. As shown in the Figure, this energy region corresponds to the peak of

the cross section for both LAr and LXe. Since existing data are for noble gases, extending the applicability for liquid phase to lower energies would require considering the liquid structure and polarization potentials [36, 37]. For electronic straggling $W$, in a first-order approximation, we assume that the electronic and nuclear straggling contributions are equal. Consequently, we adopt the expression provided by [20], in which $W$ is evaluated on the basis of the nuclear straggling. Given that $W$ appears only through the second-order derivative in Eq. (4), rescaling the electronic straggling by a factor $M_0$ induces only a small modification of $f_n(E_R)$, with the effect being more pronounced at low recoil energies.

Electronic scaling parameter

As discussed in [21], the electronic stopping power must be corrected with the fermion factor $\xi_e$, which accounts for the energy dependence of the electronic excitations and depends on the media. Then, the scaling length becomes $a=0.885a_0 Z^{-1/3}\xi_e^{2/3}$, and the electronic stopping power can be written as $S_e(\varepsilon,\xi_e)=\xi_e S_e(\varepsilon)$. Lindhard et al [38] derived the expression $S_e(\varepsilon,\xi_e)=\mathcal{K}\xi_e^{1/2}$ for bare ions, where $\mathcal{K}=0.133(Z/A)^{1/2}\xi_e$, with $\xi_e\simeq Z^{1/6}$.

For high-energy collisions, all electrons up to the Fermi level can be effectively excited (Lindhard assumption) and $\xi_e\to1$, whereas for low energies, only electrons close to the Fermi level can be excited and $\xi_e\to(5/3)^{3/2}\simeq

2.15. Instead of the fit for $\xi_e$ given in [21], here we define the scaling factor from Thomas-Fermi (TF) theory, using the plane-wave approximation for electrons:

where $E(k)=\hbar^2 k^2/2m$, $k_F$ is the average Fermi wave-number of electrons, $\Delta k$ is the maximum electron dispersion in the collision, which depends on the projectile energy, and $E(k_F)=U_0/\xi_e^{2/3}$, with $U_0=24.04Z^{1.23}$ eV for a two atom system and $U_0=20.92Z^{4/3}$ eV for a single atom.

Associated with $\Delta k$ there is an average energy available for excitation below the Fermi level, which we denote as the average Fermi energy of skin penetration $\Delta E_F$ [39]. Another interpretation of $\Delta k$ is as follows. Suppose that electrons scatter elastically, then $\Delta k=k_F\sin(\theta_c/2)$, where $\theta_c$ is the scattering angle in the CM frame, as shown in Appendix C. For $\theta_c=\pi$, we have $\Delta k=k_F$, and only radial collisions are allowed [40]. For angles above $\theta_c$ we have damping due to energy transfer from electrons to inner shells, then only electrons up to the energy $\Delta E_F$ have to be taken into account.

For noble elements, the $\xi_e$ factor has to decrease as the atomic number increases, because the TF model has been shown to be exact in the limit $Z\to\infty$ [41]. For the cascade process, the value of the momentum transferred is related to the average kinetic energy of the electrons in the atom, and consequently the electronic stopping power depends on this energy. Thus, the upper limit for $\xi_e$ can be estimated (in the gas phase) using the expected TF average kinetic energy for electrons in an atom,

For $\xi_e=1$ the average kinetic energy is overestimated, because Eq. (7) ignores the exchange effect and sudden variations in electron density near the nuclei and in the valence orbital region. The maximum value for $\xi_e$ can be calculated using the formula (7) for the hydrogen atom and $\langle T_F\rangle=\langle T_e\rangle=13.6057$ eV, which gives us $\xi_e^{\max}=1.91$.

For multielectron atoms we can determine $\langle T_F\rangle$ using atomic orbital calculations based on local density approximations [42, 43], see Table I. It can be shown that the values for $\Delta E_F$ follow the approximate relation $\Delta E_F\simeq0.29Z^{2/3}$ eV.

Element	$\langle T_F\rangle_{RLDA}$ eV	$\xi_e$	$\Delta E_F$ eV
He	37.7	1.66	1.2
Ne	348.4	1.47	26.2
Ar	797.1	1.38	85.3
Kr	2114.4	1.27	347.0
Xe	3793.9	1.19	855.6

The electron-electron interaction can be significant for low $Z$ elements, resulting in an ideal fermion system with lower $\Delta E_F$. According to the TF approach and Eq. (6) electrons with kinetic orbital energy up to $\Delta E_F$ behave as free particles. Therefore, for light elements, only the outer electrons fulfill this assumption. In the limit $Z\to\infty$, we have $\Delta E_F\to E_F$ as expected.

The estimation of $\xi_e$ for noble liquid elements can be done by adding to $\Delta E_F$ the gap energy to create a pair of electrons, $E_g$. Then, we can compute the $\xi_e$ factor for the liquid phase by means of Eq. (6), see Table II. This addition only affects those elements with an atomic number lower than $Z=36$ and, in particular, for LAr the correction is only of 3%. In the case of He, the MIMAC experiment [44] indicates a lower value than the Lindhard model ($\xi_e=1.12$), in better agreement with the present model.

Element	$(\Delta E_F+E_g)$ eV	$\xi_e$ (liquid phase)
He	25.8	1.03
Ne	47.7	1.33
Ar	99.6	1.34
Kr	358.7	1.26
Xe	864.9	1.19

Active-Inelastic channel

For elements with a high atomic number such as LXe, the Tilinin model, used to account for the contribution of the passive channel, yields a low value of $f_n(E_R)$. The discrepancies between the predictions of the model and the observed data are illustrated in Figure 3. This indicates that other electron emission processes, particularly inelastic electron emission, should be considered.

Vector (PDF)

When an electron associated with the projectile atom interacts with the inner-shell electrons of the target atom, an electronic interaction potential energy $u_A^e$ arises. This quantity is taken to correspond to the atomic electronic binding energy, $u_{TF}(\varepsilon)$, as described within the Thomas-Fermi framework. By contrast, direct interactions between the valence electrons of the projectile and those of the target predominantly drive inelastic electron emission during the cascade process. These electrons are emitted with an energy $u_A=\bar{\varepsilon}_D$, where $\bar{\varepsilon}_D$ is the average electron kinetic energy associated with the valence-orbital energy $\varepsilon_D^0$. From [48] we can estimate $\varepsilon_D^0$ using two models; Shockley $\varepsilon_D^0=3E_{\mathrm{ex}}/5$ and Alkhazov’s $\varepsilon_D^0 = E_{ex}/2$, where $E_{ex}$ is the exciton energy, see Table III.

Element	$E_{ex1}$ (eV)	$(\varepsilon_D^0)_{\mathrm{sub}}$ (eV)	$(\varepsilon_D^0)_{\mathrm{Shockley}}$ (eV)
He	19.80	9.90	11.88
Ne	17.50	8.75	10.50
Ar	12.06	6.03	7.24
Kr	10.17	5.09	6.10
Xe	8.37	4.19	5.02

Consequently, the total energy satisfies the following

where $u_P^\nu$ is the nuclear recoil energy associated with the first process, and we define the total atomic binding energy as $u = u_P^e + u_P^\nu= u_{\mathrm{TF}}(\varepsilon) + u_c$. The resulting total-quanta efficiency is

where we identify $\bar{\varepsilon}_D/\varepsilon_R = f_n^{\mathrm{AI}}(\varepsilon_R)$ as the active channel contribution. Inelastic electron emission is energetically allowed only when the projectile energy exceeds the threshold $W_i$, which is equivalently expressed as requiring an interaction distance smaller than $r_M$. At $r_E = r_M$, the emission probability $p(r_M)$ attains its maximum. For $E > W_i$ (i.e., $r_E < r_M$), electron emission may still occur, but with a reduced probability. The expected number of emitted electrons is given by $p(r_E)E/W_i$, and the corresponding mean electron energy is $\bar{\varepsilon}_D = p(r_E)\varepsilon_D^0 E/W_i$. Incorporating this inelastic energy-loss channel into the quenching factor yields:

Given that the process is governed by Coulomb interactions between electrons, the scaling relation $E/Z^2 \simeq1/r_E$ follows. Consequently, $p(r_E)$ is expected to increase proportionally to $Z^2$. The IDES formulation for $f_n(E_R)$ is determined by the underlying interatomic interaction potential, the effective scaling length $\xi_e$, and the atomic number and atomic mass of the target medium. In the present work, the effective atom–atom interaction is modeled using the Ziegler potential [16].

Active Bates-Griffing Electron Scattering

Energy-dependent $\beta$ model

In the active channel of collisions between pairs of identical atoms, one electron from each atom is released from its orbital through its energy-dependent electrostatic interactions with the electrons and nucleus of the atom it collides with. This process is modeled in terms of the Bates-Griffing (BG) scattering adapted to the present case for energies of the order of keV. We use a semiclassical approach, assuming that only valence electrons participate in the process. These are located at the radii of their valence orbitals. Since the electron clouds are not expected to change significantly during the collision, the average positions of the electrons can be considered nucleus-centered. When the two atoms are separated by a distance $x'$, ignoring spin effects, a valence electron from one atom interacts through the Coulomb potential of the valence electrons and the shielded nucleus of the other atom, see Figure 4.

The atom-atom system is described by a multi-electron wave function $\psi(r_e,R)$, where $R$ is the nucleus separation and $r_e$ denotes the set of vector positions of electrons in both atoms. In the Born-Oppenheimer approximation $\psi(r_e,R) = \psi_e(r_e)\psi_N(R)$ [50]. We define $x = R/a_Z$, in terms of the Ziegler units $a_Z = 0.885 q \xi_e^2 a_0/(2Z^{0.23})$, with $Z_1 = Z_2 = Z$. Here, $q$ is a correction factor for atom resizing, which, ignoring relativistic corrections, can be determined by the Hartree-Fock model [51]. The respective values for LNe, LAr, LKr, and LXe are 0.667, 0.92, 0.95, and 1.08. For LXe, a relativistic contraction correction is needed, giving $q = 1.05$ [52, 53], while for LAr, using data from Ref. [20], we get $q = 0.942$.

We use the semi-empirical potential energy [16],

where $\chi_Z(x)$ is the Ziegler screening function:

Consider an incoming atom with kinetic energy $E$, then the minimum approach distance, $x_m$, is obtained from the equality $E=e\phi_Z(x_m,\xi_e)$.

At the end of the cascade process, for a given incoming energy $E$, the atoms involved in the collisions may have an effective number of electrons $N<Z$ determined by the Bohr stripping criteria [54, 55], see Appendix § 2.4. Assuming that $N$ remains constant during the recombination process, the potential describing a neutral atom ($Z=N$) must be replaced by a potential characterized by an effective charge $Z-N$. This modification can be implemented using the Amaldi model. Within this framework, the electron–ion potential energy can be written as

with $x'=r/b$ and $b=q\,0.885\,\xi_e^{2/3}a_0N^{-1/3}$. For details on the validity of this potential, see Appendix § 12.

The orbital radius of the valence electron with $l\ne0$ is determined by condition $e\phi_Z^e(x'_0,\xi_e,N)/2=-U_e$, where $U_e$ is the average electron binding energy of the valence electrons: 21.6 eV, 15.82 eV, 14.1 eV, and 12.57 eV for Ne, Ar, Kr, and Xe, respectively [56, 57].

The total energy of the system is given by the sum of negative projectile binding $E_U=-U_e/2$, negative projectile-target binding $E_y=e\phi_Z^e(y',\xi_e,N)/2$, with $y'=\sqrt{x'^2+x_0'^2}$, and positive Coulomb electron-electron repulsion ($T_C=0$) $E_C=e^2/(bx')$. For an electron to be emitted from the projectile, the total energy must be zero, which gives the transcendental equation,

We can determine $x'$ by solving Eq. (13) for a given $N$. Subsequently, the associated ion-ion potential energy, denoted as $W_i$, required to ionize an electron from one of the ions can be calculated.

with $x=x'b/a_Z$. Using the definition (14),

we can obtain the exciton to ion ratio. Since $N$ is expected to decrease with $E_R$ the screening in the potential $\phi_Z^e$ will decrease, and therefore the right-hand side of (14) will increase with energy; see Figure 5, Figure 6, Figure 7 and Figure 8. Consequently, we infer that $N_{\rm ex}/N_i$ has to be an increasing function of energy, with minimum value at the limit $E_R\to0$ where $N\to Z$. $N_{\rm ex}/N_i$ is expected to also be sensitive to the scale parameter in the allowed interval of values $1\le\xi_e\le2.15$. For values near $\xi_e=2.15$ the electrons below the Fermi energy cannot be ionized and only excited by the exchange mechanism, hence it is expected that $N_{\rm ex}/N_i$ be greater than for the opposite case. For electron recoils in LAr and LXe, $(N_{\rm ex}/N_i)_{\min}=0.23$ and $(N_{\rm ex}/N_i)_{\min}=0.06$ [58, 59] may be considered the lowest limit value for the respective nuclear recoil ratios.

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Spin effects

Because of permutation symmetry between the electrons, for a multi-electronic wave function with spin there is an extra interaction energy contribution to the total Hamiltonian of nonclassical origin (exchange interaction) [61]. Consequently, a more precise description of the BG process can be given by considering the spin wave function in the binary process.

For two approaching atoms, at separation distance $x\gg1$, the total potential energy of an electron in each atom is given by the sum of the two individual potentials: $V_1(x'_1)+V_2(x'_2)$, where $x'_1$ and $x'_2$ are the corresponding electron coordinates. Taking into account the spin, the orbital wave function of the electrons becomes $\psi(x'_1,x'_2)=(\psi_1(x'_1)\psi_2(x'_2)\pm\psi_1(x'_2)\psi_2(x'_1))/\sqrt{2}$, where we can define the average Coulomb potential $\langle\psi| V_C | \psi\rangle$ for $x \gg x_i$ ($i = 1, 2$). At closer distance we have to include the perturbation potential $V_{nn} + V' + V'_{ee}$, where $V_{nn}$ is the nuclei-nuclei potential, the electron potential in atom 1 with the nuclei in atom 2, $V'$, and the electron-electron interaction between both atoms, $V'_{ee}$. Since $|V' + V'_{ee}|/V_{nn} \cong x'_i/Zx \ll1$, the potential $V_{nn}$ makes the dominant contribution.

Consequently, the perturbation potential is positive, and the orbital energy level shifts to higher values. To gain more stability, electron orbitals can lower their energies through exchange interaction, which, according to the generalized Hund rule [62], flip the spins to point in the same direction. For noble element quasi-molecules,

we can have three electron spins down for the $p_{1/2}$ (shifted) orbitals and three spins up for $p_{3/2}$. In this way, the total spin is conserved ($S = 0$) and the multiplicity is maximized.

For the singlet, the spatial wave function has to be symmetric and the Coulomb potential energy will have two terms, the direct term $e^2/(bx)$, considered previously, and the exchange term $V^{ex}_{ee}$; $V_C = e^2/(bx') + V^{ex}_{ee}$. In this situation, since the overlap of the wave function is not too strong, for distances ($x \cong2x'_0$) the active BG process is favored, where Eq. (15) is satisfied ($x \cong2x'_0$). Thus, $e^2/(b2x'_0) \gg V^{ex}_{ee}$ and Eq. (13) remains valid as a good (classical) approximation.

For the triplet case, the spatial wave function is antisymmetric, leading to a potential energy reduction $V_C=e^2/(bx')-V_{ee}^{\mathrm{ex}}$. For distances $x<2x'_0$ the overlap of the wave function becomes significant and disfavors the active BG process. Therefore, this case may lead to a passive BG process. In principle, the active channel can be computed from the interaction between two atoms and Eq. (4). This calculation requires the emission cross section, which depends on the atomic orbitals of both atoms, as well as on the energy and spin degrees of freedom. A rigorous treatment of this process would require a spinor wave function obtained from the Dirac equation [52]. Such a relativistic analysis lies beyond the scope of the present work; therefore, we adopt a semiclassical approach.

To examine the probability of electron spin exchange, we account for the Coulomb energy of the quasi-molecule and the exchange energy. We assume that the Coulomb interaction is dominant. We define the average electron exchange interaction distance $x_{\lambda}=(a_Z\sigma_tN_e)^{-1}$ and, assuming a Poisson process, determine the probability of an odd number of flip-flops during the collision, $N_{ff}=x/x_{\lambda}\geq1$:

where $P_s$ is the singlet probability ($\cong1/4$), $\sigma_t$ is the total cross section, $\sum\sigma_i+\sigma_{el}$, $\sigma_i$ is the ionization cross section for each sub-level $i$ and $\sigma_{el}$ is the elastic cross section, $N_e$ the electron density of the quasi-molecule that contributes to the dielectronic process (by exchange interaction) and $x$ the minimum separation distance of the atoms; for details, see Appendix G.

From Eq. (16) it can be argued that for low $Z$ elements $N_e$ is expected to be low, then $x\ll x_{\lambda}$, so the emission probability is suppressed $p\propto x/x_{\lambda}\ll1$. For Ne and Ar, the BG effect is negligible compared with that observed for heavier elements such as Kr and Xe.

Recombination process: Thomas-Imel box model

The recombination process is relevant only for liquid-phase detectors because the probability of its occurrence in solids is negligible. The partition between ionization and scintillation produced in a noble liquid by the initial recoil of an atom is governed by the recombination probability $p_r$ between the ions and excitons that were generated [15, 63]. The Thomas-Imel box model [63] assumes that the recombination mechanism is described by a classical process of ions and electrons recombining in a box of size $a_B$ in the presence of a uniform electric field $F_0$, given in units of V/cm.

It is generally assumed that each ionization produces a single electron and that each excited atom and each recombination emit a scintillation photon [23]. Then, for $N_i$ ion-electron pairs generated by the cascade process in the box, the number $N_e$ of electrons escaping recombination predicted by Thomas-Imel is

where

where $\mu_e$ is the electron mobility, $v_e=\mu_eF_0$ is the electron drift velocity and $\alpha$ is a constant recombination parameter with dimensions of $\mathrm{m^3\,s^{-1}}$. Scintillation is equal to exciton quanta plus the fraction of recombined ions, that is, $N_{\gamma}=N_{ex}+p_rN_i$. The remaining $(1-p_r)N_i$ electrons escape recombination. Therefore, $N_{ex}+N_i=N_e+N_{\gamma}$.

Biexcitonic process

Excitons can self-trap at a frequency similar to that of the electron-hole can recombination process. Before trapping, free excitons can diffuse and undergo biexcitonic quenching, $i.e.$, $\mathrm{Ar}^{*}+\mathrm{Ar}^{*}\rightarrow\mathrm{Ar}+\mathrm{Ar}+e$. This mechanism is proportional to the energy deposited and produces light quenching in noble liquids [64]. Consequently, the total quanta is also quenched [65]:

where $f_l(E_R)$ is the biexcitonic quenching factor, which is related to the Penning effect and can only be present in nuclear recoil ionization. Some authors [15, 66] have parameterized this process in terms of Birks’ saturation law,

where $k_B$ is called Birks’ constant. For this model to work, the timescale of exciton self-trapping must be shorter than the recombination time at this recoil energy [67].

Let us consider the total quanta efficiency measured in a dual phase detector $f_n^*(E_R)=(N_{\gamma}+N_e)W/E_R$; from Eqs. (2) and (19), it can be expressed as:

Some observations are in order. First, $f_n(E_R)$ is an upper limit for the measured efficiency, and second, $f_n^*(E_R)$ depends slightly on the electric field through $N_e$. Because $N_e/E_R$ is much smaller than the light yield for $E_R$ greater than 50 keV and $f_l(E_R)\cong1$ for recoil energies lower than this value, we can make $f_n^*(E_R)\cong f_n(E_R)f_l(E_R)$.

Using $N_i=(E_R/W)f_n(E_R)/(1+\beta)$, and taking into account Eqs. (19), (17), and (2), we obtain the following expressions for the charge and light yields:

In recent works [49, 58, 60], the recombination parameters $\beta$, $\xi$, and $\alpha$ are determined from direct measurements. In some studies, the model depends on up to five parameters for each noble liquid [68] and, particularly for LAr, the consistency with charge and light yields is not entirely satisfactory. Next, we present a comprehensive model that explains the recombination parameters and the total quanta of all noble elements as functions of a single experimental parameter. In a subsequent publication, we will provide the explicit formulas and software to implement the model.

For recoil energies below $\simeq1$ MeV ($v \ll Z^{2/3}v_0$ and $Z \ge10$) the recombination is not expected to drastically alter electronic orbitals, so that when the cascade process ends and exciton-ion recombination begins, $\Delta k$ should not change significantly. It is important to note that the previous equations strongly depend on $f_n(E_R)$, the exciton-ion ratio $\beta$ and implicitly on the length of the electronic scaling $\xi_e$. Considering the interpretation of Dahl [24], we will assume that the box size is determined by the electrostatic length scale, which implies a non-hyperbolic dependence of $\xi$ on $F_0$.

Box size description in terms of the electrostatic electron-ion interaction

In this part of the study we present a model for the box size parameter $a_B$, valid for any noble liquid detector, which explains the electric field behavior of the parameter $\xi$ in Eq. (18). The model depends on the microscopic properties of the noble liquid, the scale length $\xi_e$, and the electrostatic potential.

During recombination, electrons and ions coexist until thermalization begins to contribute significantly. In the Thomas-Imel model this is characterized by the box size parameter $a_B$, see Eq. (18), with the box assumed to be isolated and the ions and electrons uniformly distributed. Typical values of $a_B$ are of the order of micrometers [24, 69].

In many studies [70, 71], for the fit function of $\xi$, $a_B$ is taken as an average quantity and the dependence on the electric field is a power law: $\xi\propto1/F_0^\gamma$, where $\gamma\ne1$ is a nonuniversal parameter that depends on the medium. Instead, we will follow Dahl’s ideas [24] and consider that $a_B$ is determined by the electrostatic equilibrium between electron-ion pairs and the applied electric field, namely, $\phi_{\mathrm{ext}}(z,\xi_e)=-bzF_0$, where $z$ is the average electron-ion separation. To estimate the lowest electric field for which this criterion can be applied, we use the Coulomb potential for the electron-ion interaction and assume that

the electron’s time inside the box $t=2a_B/v$ is less than the thermalization time ($\simeq1$ ns), which gives us $F_0>10$ V/cm. These assumptions are supported by several observations [24, 49, 58, 65, 68].

Because the box size is expected to be much larger than the atomic scale, screening effects due to the induced charge and atom polarization have to be included. Since we are not in a relativistic regime, we can use the linear Thomas-Fermi screening theory [72], where a semi-classical approximation for a slowly varying potential is assumed. Additionally, for low energy ion tracks ($<1$ MeV), only the valence electrons need to be considered inside the box, thus the electron-ion potential becomes

where $k^2=6\pi e^2n_e/\Delta E_F$ is the Thomas-Fermi wave vector, with $n_e$ denoting the electron density in the box and, according to Eq. (6), $\Delta E_F$ is determined by the value of the scaling length $\xi_e$. In this model, the number $N$ of effective electrons for each atom inside the box is given by the average ion energy, $E_i=2Z(\Delta E_F/2)$.

We define the electron density $n_e^0$ in thermal equilibrium by $N_v n$, where $n$ is the number density of atoms in the medium at a certain equilibrium temperature and $N_v$ is the number of valence electrons. The electron density is expected to vary from $n_e$ to $n_e^0$. To establish an approximate relation between $n_e$ and $n_e^0$, we make the following specific assumptions:

The electron density is proportional to the mean free path $\lambda_e$ of electron-ion interactions, $n_e\propto1/\lambda_e^3$. Here, $\lambda_e=1/n_i\sigma$, where $n_i$ is the density of ions within the box and $\sigma$ is the elastic cross section for the electron-ion dispersion, which, for $\Delta E_F\gtrsim10$ eV, is proportional to $Z/\Delta E_F$ [28].

The density of ions is proportional to the mean separation $x_i$ of the ions inside the box, determined by the ion-ion potential $\phi_Z(x_i,\xi_e)$ in the weakly coupled dense plasma approximation [73]. This, for the interval from $x_{\min}=2x_0$ to the thermal equilibrium distance $x_{\mathrm{th}}=2r_W/a$, where $r_W$ is the Van der Waals radius.

For simplicity, we approximate the ion-ion potential by a power law, $\phi_Z(r,\xi_e)\propto Z^2/r^\Omega$, where $\Omega$ is a noble liquid dependent parameter, with $\Omega>1$ [74], see Figure 9. Then, we can establish an explicit proportionality relationship for the average ion energy for the entire range of $x_i$: $\Delta E_F\simeq\phi_Z(x_i,\xi_e)/Z\propto Z/x_i^\Omega$. This relation is valid only if $q\ll1$, with $q=(x_{\mathrm{th}}-x_{\min})/(x_{\mathrm{th}}+x_{\min})$. Then $n_i\propto(1/x_i)^3\propto(\Delta E_F)^{3/\Omega}Z^{-3/\Omega}$.

Vector (PDF)

From the previous results, for the mean electron-ion scattering length, we get

where $\lambda_0$ is a proportionality constant with units of length. Then, the electron density becomes

with $(\Delta E_F)_0$ being the average electron energy for thermal electrons, which can be estimated from the gap energy, $E_g$. Thermal equilibrium is reached when $\Delta E_F \to(\Delta E_F)_0$ for $\Omega> 3$ [75]. On the other hand, condition $q^2 \ll1$ limits the applicability of Eq. (26) for high energies ($\Omega< 3$), in which the electron density, $n_e$, increases with energy.

To determine the size of the box, $a_B$, we equal the screened electron-ion potential to the external field potential

and solve for $z$ as a function of the electric field $F_0$. Then, the size of the box is determined by doing $a_B = z(F_0)$. For an electric field of about 10 V/cm we obtain $a \simeq10$ nm for LAr and $a \simeq80$ nm for LXe. As expected, in both cases, the box size is smaller than the Onsager radius, which is about 100 nm. From Eq. (18) we see that now $\xi\propto1/(a_B^2F_0)$ and its dependence becomes more complicated than the hyperbolic one. In general, $a_B(F_0)$ is expected to be a monotonic decreasing function of the electric field, and its derivative with respect to $F_0$ is sensitive to $\xi_e$. The minimal value of $F_0$ for Eq. (27) to be valid can be estimated by equating the box size to the Onsager radius, with is related to the average thermal energy $E_T$, with the formula $F_3 = (E_T/E_g)(W_i/E_g - 1)$, where $F_3$ is the excitation Fano factor [48]. For LAr, $F_3 = 0.004$ and then $E_T = 9$ meV, which corresponds to an Onsager radius of 10 nm, and a minimal electric field of 10 V/cm. In the case of LXe, $F_3 = 0.0006$, $E_T = 1$ meV, and the Onsager radius is 67 nm, which gives a minimal electric field 20 V/cm.

For simplicity, to have some insight on the solution of Eq. (27), we can perform a power-law fit to the box size, $a_B = a'_0(F'_0/F_0)^\delta$ around a specific electric field value, where $a'_0$ and $F'_0$ are fit constants. Typical electric fields in liquid TPCs are approximately 100–1000 V/cm. As discussed above, the slope of the box size depends on the scale parameter $\xi_e$. For noble liquids considered here we obtain different values for the slope $\delta$ at different electric fields, see Table IV. The slope $\delta$ depends on $\xi_e$, the atomic number (for noble elements), and the electric field. From Eq. (27), we find $\delta_{\max} = 1/2$, reached only for a purely Coulomb interaction potential. For elements with large atomic numbers, the valence electrons are located relatively far from the nucleus, and the effective nuclear charge estimated using Slater’s rules is approximately equal to the net charge experienced by the valence shell after shielding by the other valence electrons. Consequently, the screening parameter $\delta$ is expected to increase and approach its maximal value as the atomic number increases.

For $\delta$ values corresponding to $\xi_e$ in Table II we get the fit function,

This result demonstrates that, to a very good approximation, the recombination parameter exhibits an inverse proportionality to the electric field and obeys a power law dependence, $\xi/N_i \propto1/F_0^{(1-2\delta)}$.

Element	$\delta$ ($10^2$ V/cm)	$\delta$ ($10^3$ V/cm)
Neon	0.161	0.211
Argon	0.203	0.270
Krypton	0.288	0.371
Xenon	0.385	0.449

Polarizability model for $N_{ex}/N_i$

When an electric field $F_0$ is applied to a noble liquid, the ions and atoms that recombine are expected to interact with the field. This interaction produces a small displacement $d$ of the outer electron orbitals. For an electric field applied along the $z$ direction, assuming a spherical shape of the atom in the absence of the field, to first order, the deformation due to the displacement is $x' = x'_0 + x'_d \cos\theta$, where $x'_d = d/b$ and $\theta$ is the angle with respect to the field direction.

The potential of the atom-atom system, $V(r_e,x,\xi_e)$, is given by the sum of the following contributions: nuclei-nuclei $V_{NN}(x)$, electron-electron $V_{ee}(r_e)$, and electron-nuclei $V_{eN}(r_e,x,\xi_e)$. We consider this interaction in the context of a semiclassical theory, where the quantum interference effects for distances lower than $x$ are ignored.

The average potential is given by

For an atom with $Z$ electrons, the electron-atom potential can be defined in a similar manner as in Eq. (29). For a test charge at a distance $x^{\prime}$ from the nuclei,

where $e\bar{Z}(x^{\prime})=4\pi b^3\int_{x^{\prime}}^{\infty}\rho(x^{\prime\prime})x^{\prime\prime2}dx^{\prime\prime}$ is the effective charge and $\rho(x^{\prime})=e\left|\psi_Z(x^{\prime})\right|^2$. Here, $\psi_Z(x^{\prime})$ and $V_Z(x^{\prime})$ are the wave function and the effective potential, respectively, resulting after the integration of the electronic coordinates.

The correction $\delta\phi_Z^e(x^{\prime},F_0)$, due to the external field potential $\phi_{\mathrm{ext}}(x^{\prime},\xi_e)$, modifies the total potential for electron-atom and atom-atom systems,

We determine both shifts in terms of the Stark effect [76] for a deformed charge distribution using nondegenerate first-order perturbation theory [77],

To evaluate the double integral on the right-hand side of this equation, we make the following approximations: (i) Since $x_d^{\prime}\ll x^{\prime}$, $\rho(x^{\prime}+x_d^{\prime}\cos\theta)\cong\rho(x^{\prime})+x_d^{\prime}\cos\theta\,d\rho^{\prime}/dx^{\prime}$, the radial integration is separated and, only the term dependent on $\cos^2\theta$ gives a nonzero result to the integration over $\theta$. (ii) Only interaction with valence electrons is considered and consequently $x^{\prime\prime}\gg1$; therefore, we can use the Kohn-Sham asymptotic electron density, $\rho(x^{\prime})\cong\rho_0e^{-\gamma x^{\prime}}$, with $\gamma=2(2U_{e1}/E_h)^{1/2}q\,\xi_e^{2/3}0.885Z^{-1/3}$ [78], where $E_h$ is the hydrogen binding energy ($\cong13.6$ eV), $U_{e1}$ is the second ionization energy of the atom, and $\rho_0$ is a constant associated with the electron density in the region close to the valence electrons. (iii) We assume that $x_d^{\prime}=\alpha_pF_0/eb$ and $\alpha_p$ is the polarizability constant. Using these approximations we get,

Since we use first order perturbation theory and assume $x^{\prime}\gg x_d^{\prime}$, the validity of Eq. (34) is limited to electric

fields not very intense. Since the shift is negative, the exciton to ion ratio $N_{\mathrm{ex}}/N_i$ is expected to be a decreasing function of $F_0$.

Using the same procedure as in § 3 A, we compute the ratio $N_{\mathrm{ex}}/N_i$ by means of Eqs. (13) and (15) but with the shifted potentials defined in Eq. (31). The results of this computation for LXe and LAr are shown in Figure 10 and Figure 11. From these figures, we see that the upper value for the electric field is of the order $\sim10^4$ V/cm, which is sufficiently large for the current noble liquid experiments [69].

Vector (PDF)

Birks constant

The coefficient $k_B$ in Birks’ law (Eq. (20)) is related to the probability that two excitons form a bonded system at high energies. Assuming a simple orbital bonding, the symmetric wave function (singlet) is required and $k_B\propto k'_B/4$, where $k'_B$ is another constant. From the Birks law, the constant $k'_B$ has to be proportional to $S_e(\epsilon_0,\xi_e)^{-1}$, where $\epsilon_0$ is the energy of the ions such that the effective ion charge, $Z_{\mathrm{eff}}$, allows the production of an exciton pair. The effective charge for an exciton must be in the range $0<Z_{\mathrm{eff}}<1$, thus, to estimate $\epsilon_0$ we use the Ziegler effective charge formula for ions, $Z_{\mathrm{eff}}=Z(1-e^{-v_i/v_0Z^{2/3}})$, where $v_i$ is ion velocity, such that, $Z_{\mathrm{eff}}=1/2$, see Table V. From the energies listed in Table V we can make an analytic estimate using Lindhard’s electronic stopping power, including the effect of the correction factor $q$ for the inter-atomic potential, with

Ion	$Z$	$A$	$\varepsilon_0/Q$ (keV)
Ne	10	20.180	27.92
Ar	18	39.95	35.96
Kr	36	83.798	48.1
Xe	54	131.293	56.98

$S_{\mathrm{e}}(\varepsilon_0,\xi_{\mathrm{e}}) = q^{3/2}0.133\xi_{\mathrm{e}}(Z\varepsilon_0/A)^{1/2}$, then

where we can use the approximate trend $\varepsilon_0 \simeq121.52/Z^{1.91}$ ($\varepsilon_0/Q \sim Z^{1/3}$) to give a closed formula for $\xi_{\mathrm{e}} k_B$ for all noble elements, see Table VI.

Element	$q^{3/2}\xi_{\mathrm{e}} k_B$, (35)
Ne	2.18
Ar	3.97
Kr	7.99
Xe	12.01

The results obtained for LAr and LXe are in good agreement with the experimental measurements, [15, 79]. This model predicts that for higher $Z$ the biexcitonic quenching is more active.

Recombination for noble liquids

We now show that the ratio $\xi/N_i$ can be determined for all noble liquids in terms of a single fit constant. It is convenient to write $\alpha= \alpha_0 r_0^2 v_{\max}$, where $\alpha_0$ is a dimensionless constant, $v_{\max} \simeq v_0$ is the maximum drift velocity of the electrons, and $r_0$ is the minimum box length. The electron mobility is $\mu_{\mathrm{e}} \propto\lambda_{\mathrm{e}}/v_0^{\mathrm{e}}$, where $\lambda_{\mathrm{e}}$ is given in Eq. (25) and $v_{\mathrm{e}}^0$ is the electron velocity within the box, $v_{\mathrm{e}}^0 = (2\Delta E_F/m)^{1/2}$. From the empirical formula for the atomic radius $r_0 \simeq0.18Z^{0.44}\mathring{\mathrm{A}}$, Eq. (18), and Eq. (28), we find

where $C_0 = (\alpha_0/\lambda_0)(271v_0m^{1/2}(10\,\mathrm{eV})^{1-3/\Omega}/e)$, where $\alpha_0/\lambda_0$ is an unknown constant.

To find $C_0$ we have to determine $\delta$. From the fit with the power law of Eq. (28), in the range $100\ \mathrm{V/cm} < F_0 < 2000\ \mathrm{V/cm}$ for the electric field, we obtain $1-2\delta= 0.106$ for LXe and $1-2\delta= 0.53$ for LAr, see Figure 12. Using the exact solution of Eq. (27) we obtain the curves in Figure 13 for $1-2\delta$ as a function of the electric field. We also show the mean values, $<1-2\delta> = 0.082$ for LXe and $<1-2\delta> = 0.51$ for LAr, in the range $100\ \mathrm{V/cm} < F_0 < 10000\ \mathrm{V/cm}$. We observe that the exponent behavior is nearly constant with the electric field. Except for the cases where the model is not reliable for low electric fields, $F_0 < 50\ \mathrm{V/cm}$ or near the breakdown electric field, $F_0 > 10^4\ \mathrm{V/cm}$.

Vector (PDF)

Vector (PDF)

The dimensionless parameter $4\xi/N_i$ as a function of the applied electric field is one of the main comparison tests to validate the physical consistency of our model. The dependence of this parameter on the electric field is not universal for the noble liquid elements. To determine the constant $C_0$ we fit the experimental data for LXe [49] using a power-law function of the form $P_0/(F_0/(1\ \mathrm{V/cm}))^{1-2\delta}$, which yields $(4\xi/N_i)_{\mathrm{LXe}} = (0.05961 \pm0.034)/(F_0/(1\ \mathrm{V/cm}))^{0.082}$. This fit results in a value of $C_0 = 1.88 \times10^3\ \mathrm{m^3s^{-1}(eV)^{-3/\Omega+1/2}(V/cm)^{1-2\delta}}$. Using this value in Eq. (36) for LAr, the model predicts $(4\xi/N_i)_{\mathrm{LAr}} = 1.36/(F_0/(1\ \mathrm{V/cm}))^{0.51}$, which is consistent with the direct experimental fit $(4\xi/N_i)_{\mathrm{LAr}}^{\mathrm{fit}} \simeq(1.41 \pm0.18)/(F_0/(1\ \mathrm{V/cm}))^{0.51}$ as illustrated in Figure 14.

Vector (PDF)

The predictions of the model for LKr and LNe are $(4\xi/N_i)_{\mathrm{LKr}} = (0.219\ \mathrm{m^3s^{-1}})/(F_0/(1\ \mathrm{V/cm}))^{0.26}$ and $(4\xi/N_i)_{\mathrm{LNe}} = (12.12\ \mathrm{m^3s^{-1}})/(F_0/(1\ \mathrm{V/cm}))^{0.52}$. The general trend for any noble liquid with $Z > 10$ is predicted to be $(4\xi/N_i)_Z = 2650[Z^{2.9-1.5/Z^{1/2}}]^{-1}/(F_0/(1\ \mathrm{V/cm}))^{1-2\delta}$. For the sake of completeness we present the results for total quanta using Eq. (9) for LXe, LAr, GNe and LNe. The BG effect for neon is negligible and is not shown, see fig. 15.

The constant $C_0$ allows us to estimate the recombination probability constant, $\alpha_0$, providing deeper physical insight into the underlying modeling. From Eq. (36), we derive the ratio $\alpha_0/\tau_c = 2.2 \times10^7\ \mathrm{s}^{-1}$, where $\tau_c$ is the mean collision time within the recombination box. This quantity is related to electron mobility by $\mu_e = e\tau_c/m$ and can be estimated in terms of the plasma frequency for LXe ($\omega_e^2 = 4\pi e^2n/m\epsilon_0$), giving us $\tau_c \simeq2\pi/\omega_e = 3.4 \times10^{-16}\ \mathrm{s}$, which yields $\alpha_0 \simeq7.5 \times10^{-9}$. By applying the Jaffe equations to the stationary regime for

both ions and electrons [63], the recombination constant can be approximated as $\alpha_0 \simeq(\tau_r N_- r_0^2 v_F)^{-1}$, where $\tau_r$ is the recombination time for nuclear recoils and $N_-$ denotes the electron density, estimated as $N_- \simeq N v_n$. Proceeding in such a way, we obtain $\tau_r \simeq0.033\ \mathrm{ns}$, which is in good agreement with the results in the literature $\tau_r \ll1\ \mathrm{ns}$ [24, 49]. Finally, in Figure 17-Figure 21, we present the light and charge yields for LAr and LXe predicted by our model and compare them with the data and NEST predictions. We further compare our model predictions with the mean light-yield efficiency, $L_{\rm eff}$, in liquid argon (LAr) at 36.1 keV under varying applied electric-field strengths, see Figure 16.

Vector (PDF)

Vector (PDF)Vector (PDF)Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)Vector (PDF)Vector (PDF)Vector (PDF)Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Vector (PDF)

Conclusions

In this work, we have developed a first-principles theoretical formalism to describe total quanta production in noble liquids over an energy range spanning from the sub-keV regime to 1 MeV. The commonly used Lindhard formula is only valid for $\varepsilon> 1$ and assumes an electronic stopping power proportional to the velocity (only the parameter $k$). In contrast, our total quanta model is elaborated from the second-order integro-differential equation, which accounts for a general electronic stopping power, electronic straggling, and nuclear stopping power. The electronic stopping is characterized using Tilinin’s non-linear analysis of the momentum transfer cross-section, incorporating recent experimental data for liquid argon and xenon, as well as low-energy Coulomb repulsion and high-energy Bohr stripping criteria.

We demonstrated that elastic momentum transfer effects alone are insufficient to characterize ionization measurements in noble elements with $Z > 20$. A significant contribution of this work is the inclusion of the Bates-Griffing dielectronic process to account for inelastic electron emission during the atomic cascade. The explicit mechanism is formulated in terms of the electrostatic potential, the spin configuration of the diatomic system, the exchange energy, and the kinetic energy of the projectile atom.

Furthermore, we introduced the electronic scaling length, $\xi_e$, as a new fundamental physical parameter for nuclear recoil interactions. This parameter incorporates the Pauli exclusion principle within the Thomas-Fermi model, providing a more precise description of electronic orbitals during the nuclear cascade process. Within the recombination framework, $\xi_e$ gives a better characterization of the exciton to ion ratio, implemented as a function of the energy and the electric field, rather than as a constant value. The recombination process is described by an improved version of the Thomas-Imel box model, where the box size is derived from the electrostatic equilibrium between the external electric field and the field associated with the electrons and ions inside the box. This provides a physically grounded description of the field dependence, without the need for an empirical power-law fit.

Finally, the model shows overall consistency for the noble liquids most commonly used in experiments: neon, argon, and xenon. All of them can be described by a single fitting parameter, while other relevant parameters are derived from fundamental atomic properties. The results presented here show that our approach outperforms empirical methods, such as NEST, particularly in the low-energy regime, predicting substantially higher light and charge yields. The model can provide an even more accurate description if we add an additional fitting parameter to describe an energy dependence of $\xi_e$. These findings highlight the qualities of the model as a robust tool for the next generation of ultra-low threshold experiments, including the search for dark matter and the detection of coherent elastic neutrino-nucleus scattering.

Acknowledgments

This research was supported in part by SECIHTI grant CF-2023-I-1169, a SECIHTI postdoctoral grant, and DGAPA-UNAM grant PAPIIT-IN104026. The authors also acknowledge SNII support.

Energy equilibrium for the cascade process

Consider the collision between two atoms, with the target and projectile having $Z_1$ and $Z_2$ electrons and masses $M_1$ and $M_2$, respectively. The total kinetic energy of the target atom for each atomic collision in the CM frame is, $T_N = E^\prime_R M_2/(M_1+M_2)$, where $E^\prime_R$ is the nuclear recoil energy for each atomic collision in the lab frame. When the quasi-molecule is formed we assume that all electrons remain in a bound state, such that all the energy released to the system is absorbed by them,

with $T_{ei}$ being the energy of the electron $i$. We can define an average electron energy, $\overline{T}_e$, such that $T_e = (Z_1+Z_2)\overline{T}_e$,

In the CM frame, $\overline{T}_e \simeq\overline{T}^{\mathrm{lab}}_e \equiv E_e$, since the electron mass is much lower than the atom mass; hence,

The energy $E_e$, can be seen as the average energy of the electron-atom interaction during atomic collision. For identical atoms, $E_e = E_R/4Z_1$, which coincides with the formula obtained using the Ziegler scale and can be compared with the Tilinin formula with the Lindhard scale, $E_e = E_R/2\sqrt{2Z}$ [28].

For comparison, let us test the group IV of the periodic table. In particular the minimum energy to excite an electron from the valence band to the Fermi level is $E_g/2$. For an electron energy below this value, no ionization can occur, then $E_e \approx E_g/2$. The minimal recoil energy, $E_{Rm}$, can be calculated from, $E_{Rm} = 2Z_1E_g$. For energy gaps of 5.47, 1.17, 0.67 eV, we obtain $E_{Rm} = 65.7$, 32.8, 42.9 eV, for C, Si, and Ge, respectively. In the case of Si, the obtained value is close to those obtained in other studies [80]. According to the solution of the integral equation for the atomic motion (Eq. (4)), the recoil threshold is approximately equal to $2U_F$, where $U_F$ is the Frenkel pair energy. From the previous results it can be directly obtained $U_F = E_{Rm}/2$, which compares well with previous studies [81, 82, 83].

In noble liquids, the minimum energy to produce a metastable state is given by the spin-exchange, where the energy range to form this state is about 0.1 − 1.0 eV [84, 85]. Using eq. (A3) this gives, $4\ \mathrm{eV} < 2U_c < 216\ \mathrm{eV}$, since the energy range under study is above 500 eV, as a first approximation we can take the average value, $U_c = 50\ \mathrm{eV}$ for all noble liquids.

Free-electron model for the mean energy required to create an electron-hole pair

To estimate $W$, we employ the free electron gas model in the limit $T \to0$. In this framework, the change in internal energy, $\Delta U$, of the valence band electrons due to an energy transfer ($\Delta E \propto k_B\Delta T$) from an ionizing particle is given by the following relation:

where $E \ll E_F$ is the energy transferred to the Fermi gas and $E_F = E_V + E_g/2$, with $E_V$ denoting the valence-band energy. The energy $\Delta U$ corresponds to the average energy required to excite a valence electron, $W_V$, multiplied by the total number of excitations. Thus, assuming that the derivative $\partial U/\partial E$ remains constant, we obtain $E/W_V = \pi^2E/2E_F$, where it follows that $W_V = 2E_F/\pi^2$. To account for the additional energy $E_g$ required to promote an electron from the valence band to the conduction band, the total average energy $W$ is given by:

The width of the valence band can be estimated from the free-electron model as $E_V \simeq\hbar^2/2m(3\pi^2n_e)^{2/3}$, where $n_e$ is the number density of electrons ($N_v n$), see Table VII and Table VIII.

Material	$\rho$ (g/cm$^3$)	$n_e$ (m$^{-3}$)	$E_V$ (eV)	$E_g$ (eV)
LNe	1.207	$2.88 \times10^{29}$	15.91	21.50
LAr	1.395	$1.68 \times10^{29}$	11.11	14.30
LKr	2.413	$1.38 \times10^{29}$	9.77	11.55
LXe	2.940	$1.08 \times10^{29}$	8.26	9.30
Si	2.330	$1.99 \times10^{29}$	12.47	1.17
Ge	5.323	$1.76 \times10^{29}$	11.48	0.66
Diamond	3.510	$7.04 \times10^{29}$	28.86	5.48

Material	$W_{\mathrm{exp}}$ eV	$W$ eV	Eq. (B2)
$CDiamond$	13.0-13.3	11.88
Si	3.66-3.75	3.70
Ge	2.8-3.0	2.99
LNe	25-27	26.9
LAr	19.5	18.0
LKr	14.2-15.9	14.7
LXe	11.5-13.7	11.92

Free electron plasma dispersion

Consider the average Fermi energy for a free electron gas, $\bar{E}_{F} = \frac{3}{5}(\hbar^{2}k_{F}^{2}/2m)$, where $k_{F}$ is the maximum wave vector that an electron can have. The variation of this energy,

where $\mathbf{k}_{F}\cdot\Delta\mathbf{k}_{F}=k_{F}\Delta k\cos\psi$, with $\psi$ the collision angle between the electrons in the laboratory frame. On the other hand, to estimate the maximum value of $\Delta E_{F}^{M}$, we can use Lindhard dielectric function [38], to compare with (C1),

where $\hbar\omega_{F}=mv_{F}^{2}/2$. By equating (C1) and (C2) [86], we obtain the critical collision angle $\psi_{c}$,

where $\theta_{c}$ is the critical angle in the CM frame.

For angles $\theta<\theta_{c}$, then $\Delta E_{F}<\Delta E_{F}^{M}$ and there is no significant energy transfer from free electrons to bonding electrons. In the opposite case, the free-electron gas must be damped, due to direct transfer of all its energy to a single inner orbital electron. This can be interpreted as meaning that for electron energies lower than $\Delta E_{F}$ we can approximate the electrons as free. In the context of TF theory, inner electrons with energies larger than $\Delta E_{F}$ do not contribute significantly to momentum transfer in atomic collisions.

Effective number of electrons in atomic binary collisions

From appendix A we see that the total electron energy given to the system of identical ions during a collision is $E_{e}=E_{R}/4Z$. Using the adiabatic criteria [55], $v\leq\sqrt{3}v_{F}$, we have $E_{F}\approx E_{e}/3$. Then calculate the effective electron distance for the electron to be stripped, $E_{F}=\hbar^{2}(3\pi^{3}n_{F}^{e})^{2/3}/2m$, where $n_{F}^{e}(x')=$

Amaldi model for ionized atoms

According to the Thomas-Fermi model for the atom, the screen function $\chi_{Z}(x')$, $x'=r/a$ with $a=a_{0}0.885/Z^{1/3}$, obeys the equation

with the boundary conditions $\chi_{Z}(0)=1$ and $\chi_{Z}(x')=0$. For an ionized atom with $N$ electrons, the boundary condition at $x'_{0}$ becomes $\chi_{Z}(x'_{0})=(Z-N)/Z$. Amaldi’s approach [87] approximates the ionized boundary condition, imposing the ansatz,

with $(Z-N)/Z\ll1$ and $\kappa$ a correction factor for screening length to be determined. Let us notice that this ansatz satisfies the boundary conditions for an ionized atom mentioned above.

Requiring that the expression Eq. (E2) satisfy (E1), we find

if we set $\kappa^{-3/2}(N/Z)^{1/2}=1$, we deduce Amaldi’s screening length, $a_{A}=a_{0}0.885/N^{1/3}$ screening length.

Polarizability for free electron gas during recombination

Here, we apply the Lindhard dielectric function theory [38] for a gas of free electrons in the presence of a static electromagnetic field. In his semiclassical model the dielectric function is,

where $\omega_{e}^{2}=Ze^{2}\bar{n}_{e}/m\varepsilon_{0}$ is the electron plasma frequency (in standard units), with $\bar{n}_{e}$ given by the Thomas-Fermi model:

is the TF electron density. Having the distance $x'$ we can compute the effective number of electrons in the binary collision by

where $\chi_{Z}$ is the screen function associated with $\phi_{Z}(x',\xi_{e})$. where $U_{TF}=20.9\,Z^{4/3}/\xi_e^{2/3}$ eV. In Eq. (48) $v_F$ is the Fermi velocity associated with the average energy of the valence electrons during recombination, $k_Fv_F=(2/\hbar)(\Delta E_F/N_v)$.

The polarizability is given by

where $\bar{N}_i$ is the average ion density during recombination. Thus,

from the previous result the polarizability is given by,

The average ion density can be estimated as $gN$, where $N$ is the atomic number density, and $g=(r_W/\bar{r})^3$, with $\bar{r}=(r_W+r_0)/2$. For LXe $g=2.37$, $\alpha_p/(4\pi\varepsilon_0)=1.30\times10^4\,\mathring{\mathrm{A}}^3$, and for LAr $g=6.07$, $\alpha_p/4\pi\varepsilon_0=7.77\times10^3\,\mathring{\mathrm{A}}^3$. To test Eq. (50) we estimate the polarizability of LXe and LAr in thermal equilibrium, in which case, instead of $U_{TF}$ we have to use the first ionization energy of the atom: 12.8 eV for LXe and 15.9 eV for LAr. Thus, we obtain $\alpha_p\simeq4.5\,\mathring{\mathrm{A}}^3$ for LXe and $\alpha_p\simeq2.4\,\mathring{\mathrm{A}}^3$ for LAr, which are similar to the measured values of $4.0\,\mathring{\mathrm{A}}^3$ and $1.7\,\mathring{\mathrm{A}}^3$, respectively [88].

Mean free path for electron ionization in atomic collisions

In § 2 A, we discussed the dielectronic emission process as a result of the collision between two atoms whose electrons interact through the electron-electron potential and the screened nucleus potential. At the threshold energy $W_i$, the electrons are emitted with an average energy $\bar{\epsilon}_D$. For energies larger than $W_i$ electrons feel the attractive potential of both approaching atoms and, as a consequence, there is a probability of no emission. The probability of electron pair creation is related to the mean free path of the sum of the elastic and inelastic cross sections times the electron density in the binary collision.

For electrons with energies above $27.2\,Z^{2/3}$ eV, the elastic cross section is given by [28]

where $E_0=13.6$ eV. The inelastic channel is given by the sum of all ionization (or excitation) processes [89]:

where $t_i=E/U_i$, $S_i=4\pi a_0^2n_i(E_0/U_i)^2$, $u_i\simeq1/2$, and $Q$ is the dipole constant defined in terms of the continuum dipole oscillator strength. Here, $U_i$ and $n_i$ are the corresponding binding energy and orbital occupation number.

The electron density, $N_e$, relevant for the BG process in atomic collisions, is related to the exchange energy. Using the TF model, we estimate as

where $\bar{E}_{ex}$ is the average exchange potential energy in eV, calculated in terms of the sum of all orbital contributions. According to the electron free gas model [90], we obtain

Here, $N_i$ is the electron density in orbital $i$ and

with $\eta_i=k_i/\Delta k_F=\sqrt{U_i/2\Delta E_F}$, where by the virial theorem $(T_i=U_i/2)$. The exchange energy contributes only if $\ell\ne0$ and $\eta_i\le1$, for $k_i\ge\Delta k_F$ electrons are above the Fermi level and, if $F(\eta_i)=0$, consistent with the interpretation of $\Delta E_F$. To compute $N_i$ we use the Thomas-Fermi model with the AVG potential [20] evaluated at the corresponding radii, $x'_i=a_0n^2/bZ^{i}_{eff}$, where $n$ is the principal quantum number of the orbital and $Z^{i}_{eff}$ effective charge of orbital $i$, see Table IX and Table X.

$i$	Shell	$Z^{i}_{eff}$	$U_i$ eV	$n_i$
1	(5s)	14.22	23.3	2
1	(5p)	12.42	12.75	6
2	(4s)	26.17	213.2	2
2	(4p)	24.96	146.1	6
2	(4d)	21.89	68.5	10
3	(3s)	35.58	1148.7	2
3	(3p)	35.67	971.4	6
3	(3d)	39.95	682.7	10
4	(2s)	39.8	5453	2
4	(2p)	49.83	4946.5	6
5	(1s)	52.92	34561	2

$i$	Shell	$Z_{\mathrm{eff}}^{i}$	$U_i$ eV	$n_i$
3	(3s)	7.76	29.3	2
3	(3p)	6.76	15.8	6
4	(2s)	12.23	326.3	2
4	(2p)	14.01	249.5	6
5	(1s)	17.51	3206	2

For Ar, we need only consider the $3p$ orbital with electron density $8.68\times10^{29}$ e/m$^3$, which gives $\bar{E}_{ex}=1.46$ eV. In the case of Xe, the orbitals $5p$, $4p$, $4d$, $3p$, and $3d$ must be included. Their electron densities, in units of $e/\mathrm{m}^{3}$, are $3.86 \times10^{29}$, $2.94 \times10^{31}$, $2.90 \times10^{31}$, $5.91 \times10^{32}$, and $5.98 \times10^{32}$, respectively, and $\bar{E}_{\mathrm{ex}} = 10.8$ eV.

For the two-atom system, the mean free path is $\lambda^{-1} = N_e\sigma_t(\Delta E_F)$, where $\sigma_t(\Delta E_F) = \sigma_{\mathrm{el}}(\Delta E_F) + \sigma_{\mathrm{in}}(\Delta E_F)$. By assuming a Poisson distribution, the probability of electron emission is $1-P$, with $P=e^{-N_\lambda}$, where $N_\lambda=r/\lambda$ is the number of mean free paths for a separation $r$ of the approaching atoms. In Ziegler units, $N_\lambda= \xi_e^{2/3} 0.885\,q\,a_0 Z^{-0.23}N_e\sigma_t(\Delta E_F)x/2$, where we considered the dependence of the scaling length on the electron density. The probabilities for no emission are $P=e^{-0.068x}$ and $P=e^{-0.5x}$ for LAr and LXe, respectively. Note that the exponential parameters are approximately proportional to $Z^2$, thus for light noble liquids this effect becomes negligible [91].