Migdal effect and photon Bremsstrahlung: improving the sensitivity to light dark matter of liquid argon experiments

The search for dark matter weakly interacting massive particles with noble liquids has probed masses down and below a GeV/c^2. The ultimate limit is represented by the experimental threshold on the energy transfer to the nuclear recoil. Currently, the experimental sensitivity has reached a threshold equivalent to a few ionization electrons. In these conditions, the contribution of a Bremsstrahlung photon or a so-called Migdal electron due to the sudden acceleration of a nucleus after a collision might be sizable. In the present work, we use a Bayesian approach to study how these effects can be exploited in experiments based on liquid argon detectors. In particular, taking inspiration from the DarkSide-50 public spectra, we develop a simulated experiment to show how the Migdal electron and the Bremsstrahlung photon allow to push the experimental sensitivity down to masses of 0.1 GeV/c^2, extending the search region for dark matter particles of previous results. For these masses we estimate the effect of the Earth shielding that, for strongly interacting dark matter, makes any detector blind. Finally, we show how the sensitivity scales for higher exposure.

Direct detection experiments generically are insensitive to nuclear recoils with energy below the keV, corresponding to sub-GeV/c 2 dark matter scattering. This relies on the assumption that the electron cloud around the nucleus follows instantaneously the nucleus itself, keeping the atom neutral. However, the sudden acceleration of a nucleus after a collision may lead to excitation and ionization of atomic electrons. This is an old idea from neutron-nucleus scattering experiments [40][41][42][43][44]. Furthermore, the electron will get accelerated in order to follow the nuclear recoil trajectory, resulting in a finite probability that a photon will be emitted via Bremsstrahlung. Therefore, this new process may lead to energetic photons and ionization electrons produced from the primary interaction. The first process is the Bremsstrahlung photon emission from a nucleus [32], the latter is the Migdal effect [33][34][35], and they both have been already exploited by experimental collaborations [38,[45][46][47].
In the present work we examine both the Migdal effect and photon Bremsstrahlung from the nucleus in experiments exploiting LAr detectors, and we estimate the effect of the Earth atmosphere and crust that makes experiments blind to large cross sections. Previous DM searches of this kind include DarkSide-50 [9] and DEAP-3600 [12], and have mainly focused on masses greater than 10 GeV/c 2 . DarkSide-50 has published a low-mass analysis [28], exploiting the ionization-only signal, sensitive down to masses of 1.8 GeV/c 2 . In this article, we show how this analysis could be extended down to masses of 0.1 GeV/c 2 including signals from the Migdal effect and the photon Bremsstrahlung.
The rest of the paper proceeds as follow. In section 2 we review the Migdal effect and photon Bremsstrahlung process, show their differential rates in the LAr detector, and describe the effect of the Earth attenuation. Section 3 reviews the analysis with LAr. We describe the probabilistic inference and all the details of our simplified treatment of systematic effects. In section 4, we show the expected sensitivity and the projections for higher exposure. Finally we conclude in section 5. In addition, we provide the numerical codes used to evaluate the Migdal and Bremsstrahlung rates [48], and to perform the statistical analysis [49].

Migdal effect and photon Bremsstrahlung
We start this section describing our notation and our assumptions for the elastic DMnucleus scattering rates, and continue presenting the computation of the differential rates for the Migdal effect and the photon Bremsstrahlung process.
The elastic DM-nucleus differential rate with respect to the nuclear recoil energy E R , per unit detector mass, is here N T is the number of target nuclei per unit detector mass, E R is the recoil energy given by an incoming dark matter particle with velocity v > v min = (m N E R )/(2µ 2 N ), m χ and m N are the DM and nucleus mass, respectively, while µ i = m i m χ /(m i + m χ ) is the reduced mass of the nucleus or nucleon-DM system (with i = N, p). The rate depends on our assumptions on the local dark matter density, ρ χ , and the dark matter velocity distribution f ( v). In this work, we use the value 1 ρ χ = 0.3 GeV/c 2 /cm 3 , and the Standard Halo Model [54] with a Maxwell-Boltzmann velocity distribution with dispersion velocity v 0 = 220 km/s and escape velocity cut off of v esc = 544 km/s. The differential elastic DM-nucleus cross section depends on the recoil energy and the DM velocity and is given by 1 New determinations of the local dark matter density give results that fall in the range ∼ (0.3 − 0.4) GeV/c 2 /cm 3 , see [50][51][52][53].
where σ p is the DM-proton cross section (assumed equal to the DM-neutron cross section). The nuclear form factor F (E R ) ∼ 1 for small momentum transfers, while A is the atomic mass, leading to the coherent enhancement of the cross section.

Migdal effect
The rate of ionization due to the Migdal effect for a nuclear recoil energy E R accompanied by a ionization electron with energy E e is given by the standard DM-nucleus differential recoil rate in eq. (1) multiplied by the ionization rate [33] where the ionization rate is given by Here, n and are the initial quantum numbers of the emitted electron, q e = m e 2E R /m N is the electron momentum in the nucleus rest frame immediately after the DM collision, m e is the electron mass, m N is the nucleus mass, and p c qe (n → E e ) is the probability to emit an electron with final energy E e . An approximate estimate of the total energy deposited in the detector is given by E d = E e + E n , where E n is taken to be the binding energy of the (n, ) state. This takes into account the fact that the emitted electron may come from an inner orbital and the remaining excited state will release further energy in the form of photons or additional electrons in order to return to the ground state.
The differential probability rates were computed in Ref. [33] without taking into account the shifts in electronic energy levels because atoms are actually in a liquid (such as in argon or xenon targets) or crystal state (such as in germanium, silicon or sodium iodide detectors). This effect should decrease the ionization energy, and thus, if neglected, should lead to a conservative ionization yield estimate [55]. Figure 1 shows the differential ionization probabilities as a function of the detected energy E d = E e + E n for isolated argon atoms. We use the probabilities computed in Ref. [33].
The accuracy of the differential ionization probabilities relies on the fact that the computed wave functions can reproduce the binding energies for the different levels with an accuracy of ∼ O(20%) and it should provide a correct estimate of the expected signal rate for inner-shell electrons. On the other hand, the prediction for the valence electron shells should be taken as an order of magnitude estimate. 2 The different curves of Fig. 1 show the contributions for different principal quantum number n, where the contributions for different orbital angular momenta in the initial state and all possible final states are summed, and q e = 1 eV. This Figure shows that given the thresholds of 4 or 7 electrons of a hypothetical LAr experiment, the contribution of the valence electrons can maximize the sensitivity to nuclear scattering. This is in contrast with the reported results presented for xenon [38,47] and germanium detectors [45,46], where the detector thresholds are higher than the one needed to see a dominant signal from outer shells. As a consequence, it would be crucial to have a reliable computation of the transition probabilities in the case of LAr.
We stress that the same considerations discussed here are applicable to neutron scattering and should be taken into account by the experimental collaborations when estimating Figure 1: Differential ionization probabilities and related uncertainties as a function of the detected energy E d for isolated argon and different principal quantum number n. We show also the 4 and 7 electron thresholds for DarkSide-50.
In principle, the Migdal effect results can also be compared to DM-electron scattering bounds [34,56,57], in scenarios where the DM couples with equal strength to protons and electrons, as in models with interactions mediated by a dark photon with kinetic mixing [58,59]. Such a connection would need some model dependent assumptions that will impact the generality of our results. In addition, making this connection is beyond the scope of this work and as such requires a dedicated study.

Photon Bremsstrahlung
The displacement of the charges of the nucleus and of the electron after the DM-argon scattering leads to photon emission from the polarised argon atom. Therefore, the elastic nuclear recoil R is the nuclear recoil energy, while ω is the photon energy. This process can provide a detectable signal for dark matter masses that produce elastic nuclear recoils below the detector threshold.
The Bremsstrahlung cross section can be written in terms of the factorised elastic 2 → 2 cross section [32] where E R is the nuclear recoil energy, ω is the photon energy, α is the fine structure constant, and m N is the mass of the target nuclei. The atomic scattering factors f (ω) = f 1 (ω) + if 2 (ω) are tabulated in the NIST Standard Reference Database [60]. Notice that for large energies f 1 → Z f 2 (at 4 keV f 1 ∼ Z). We can then derive the Bremsstrahlung rate as

Rates in argon detectors
We can now show the rates associated with the Migdal effect and photon Bremsstrahlung for dark matter-nucleus scattering in argon detectors. The rates can be obtained by integrating eq. (3) and (6) over those combinations of E R , E e (or ω) and v that satisfy momentum and energy conservation. In the limit of low momentum transfer both the Migdal effect and the photon Bremsstrahlung process share the same kinematics of inelastic dark matter models [61], where the DM mass splitting δm is replaced by the total electronic energy E d or photon energy ω. In particular, we have that where δ correspond to E d or ω for the Migdal or Bremsstrahlung processes, respectively. The maximum nuclear and electronic recoil energy for a given DM mass are This shows that δ max > E R,max for m χ m N due to the suppression factor µ N /m 2 N . Indeed, for v max ∼ 800 km/s ∼ 2.7 · 10 −3 c, m N 40 GeV/c 2 (the approximate argon mass) and a DM mass of 0.5 GeV/c 2 , we find E R,max ∼ 0.09 keV, while δ max ∼ 1.8 keV. As a result, there is a range of DM masses for which it is easier to detect the electronic energy originating from the Migdal or the photon Bremsstrahlung processes rather than nuclear recoils, as a consequence of the fact that more energy can be carried off by light or massless particles for a given momentum transfer.
In particular, eq. (8) shows that experiments exploiting argon detectors (such as DarkSide-50) lose sensitivity for DM masses below 1.8 GeV/c 2 , where E R,max 1 keV. On the other hand, when considering the Migdal effect or the photon Bremsstrahlung process, the DarkSide-50 experiment is sensitive down to m χ ∼ 0.02 GeV/c 2 or m χ ∼ 0.04 GeV/c 2 for electron thresholds of N e − = 4 or N e − = 7, respectively. 3 In Fig. 2 we show the Migdal and Bremsstrahlung rates as a function of the number of detected electrons, induced by a DM particle scattering on argon with a cross section σ SI = 10 −35 cm 2 (σ SI = 10 −33 cm 2 for the Bremsstrahlung process) and mass m χ = 1 GeV/c 2 . In order to compute the rates as a function of the number of detected electrons in a LAr detector, we need to transform eq. (3) and (6) with dE e /dN e − . We find this information using the calibration curve used by the DarkSide-50 collaboration to convert electron recoil spectrum to an average ionization spectrum (Fig. 2 of Ref. [65]). 4 In this context we neglect the fluctuations associated to the detector response which are discussed in Sec. 3.2.2. In addition, in the same plot we show the total background (blue histogram) and the measured spectrum (black histogram) from the DarkSide-50 experiment, taken  Figure 2: Number of events per kg per day for the Migdal effect signal for m χ = 1 GeV/c 2 and σ SI = 10 −35 cm 2 , for n = 1, 2, 3 (red dashed) and n = 1, 2 (solid red). The green curve shows the number of events for the photon Bremsstrahlung signal for a DM mass of 1 GeV/c 2 and σ SI = 10 −33 cm 2 . The histograms show the DarkSide-50 spectrum (black) and background (blue) for an exposure E = 6786 kg d [28].
from Fig. 7 of Ref. [28]. The figure shows also that the Migdal effect dominates over the Bremsstrahlung rate across all energies.
The python code we used for the evaluation of the nuclear recoil, Migdal and Bremsstrahlung rates is publicly available on GitHub [48].

Effects of the Earth attenuation
Direct detection experiments generally lose sensitivity to strongly interacting DM 5 because the same interactions that happen in the detector occur also in the Earth atmosphere and crust. As a consequence, DM particles are slowed down and deflected, reducing the flux of DM particles that arrives at the underground detector. Therefore, there is a critical value of the cross section for which the effect of the Earth attenuation is large enough to make any detector blind to DM interactions. The average energy loss for a DM particle traversing the atmosphere or the Earth crust due to elastic scattering is given by [66,67] where x is the distance traveled, i denotes the different nuclei species encountered and n i (r) is the corresponding number density at position r. Given that E χ = m χ v 2 /2, a change in the DM energy influences the DM velocity distribution at the target. In particular, assuming that all the DM particles move on a straight line from the atmosphere to the detector, the particle flux must be conserved and one can write the DM velocity distribution at the detector as where d is the depth of the underground detector and where µ i is the reduced mass for the nuclei species i, A i its atomic mass and m i its nuclear mass. In this equation we have neglected the effect of the form factor F (E R ) ∼ 1 for light DM.
In order to give an estimate of the Earth's attenuation effect, we compute the DM velocity after the interaction with the Earth's atmosphere and crust using the VERNE code [68,69]. We also assume that the DM travel in a straight line from the surface to the detector, in the Laboratori Nazionali del Gran Sasso (LNGS), Italy. For simplicity, we take the maximal DM velocity in the laboratory reference frame, disregarding the daily and annual modulation. More specifically, we compute the maximal DM velocity at the detector depth and set an upper limit on the detector sensitivity on the cross section solving numerically the following equation for m χ and σ Here the energy δ max is set by the threshold of the experiment. Notice that in this way we overestimate the sensitivity of the experiment to large cross section, but it works as a order of magnitude estimate. We show the results in Section 4.

Sensitivity calculation and LAr simulated experiment
In this work, we adopt a Bayesian approach to infer the upper bound and to estimate the expected experimental sensitivities to the interaction of DM candidates with LAr. Similar approaches for the analysis of DM experimental data have already been deployed [70][71][72][73][74][75][76], although they are not frequent among analyses carried out by the experimental collaborations, as for example [28,39,[45][46][47]65]. Within this approach we can compute, at least in principle, the probability of any specific proposition given some state of information. Theoretical models, parameters of interests, and results of experiments before they are carried out are intended as uncertain propositions connected by the rules of probability. Exploiting the Bayes theorem we can update the initial probability for a model or a parameter after new information is available in the form of experimental observations. The experimental information is fully contained in the so-called likelihood function. This term refers to the conditional probability for the data given the model when it is regarded as a function of the model's parameters while keeping the data fixed to the experimental observations.
We use the following notations: -D = {x i } represents the data, possibly organised in different classes i; -E is the experimental exposure given in terms of the duration time T of the data-taking period and the fiducial mass M det of the detector.
-H r represents a specific hypothesis: H 0 is the background-only hypothesis according to which the known physics processes (backgrounds) are enough to explain the observations; H r S is the background-plus-signal hypothesis for which some DM signal with rate r S is required to explain the data. We note that the two hypotheses H 0 and H r S are nested since H 0 can be obtained for H r S by setting r S = 0. For what concerns our inferential problem of constraining r S , we will always work within the hypothesis H r S , assuming its validity. -r S indicates the expected rate of DM interaction for a given σ SI per unit mass and time expressed in evt/kg/day. It is also a function of the mass m χ of the DM candidate through the cross section. We take as reference cross section the following values: σ SI (ref.) = 10 −41 , 10 −37 , 10 −34 cm 2 for the nuclear recoil, Migdal electron, photon Bremsstrahlung analyses respectively. -r B indicates the rate of the total background events expressed in evt/kg/day; -π(r) is the prior probability density function (p.d.f.) for the generic parameter r and encapsulates all the available knowledge on the parameter r before the experiment is carried out; -L(r; D), or simply L(r) is the likelihood for the generic parameter r of the hypothesis H r , and coincides with p(D | r, H r ); -p(r | D) is the posterior p.d.f. for the generic parameter r given the data D; -θ = (r B , . . . ) is the list of parameters necessary to describe the experimental conditions or theoretical assumptions which are not exactly known but can vary according to they prior p.d.f. π(θ), these are the so-called nuisance parameters as we are not explicitly interested in inferring their posterior values; -Ω is the nuisance parameters space.
We recall that the posterior p.d.f for the parameters of the model can be computed by means of the Bayes theorem as: with The marginal p.d.f. of the parameter of interest r S is given by: and similarly for any other parameter of the model.

90% Credible Interval upper bound
We compute the upper bound for the DM signal as the 90% Credible Interval (C.I.). This is defined as the value of σ SI (m χ ) corresponding to the 90% quantile of the posterior p.d.f. for r S : In the Bayesian approach the upper bound is a statement on the true value of the parameter of interest. The quantity r S (90% C.I.) has to be interpreted as the value below which we believe at 90% probability level the true value of r S lies, given the present experimental information.

Prior choice
It is evident from eq. (13) that the posterior p.d.f. depends on the priors on all parameters. However, we have to distinguish the effect due to the priors on nuisance parameters from the one due to the prior on the parameter of interest. The former has the effect of averaging the posterior over the nuisance parameters space, that is an elegant way of propagating systematic effects on the parameter of interest. In addition, the prior of nuisance parameters is often a parametrization of calibration measurements. The latter, although indispensable to invert the probability and get the posterior, has a degree of 'subjectivity' with potentially a significant impact on the posterior. The prior π(r S ) represents the knowledge on r S before the experiment is carried out, and gets updated by a factor proportional to the likelihood of the observed data. It is a critical term in many respects, and it should reflect the researcher state of knowledge. Especially for searches where the sought quantity is unknown and the search is pushed to the limit of the experimental sensitivity, the input from π(r S ) might have sizable effect on the posterior. The prior has thus to be well justified and the posterior sensitivity to different prior choices needs to be explored. In our case, r S depends on m χ and on σ SI . For m χ we chose a flat prior. For σ SI , to explore the sensitivity of the upper bound to the prior choice, we studied its behaviour for a Migdal-only signal at mχ = 1 GeV/c 2 and with n = 1, 2. We generated a pseudo-dataset from the background template and performed a fit using different priors. We tested, in two possible domain ranges, namely D 1 = [10 −40 , 10 −37 ] cm 2 and D 2 = [10 −44 , 10 −37 ] cm 2 , four different prior choices: a uniform prior, a wider gamma prior with a shape k = 1 and a scale θ = 10 −37 cm 2 , a narrower gamma prior with a shape k = 1 and a scale θ = 10 −38 cm 2 and a uniform prior in log(σ SI /cm 2 ) (we will call this prior "loguniform"). In principle, there is no need to restrict the domain to a finite interval, but both the uniform and the loguniform distributions are not normalizable otherwise. In addition, there could be physical motivations that define a reasonable interval. In our opinion, for our problem, from above the natural constraints come from upper bounds imposed by previous experiments as for example Xenon1T [47] or CRESST-III [39], which for the chosen configuration exclude at 90% C.L. cross sections of the order 10 −38 cm 2 . From below we can use two arguments: the first is that below 10 −44 cm 2 the rate would be dominated by neutrino coherent scattering [77] (the so-called 'neutrino floor'), the second is that the experimental sensitivity does not extend below 10 −40 cm 2 , and will be discussed in Sec. 3.1.3. Therefore the choice of the domain D 1 is driven by physical considerations about a LAr experiment with features similar to DarkSide-50 while the choice of the domain D 2 extents up to the maximum experimental sensitivity that an experiment of this kind can reach before hitting the neutrino floor.
The results, in terms of posterior p.d.f. for r S and σ SI are reported in Figure 3: here we show, as an example, all p.d.f. in the D 1 domain (however the posterior p.d.f. using the domain D 2 are very similar to the one showed in Figure 3 and the differences in terms of the 90% C.I. upper bounds are reported in Table 1). There is no much difference between the uniform and the wider gamma cases, and that is because, as one can see from the left part of Figure 3, these priors are quite flat in the region where the likelihood (the red line in the plot) is mostly informative; on the other side, for the loguniform and the narrower gamma, this is not true in the range where the experiment sensitivity is lost (below 10 −39 cm 2 ), and this is reflected both in the posterior p.d.f. and the 90% C.I. upper bound, as reported in Table 1. We can therefore quantify the dependence of the bound from the prior choice in a factor as big as 10, and this confirms the importance of choosing the prior in a reasonable and coherent way. For simplicity and for reason that would be clear in the Sec. 3.1.3 in the rest of the paper we will report upper bounds

Experimental sensitivity and Bayes factor
A meaningful way to report the experimental sensitivity to the sought phenomenon which is as much as possible independent form the priors is given by the Bayes factor. The posterior p.d.f. for a signal rate of r S given a background rate of r B , and x observed events can be normalized to the posterior for r S = 0, obtaining: where the first factor on the right hand side is called Bayes factor. It is independent on priors and it is simply given by the likelihood ratio of the two hypotheses: .
The properties of the R function have been discussed in great detail for a similar case study in ref. [78]. Here we only mention that R has the probabilistic interpretation of hypotheses belief updating ratio. It is equal to 1 in the limit r S → 0, in this limit the experimental sensitivity is lost and thus the experiment does not change the relative belief. While R → 0 for large r S , where the posterior density for r S vanishes no matter how strong it was before. In addition, the quantity R is used as test statistic in the frequentist approach to limit settings, for details see Ref. [79] (section 39-Statistics) and Ref. [80,81].
In the simple case of a Poisson process of intensity (r S +r B )E, where E is the exposure, and observed counts x the likelihood is proportional to e −(r S +r B )E [(r S + r B )E] x , thus: We evaluated the R function for the same configuration used in the previous section. To explore how R changes when data differ from the expectation (or pseudo-data in this case) due to statistical fluctuations, we consider two additional pseudo-dataset obtained from the previous one by letting pseudo-data fluctuate by ± 1 (Poisson) standard deviation. The results are reported in the left side of Figure 4 (black lines), where the green lines represent the corresponding results for a greater value of the exposure (i.e. 20 ton yr which is roughly 10 3 times the current exposure E = 6786 kg d). From this figure we see that with the current exposure the informative region where the experiment has sensitivity and R → 0 is starting from r S ∼ 10 −39 cm 2 . Any conventional value of r S in this region would be representative of the experimental sensitivity. The statistical uncertainty associated to this value is given by the band contained between the dashed lines that represent the R function computed using a pseudo-data generated varying the expected rate by ±σ. A possible value of r S representative of this region could be such as that R(r S ) = 0.10, which corresponds to a probability update ratio of 10% with respect to the null hypothesis.
As well described in ref. [78], in order to extract any probabilistic statement on r S from R one has to add the information about the prior. We would like to stress that there's a conceptual difference in using the 90% C.I. upper bound or taking the r S such that R = 0.1: the former is the cumulative of the posterior p.d.f. and then it takes into account all the possible values of r S from 0 to r S (90% C.I.) as well as the prior choice; the latter is the likelihood ratio and it is a punctual comparison, namely it takes into account only one single possible value of r S (the one which solves R(r S ) = 0.1), and it is prior independent.
In Table 2 we show how the two methods gives very similar results. For the r S (90% C.I.) method, the results reported in this table are obtained using the procedures described in the next subsection. We also show the results of the projection to higher exposures, reproducing the expected 1/ √ RE scaling, with RE = E/E 0 and E 0 = 6786 kg d.  Table 2: Upper bound results in two different methods for a Migdal-only signal at m χ = 1 GeV/c 2 and n = 1, 2 using pseudo-datasets generated from the background-only likelihood with an exposure E = 6786 kg d.

Expected sensitivity
The expected sensitivity can be quantified either with the upper bound or with the R function using pseudo-data generated with the likelihood under the background-only hypothesis. In the following of this work we quote only upper bounds at 90% C.I. obtained with a flat prior on r S . This choice is motivated by the pragmatic argument that a bayesian limit is easily comparable with results obtained with different statistical strategies. This approach allows also a straightforward procedure to fold in different priors. Pseudo-data allow to explore the sensitivity to a given signal rate for a specific exposure. For the same exposure, we can generate pseudo-data including also the statistical fluctuation of the observed number of events allowing to study the dependence of the sensitivity to the expected statistical fluctuation if no signal were present. In principle with the same technique we could also explore the impact of systematic effects.
To compute the expected sensitivity and produce the results reported in the right column of Table 2, we generated 1024 pseudo-experiments each of which consisting of a set of D = {x i } pseudo-data simulated from the likelihood under the background-only hypothesis. We then compute the posterior p.d.f. and the upper bound r S (90% C.I.). As central value for r S (90% C.I.) we take the sample mean, and as uncertainty the sample standard deviation. An example of the histogram of r S (90% C.I.) is given in the right side of fig. 4.

The tea-lab simulated LAr experiment
In this work we are interested in studying the impact of the Migdal electron and photon Bremsstrahlung to the sensitivity of LAr experiments to light DM. We expect sizable effects for DM masses below ∼ 2 GeV, where the sensitivity of current experiments is lost.
We develop a toy simulation of a LAr experiment that we call tea-lab (Toy Experiment Analysis of Liquid Argon Behaviour) loosely inspired by DarkSide-50. DarkSide-50 is a Liquid Argon Time Projection Chamber (LAr TPC) [82] operated in the LNGS in Italy. The LAr TPC is red-out by Photomultipliers (PMTs) sensitive to the scintillation light produced by the ionizing events in the active LAr target, the so-called 'S1' signal.
The ionization electrons produced at this stage, and surviving the recombination process, are drifted by the TPC electric field to the liquid surface, where they are extracted into an argon gas layer. The electric field in the gas is large enough to accelerate the electrons which excite the argon such to generate a secondary scintillation signal, 'S2'. The lowest threshold is achieved by exploiting the high gain of the S2 signal and corresponds to a number of primary ionization electrons N e − = 4. This result has been achieved by the DarkSide collaboration thanks to a detailed understanding and calibration of the detector response [9,28,[83][84][85][86]. The DarkSide-50 spectra, as given in Fig. 7 of Ref. [28], refer to S2-only events and correspond to a 6786.0 kg d exposure (corrected for the fiducialization cut). For a detected energy E d > 0.05 keV ee , well below the analysis threshold, the LAr TPC is fully efficient [65], thus no efficiency correction is needed.
tea-lab assumes the DarkSide-50 total background spectra and includes some relevant experimental effects as described below.

The likelihood function
The likelihood function represents the connection between the parameter of the model, both theoretical and experimental, with the observed data D. The experimental observable is the number of nuclear recoils x i that occur producing i primary ionization electrons. We factorise the likelihood 6 in three terms: with the following meaning: -L C : describes the probability that in each observation bin i the number of counts is x i . This is assumed Poisson distributed given the expected number of counts λ i . The counts in different bins are taken as independent. Under these assumptions: with: where E = T M det is the experimental exposure. The quantities S i , B i , and LowN e i are associated respectively to the DM signal, to the total background, and to a possible background source at low N e − to account for the excess of events, assumed to be due to some not completely understood experimental effect [28,65], visible in the published spectrum below N e − = 7. -L B : describes how the background template is affected by systematic uncertainties. We don't have the information about how to include the different systematic effects, and a thorough implementation goes beyond the scope of this paper. However, we decided to include bin dependent Gaussian uncertainties with a standard deviation compatible with a simulated sample size of 10 5 events: this corresponds to a 3-6% uncertainty for N e − ≥ 10. We don't account for any systematic effect on the background spectra. For the LowN e background we proceeded differently. We parametrized its contribution in the range 4 ≤ N e − ≤ 7 with a 2-parameter function, and use the parameters to control its contribution ranging from no contribution at all to something similar to what is visible in the DarkSide-50 spectra. This likelihood term is given by with bkgd = {B, LowN e(p 0 , p 1 )}. For the LowN e we introduced explicitly here the dependence on the two parameters used to model it. Here, these parameters are assumed as given, in the next section we discuss how we deal with their uncertainties. -L S : this factor depends on the systematic uncertainties on the signal template S i . We are considering 2 effects: one that parametrises the uncertainty in the emission probabilities as discussed in section 3.2.2, and the other that describes the experimental efficiency to convert the energy of the Migdal electron in primary ionization electrons. L S is given by where S i (f , N max e − , ) represents the probability of the signal to give rise to an event. It depends on the maximum number of ionization electrons that can be in principle produced N max e − and on the efficiency of the production mechanism and detection as it is explained in section 3.2.2. Finally, the parameter f = (f, f val ) controls the contribution of the emission probability due to the inner (f ) and valence (f val ) shell and are used to parametrise the systematic uncertainties associated to their calculation. The central value of the calculation is obtained for (f, f val ) = (1, 1). Independently of N e − and , the contribution of the different shells can be expressed as: where n indicates the electron shell(s) considered.

Simplified treatment of systematic effects
The likelihood function described above is quite general and can be used to parametrise several systematic effects. However, any complete description of such effects requires a detailed knowledge of the detector which is beyond the scope of this work. For this reason, we leave the description of the systematic effects to the experimental collaboration except for a simplistic treatment of the following few relevant effects.
-Background rate normalization. The total background rate is controlled by the nuisance parameter r B . Although it can be predicted by an accurate simulation, we leave it float with a uniform prior and then constrain it to few percent (see Fig. 7) in the fit with high N e − spectrum. -Low N e − excess. The most conservative way to deal with this not completely understood effect is to remove from our fit the region where effect emerges. By setting a threshold N e − = 7, this region is removed. By lowering the threshold to N e − = 4, the LowN e contribution becomes important. For this configuration we explored 2 options. The first, and more conservative option, is to let the fit account for the excess with a signal contribution, and thus weaken the limit. The second is to model this effect with a 2-parameter function and assign it to an unknown background contribution. We decided to report results also with this configuration as it gives the level of sensitivity one may reach if the excess were understood. In this case we assign normal probabilities π(p 0 ) and π(p 1 ) to the parameters as given by the fit and then we marginalise these parameters in the limit computation.
-Contribution of the electron shells: The contribution of the outermost electron shell to the Migdal effect in LAr is affected by large theoretical uncertainties, and the result given in [33] can safely be taken as an order of magnitude estimate (see Sec. 2.1). For this reason we decided to explore its impact by setting f val = 0, 2 in eq. (25). These values correspond to a variation of ±100% around the estimated contribution of the valence electron give by f val = 1. The additional configuration explored, is to consider f val as a nuisance parameter with a flat prior π(f val ) in the range [0,2] and marginalize this parameter in the limit evaluation. For the inner shells (n = 1, 2) we included a gaussian prior π f = G(mean = 1, std = 0.2) to account for the O(20%) theoretical uncertainty in their calculation [33]. -Fluctuation induced by the detector response: The conversion of the Migdal electron energy (E e ) into a number of ionization electrons (N e − ) is a stochastic process which depends on the details of the liquid Ar ionization and excitation processes and on the detector response 7 . To model these fluctuations we proceed with a very crude approximation of what is done in Ref. [28,65,87]. The average number of ionization electron N e − per keV is taken from Fig. 2 of Ref. [65], while the maximum number of ionization electron can be estimated as N max e − = E e /W , where W 19.5 eV [84,88,89] is the effective LAr ionization work function required to produce an electron-ion pair. The final efficiency can be estimated as = N e − /N max e − . Thus the probability of having a certain number N e − of ionization electrons to produce the detected signal is given by the binomial distribution: As a reference the width for N e − = 10(30) is σ = 1.2 (3.8). This probability depends on the energy E e , and thus it is convoluted with the energy spectrum of the Migdal electron emission.

Analysis model implementation
The computation of the posterior p.d.f., even for models relatively simple as those described in the previous section, is often only possible by Monte Carlo integration. The most common way to solve problems of this kind is by sampling the not normalised posterior distribution by a Markov Chain Monte Carlo (MCMC). For our study we used the general analysis framework R [90] and the MCMC algorithm called Gibbs Sampler as implemented in JAGS [91] and interfaced with R in the package rjags [92]. The details of the implementation and the source code for the analysis are publicly available on GitHub [49]. The Monte Carlo simulation gives the unnormalised posterior p.d.f. of the parameters of interest sampled using the Gibbs algorithm. The results reported in this work are obtained with a single Markov chain with 10 5 steps, which is enough to guarantee the chain thermalization.

Sensitivity to Migdal electron and photon Bremsstrahlung
Having computed the rates for the Migdal effect and the photon Bremsstrahlung process, and described the tea-lab simulation, we are able to study the expected sensitivity to low mass dark matter using tea-lab as a case study. We assume no isospin violation (the neutron and proton cross sections are equal). The nuclear recoil contribution was ignored in the Migdal and Bremsstrahlung signal models because it is small if compared with the electron recoils, for masses below 1.8 GeV/c 2 . For our simulated tea-lab experiment we generated, as already explained in Sec. 3.1 and 3.2, a dataset from the background-only template, including this time the LowN e excess. We take always as a reference exposure E = 6786 kg d.
In order to validate our analysis, we compute the bounds for the nuclear recoil for a DM mass in the range 1.8-15 GeV/c 2 and compare with the published results of the DarkSide collaboration [28]. We compute the bounds with our bayesian approach reporting the (90% C.I.) lower bound for both thresholds at N e − = 4 and N e − = 7 electrons. The result is shown in Fig. 5 and 6 (brown dashed line), together with the DarkSide-50 constraint (solid red line). Here we also validate our simplistic binomial model (brown dotted line) for the experimental response fluctuation. Our result is in good agreement with the experimental bounds, considering also the fact that the latter are calculated using the frequentist approach known as CL s [80]. Fig. 5 and 6 show the impact of the Migdal effect and the photon Bremsstrahlung for tea-lab in extending the sensitivity region of LAr experiments from m χ ∼ 2 GeV/c 2 to m χ ∼ 0.05-0.08 GeV/c 2 . For m χ ∼ 1 GeV/c 2 the sensitivity based on the Migdal effect is σ SI ∼ 10 −37 cm 2 ; for m χ 0.11 GeV/c 2 the sensitivity is comparable with the Xenon1T exclusion limits [47] and extends up to m χ ∼ 0.05-0.08 GeV/c 2 . For the photon Bremsstrahlung the sensitivity spans in the range m χ ∼ 0.12-2 GeV/c 2 , exploring a region inaccessible to Xenon1T.
In addition, we estimate the effect of the Earth attenuation, as discussed in Section 2.4. We assume for the tea-lab estimate a shielding due to the Earth crust compatible with the LNGS underground cavern. Our upper limit estimate corresponds to the largest cross section that can be probed by tea-lab using only kinematics assumptions. This implies that after a thorough accounting of background contributions, experimental effects, and signal time modulation, the corresponding upper bound might result below the one quoted here. In particular, we computed the bound for the two thresholds under consideration, shown as a black dashed curve in Figures 5 -9, 11 and 13. In these figures the gray shaded area, defined from below by the Migdal 90% C.I. upper bound, and from above by the Earth attenuation lower bound, has to be intended as the region where a generic LAr experiment as simulated in tea-lab has the sensitivity to exclude a DM signal exploiting the Migdal signal 8 .
In the next subsections we discuss in details our results exploring the impact of the various systematic effects.

Impact of theoretical uncertainties
As a first step of our study, we focused on the determination of the impact of the theoretical uncertainties on the outermost shell contribution to the Migdal cross section. As already explained in Sec. 3.2.2, since these are O(1) uncertainties, we decided to compute the limit in four possible scenarios: assuming f val = 0, 1, 2, or considering f val as a nuisance parameter with a flat prior in the range [0, 2]. In order to isolate this effect, we produced a tea-lab simulated dataset based on the background-only template, without the LowN e excess.
To make sure that the introduction of additional nuisance parameters does not create The tea-lab bound for the photon Bremsstrahlung signal is reported (orange line). The estimate of the Earth shielding effect for tea-lab is also reported (black dashed curve). The gray shaded area, has to be intended as the region where tea-lab has the sensitivity to exclude a DM signal exploiting the Migdal electron. The upper limits of the Xenon1T [47], CRESST [39], and DarkSide-50 [28] are reported. As a cross check of the tea-lab simulation we report also our calculation for the DarkSide-50 bounds on the NR signal with binomial fluctuation (brown dotted line) and without (brown dashed line).
instabilities in the fit, we studied the global posterior p.d.f. which is represented in Fig. 7. The sensitivity results are reported in Fig. 8 and 9 where we considered spectra starting from N e − = 4 and N e − = 7, respectively. With respect to sensitivity with N e − = 4 threshold, for m χ 200 MeV/c 2 the effect of the valence shell is most clear: if f val 1 the sensitivity is lost due to the fact that for small masses the greatest contribution to the Migdal signal comes from the outer shell, as depicted in Fig. 10, where we show the Migdal The tea-lab bound for the photon Bremsstrahlung signal is reported (orange line). The estimate of the Earth shielding effect for tea-lab is also reported (black dashed curve). The gray shaded area, has to be intended as the region where tea-lab has the sensitivity to exclude a DM signal exploiting the Migdal electron. The upper limits of the Xenon1T [47], CRESST [39], and DarkSide-50 [28] are reported. As a cross check of the tea-lab simulation we report also our calculation for the DarkSide-50 bounds on the NR signal with binomial fluctuation (brown dotted line) and without (brown dashed line).
signal for a mass of 130 MeV/c 2 with respect to the Migdal signal for a mass of 1 GeV/c 2 . As f val increases, the contribution of the third shell becomes more important, and therefore the bound becomes much stronger, even of an order of magnitude. For m χ 200 MeV/c 2 this is no more true because the contribution of the inner shells are now of the same intensity of the contribution of the outer shell. Then, for m χ 200 MeV/c 2 , we can assert that our limit is solid against systematic uncertainties on the contribution coming from the valence shell to Migdal signals in LAr. For the same reasons, if we consider spectra starting from N e − = 7 the contribution of the outer shell is almost completely cut off, and the dependence on f val is weakened, leading to a departure of the sensitivity in the various cases from m χ 100 MeV/c 2 . With respect to the photon Bremsstrahlung signal, the comparison between the two figures shows that lowering the threshold does not produce a significant change in the sensitivity. The same considerations hold for Fig. 5, 6 in which the overall sensitivity is reduced, due to the presence of the LowN e excess, by a factor as big as ∼ 3.

Impact of the experimental effects
As we have seen in the previous paragraph a threshold N e − = 4 allows to get an expected sensitivity to masses down of 0.06 GeV/c 2 which is much stronger that the one obtained with a higher threshold. To study the impact of LowN e excess in more specific way, as detailed in Sec. 3.2.2, we used two strategies. The most conservative is to let the fit The tea-lab bound for the photon Bremsstrahlung signal is reported (orange line). The estimate of the Earth shielding effect for tea-lab is also reported (black dashed curve). The gray shaded area, has to be intended as the region where tea-lab has the sensitivity to exclude a DM signal exploiting the Migdal electron. The upper limits of the Xenon1T [47], CRESST [39], and DarkSide-50 [28] are reported. account for the excess with a DM signal in the limit calculation procedure. This results in a lower bound on DM and corresponds to the blue lines in Fig. 5. We produce results also in the case where the excess would be understood and modeled as a background. To emulate this circumstance we fit the excess as described in Sec. 3.2.2 and assign it to the background component. The new bounds are plotted in Fig. 11 as a solid blue line (for f val = 1, i.e. including the contribution of the valence shell) and yellow line (for f val = 0). For comparison the dashed lines give the limit in the default scenario where the excess is not accounted for as a background. We notice that, if the contribution of the valence electrons is neglected, there is no difference between the two approaches. This is because the inner shells do not give any signal in the region where the excess is present (see Fig. 10). On the contrary, since the valence shell contribute significantly to the region where the excess appears, accounting for it with an additional background contribution leads to a stronger bound, roughly independent of the mass. In Fig. 11 we also show how the sensitivity bounds would improve extending our analysis to N e − ≥ 2, always assigning The tea-lab bound for the photon Bremsstrahlung signal is reported (orange line). The estimate of the Earth shielding effect for tea-lab is also reported (black dashed curve). The gray shaded area, has to be intended as the region where tea-lab has the sensitivity to exclude a DM signal exploiting the Migdal electron. The upper limits of the Xenon1T [47], CRESST [39], and DarkSide-50 [28] are reported.
the LowN e excess to the background component. Finally, as for Sec. 4.1, for a threshold N e − = 7 the dependence on f val is weakened. Fig. 12 shows also the impact on the signal spectra due to the binomial fluctuation of the detector response. For a NR signal at m χ = 1.8 GeV/c 2 , it is clearly visible how these fluctuations let some of the events spill above the N e − = 4 threshold. This effect is not present for the Migdal electron spectra since they extend well above the analysis threshold independently of the fluctuations model. We checked that including the binomial model does not change the results for the Migdal electron and photon Bremsstrahlung.
In conclusion, we point out that to exploit at the best the contribution of the Migdal effect and the photon Bremsstrahlung, it is crucial to have a precise description of the backgrounds in the low N e − region.

Projected sensitivity for future experiments
To illustrate how the sensitivity to DM particles enhanced by the Migdal electron and the photon Bremsstrahlung scales with the experimental exposure, we show in Fig. 13 the tea-lab 90% C.I. bounds for an exposure E = 5 ton yr. This exposure corresponds roughly to an increase of a factor RE = 270 with respect to the results presented in Fig. 8 for the exposure E = 6786 kg d. As expected, the improving factor in the sensitivity is about 1/ √ RE 16, which corresponds to the square root of the increased exposure. For this projection we also scaled the uncertainties by a factor of 1/ √ RE in order not to be dominated by systematic effects. We point out that the best achievable improvement with the exposure is attainable only if the bounds are limited by the sample size. When the statistical uncertainty becomes comparable with the systematic uncertainty, the limit stops improving with an augmented data sample. This effect is particularly relevant for the background subtraction. In fact, already with an exposure factor RE 10, the limits become dominated by systematic effects due to the large contribution of the backgrounds and their limited knowledge. This effect can be canceled up to a certain RE by having a more precise knowledge of the background contamination; however, already at RE = 270, the background uncertainty has to be smaller than 1‰. This clearly imposes the necessity to define strategies to reduce the background contamination to push the limit down with large exposures. Fig. 13 shows also the neutrino floor for a liquid argon experiment as given in ref. [28]. For m χ ∼ 6 GeV/c 2 the projected sensitivity for the nuclear recoil spin independent signal starts to approach the neutrino floor. For the Migdal effect this is not the case as we are few order of magnitude above. However, we point out that it will become important to have a reliable estimate of the neutrino floor in liquid argon experiments for The estimate of the Earth shielding effect for tea-lab is also reported (black dashed curve). The gray shaded area, has to be intended as the region where tea-lab has the sensitivity to exclude a DM signal exploiting the Migdal electron. The upper limits of the Xenon1T [47], CRESST [39], and DarkSide-50 [28] are reported. masses below m χ 0.5 GeV/c 2 .
Bearing these considerations in mind, it appears clear that any reasonable extrapolation of our result to higher exposure as for example DS20k [93], which is foreseen reaching values as big as E = 3000 ton yr with a completely redesigned detector, can only be done by the experimental collaboration after a thorough assessment of the relevant systematic effects. Figure 12: Migdal electron signal for m χ = 1 GeV/c 2 (blue lines) and m χ = 1 GeV/c 2 folded with the binomial model for the detector response (orange lines), both evaluated for σ = 10 −37 cm 2 , considering only n = 1, 2. The dashed lines represent the same signals including also n = 3. NR for m χ = 1.8 GeV/c 2 and σ = 10 −41 cm 2 with the binomial model for the detector response (green line), without (brown line). Figure 13: 90% C.I. upper bounds on the σ SI exploiting the Migdal electron and photon Bremsstrahlung signals for the tea-lab simulated experiment loosely inspired by DarkSide-50. These results are obtained using and exposure E = 5 ton yr, and simulating a pseudo-dataset with the background-only template, not including the LowN e excess (see Sec. 3.2, 3.1). The tea-lab bounds are computed for a threshold N e − = 4. The estimate of the Earth shielding effect for tea-lab is also reported (black dashed curve). The gray shaded area, has to be intended as the region where tea-lab has the sensitivity to exclude a DM signal exploiting the Migdal electron. The upper limits of the Xenon1T [47], CRESST [39], and DarkSide-50 [28] are reported. We report also our calculation for the DarkSide-50 bounds on NR signal extrapolated at E = 5 ton yr (brown dashed line). The neutrino floor for LAr experiments [28] is reported as a dark-gray shaded area.

Conclusions
We considered both the Migdal effect and photon Bremsstrahlung from the nucleus hit by a DM candidate in experiments based on LAr detectors. Previous DM searches exploiting LAr detectors include DarkSide-50 [9] and DEAP-3600 [12], and have mainly focused on masses greater than 10 GeV/c 2 by considering only the nuclear recoil signal. We decided to develop the tea-lab simulated experiment taking inspiration from the DarkSide-50 experiment, for which there is a published low-mass analysis [28] sensitive down to masses of 1.8 GeV/c 2 . Since at low masses relative large DM cross sections are probed, we estimated the effect of the Earth attenuation, which results in a lower limit on σ SI above which the experiment becomes blind.
In order to study the sensitivity to low masses exploiting the Migdal effect and photon Bremsstrahlung, we built a simplified likelihood describing the experimental response of tea-lab. In particular we payed attention to background components, to the description of the low N e − excess, and to the detector response including the fluctuation in the number of ionization electrons. We also considered the systematic uncertainties on the theoretical computation of the Migdal spectra, among which the contribution of the valence shell of argon plays a critical role, being estimated with an uncertainty of O(1). We therefore used this likelihood in the context of a Bayesian approach to compute the expected sensitivity, the exclusion limits at the 90% C.I., and the projections at higher exposure. We show that the limits are mildly dependent on reasonable prior choices describing the present state of knowledge on the sought signals.
As a preliminary cross check, we computed the limit for the nuclear recoil signal using the DarkSide-50 measured spectrum obtaining a result compatible with the expected one.
We then studied the tea-lab expected sensitivity at the DarkSide-50 exposure and show the impact of the systematic effects on the final result. The sensitivity spans masses in the range 0.1-2 GeV/c 2 . For m χ 0.2 GeV/c 2 , the results are stable and insensitive to significant theoretical uncertainties on the argon valence shell contribution to the Migdal effect. For m χ 0.2 GeV/c 2 instead these effects can substantially modify the limit, even by two orders of magnitude for the lowest masses (see Fig. 8).
We studied the interplay between the poorly understood LowN e excess and the Migdal effect showing that the bounds change only if the contribution of the valence shell is considered. In this case, a better understanding of the excess results in a stronger limit. For m χ ∼ 0.11 GeV/c 2 the limit is comparable with that of the Xenon1T experiment obtained with a much larger exposure, and it reaches m χ ∼ 0.06 GeV/c 2 extending the sensitivity of noble liquids to lower masses.
Finally we projected the sensitivity analysis up to an exposure E = 5 ton yr, finding that the current knowledge of the background systematic uncertainties prevents the limit to scale with the statistics for exposure larger than ∼ 0.2 ton yr. This shows how any reasonable extrapolation of our results can only be done by the experimental collaboration after a thorough assessment of the relevant systematic effects of the foreseen experiment.
We emphasize that the large uncertainty on the contribution of the valence shell for the computation of the Migdal effect already plays a crucial role in setting a bound from argon detectors for DM masses below 0.2 GeV/c 2 . As a consequence, a dedicated theoretical effort to compute the transition probabilities with higher precision is needed.
We conclude by stressing that the contribution of a Bremsstrahlung photon and a Migdal electron is sizable also for LAr experiments, and it is a powerful tool to explore mass regions below the GeV/c 2 scale, inaccessible by modeling the signal only with the traditional nuclear recoil interaction. We look forward to an update of low masses analyses in LAr experiments that includes such effects.