Isospin-violating dark matter at liquid noble detectors: new constraints, future projections, and an exploration of target complementarity

There is no known reason that dark matter interactions with the Standard Model should couple to neutrons and protons in the same way. This isospin violation can have large consequences, modifying the sensitivity of existing and future direct detection experimental constraints by orders of magnitude. Previous works in the literature have focused on the zero-momentum limit which has its limitations when extending the analysis to the Non-Relativistic Effective Field Theory basis (NREFT). In this paper, we study isospin violation in a detailed manner, paying specific attention to the experimental setups of liquid noble detectors. We analyse two effective Standard Model gauge invariant models as interesting case studies as well as the more model-independent NREFT operators. This work demonstrates the high degree of complementarity between the target nuclei xenon and argon. Most notably, we show that the Standard Model gauge-invariant formulation of the standard spin-dependent interaction often generates a sizeable response from argon, a target nuclei with zero spin. This work is meant as an update and a useful reference to model builders and experimentalists.


I. INTRODUCTION
After decades of dedicated searches, the nature of dark matter has only been probed through its gravitational interactions.Although not necessary, many theories produce the required abundance of dark matter in the Universe through interactions with the Standard Model (SM) which are much stronger than gravity.Many such theories produce dark matter through thermal freeze-out, which has appealing features such as insensitivity to initial conditions as well as the prediction of detectable signals for experimental searches, such as colliders, telescopes, and direct dark matter detectors.So far, nature seems more complicated, and large regions of parameter space for the simplest particle models of dark matter are highly constrained [1,2].
Nevertheless, the next generation of direct dark matter detection experiments is poised to make further progress in the next decade [3][4][5][6].At the same time, the theoretical community has shown that there are numerous possibilities that are still viable.As progress is made both experimentally and theoretically, it is timely to ensure we understand any potential blind spots in our search techniques.Sometimes these blind spots are difficult to circumvent, such as p-wave annihilation for indirect detection searches and momentum or velocity-suppressed scatterings in direct detection.In this paper, we focus on the possibility that interactions in direct detection experiments are suppressed by a conspiracy of couplings between dark matter and quarks, such that nucleons no longer cumulatively sum, and instead cancel.This can happen when couplings between the dark matter and neutrons and protons (c n and c p ) are not equal and is known as isospin violation [7].
The phenomenological motivation for considering isospin violation is clear, there are numerous theoretical models of dark matter, many of which violate isospin [8][9][10][11][12][13][14][15][16].One simple situation where the isospin-conserving case happens is when dark matter is coupled to SM quarks universally via a pure vector interaction.This is a specific case, and for a spin-1 mediated model, there are a sufficient number of free parameters to achieve any relation between c n and c p , the couplings of dark matter to the neutron and proton respectively.However, for UV-consistent models that do not contain anomalies and therefore indicate some unitarity violation, the situation is often more complicated [17][18][19][20], leaving one with a limited number of simple models.Amongst them, is the U (1) Lµ−Lτ model with a fermionic dark matter candidate.Here, direct detection with nuclei proceeds from kinetic mixing between the photon and the U (1) Lµ−Lτboson.Any vector interaction with the photon results in c n ̸ = c p since their electric charges are different.
Isospin violation in direct detection has been explored previously [21][22][23][24], and most commonly in the context of spin-dependent (SD) and independent (SI) interactions.Such interactions are independent of incident velocity v and transfer momentum q and therefore it is often a good approximation to estimate the effect of isospin violation in the q, v → 0 limit.This is obviously not possible for interactions that are v and q dependent.In this paper, we use the non-relativistic effective field theory (NREFT) basis [25][26][27][28] to explore isospin-violating effects on these momentum and velocity dependent operators.
Through analysis of scintillation and ionisation atomic responses to dark matter scattering, the liquid xenon and argon-based TPC experiments XENON [4,29,30], LZ [6,31,32], PandaX [33][34][35], DEAP-3600 [36] and DarkSide-50/20k [3,[37][38][39] have recently provided or will in the future provide the greatest sensitivity to dark matter for a large number of possible candidates in the GeV-TeV regime.To study the effect of isospin violation, we carry out a realistic modeling of a representative subset of current and future xenon and argon experiments -specifically XENON1T, DarkSide-50, DEAP-3600, LZ, and DarkSide-20k -using publicly available data, and from this are able to reappraise the current constraints on dark matter candidates and the future prospects for detecting isospin-violating dark matter in direct detection experiments.We evaluate the importance of certain experimental inputs such as detector efficiency on the levels of signal suppression and exemplify how the signal region can affect the interpretation of spectral shapes that arise in isospin-violating scenarios, which could lead to operator misidentification.The main results of this paper are our projections for the coming generation of argon and xenon detectors, how their limits are affected by isospin violation, and how the two targets will complement each other.We present our results by way of multiple case studies, two of which are motivated by gauge invariant dark matter models that lead to the canonical spin-independent (SI) and spin-dependent (SD) interactions.We show that, somewhat surprisingly, argon will be very competitive in the SD case.Furthermore, we present results for numerous NREFT operators in a way that accounts for the effects of m χ , such that others can use this work to estimate how other configurations of c n and c p will affect detectability.
This paper is organised as follows, in Section II we introduce the general formalism used to calculate predicted signals of particle dark matter in direct detection experiments as well as the two specific models we focus on.In Section III we describe the setups we use to emulate existing and future direct detection experiments, we also describe the statistical methods we employ to generate our results, which are presented in Section IV.The results are separated into two parts, subsection IV A shows our results for the spin-(in)dependent models we consider, and subsection IV B describes how to interpret our model-independent results on the NREFT basis: a comprehensive suite of results are shown in the Appendix.Finally, we conclude in Section V.

II. THE NON-RELATIVISTIC EFT AND GAUGE INVARIANT MODELS
In direct detection, the incident dark matter from the halo is non-relativistic and one can use the aforementioned NREFT.In this approach, one takes the nonrelativistic interaction Lagrangian, where, χ and N are the dark matter and nucleon fields respectively [26].The O i is a Galilean operator built from 1, the transverse initial velocity of the dark matter, v ⊥ , transfer momentum, q, and the spin of dark matter and nucleon, S χ and S N respectively.A non-exhaustive list of these operators is shown in Table I.The coefficients c i are model dependent and determined by the non-relativistic limit of whatever fundamental interactions take place between dark matter and quarks [40].
For these coefficients, we take the convention that normalises them with respect to the Higgs vacuum expectation value, ⟨v⟩ Higgs = 264 GeV, as in Refs.[27,41,42].The list in Table I only shows operators that we explicitly consider in this work.Additional operators have been derived and studied, see Refs.[27,28,40,43], but we have decided not to include them in this work.Table I covers all operators that arise for spin-0 and spin-1/2 dark matter models with scalar, fermion, or vector mediators.For the operators above O 11 , known scenarios include a spin-1 dark matter candidate which is mediated by a vector or charged fermion.Often these operators are present along with those shown in table I and are frequently subdominant.Nonetheless, they can be important for specific cases [43] and a study of operators beyond O 11 will be the subject of future work.The differential recoil rate, w.r.t.nuclear recoil energy E R is given by, where ρ 0 is the local dark matter density, m T is the target nuclei mass, and v is the velocity of dark matter, which is distributed according to the halo distribution f (v).The integral limit v min refers to the minimum velocity required to induce a recoil of energy E R .The differential cross-section, dσ χT /dE R is related to the coefficients in Eq. ( 1) by where are the nuclear form factors.We follow the convention of Ref. [26], where the operator interference is explicitly included.These are decomposed into specific nuclear response functions which have been parameterized and made publicly available [26,27,44,45].
Although the NREFT approach allows for a powerful model-independent analysis of direct detection data (see Refs. [46][47][48][49] for recent experimental analyses), specific models only require a small number of them to be treated.Additionally, due to the hierarchy in scales for certain operator responses, and the ubiquity of some operators (O 1 , for example), not all operators are of equal interest theoretically.For this reason, we focus on two effective particle models that facilitate isospin violation in a realistic and general way.Namely, we consider a fermionic dark matter particle with either vector-like couplings ( χγ µ χ) or axial vector-like couplings ( χγ µ γ 5 χ).Of course, both can be present simultaneously, but we consider them separately because the former leads to the dominant O 1 response and if dark matter is Majorana, the later is the only non-vanishing vector current.
We want to ensure that SM gauge invariance is always manifest, in order to do that we follow the parametrisation of Refs.[50,51], where we start with the relativistic effective operators where Q L , u R and d R are the quark left-handed doublet, up-type and down-type right-handed singlets respectively.The index i denotes the quark generation.
The operator numbering we use follows ref. [51], which contains the dimension in the superscript (for this paper always 6), the generation i, and an arbitrary number that labels the operators by their form.After electroweak (EW) symmetry breaking, the three operators above match onto 1,q = ( χγ µ χ)(qγ µ q) and Q 3,q = ( χγ µ χ)(qγ µ γ 5 q) (5) to ensure that the Lagrangian after EW symmetry breaking remains gauge invariant under the SM, one has to ensure that the couplings in the chiral basis of Eq (4) gives the same results as that of the Dirac basis in Eq (5).The matching conditions are [18,52,53] and and where the C's and C's represent coefficients to the effective operators before EW symmetry breaking, Eq. ( 4) and after, Eq. ( 5), respectively.With this matching, the four coefficients that are relevant to direct detection contain some redundancy.We eliminate this redundancy by following the parametrisation of Ref. [50], which focused on the dark axial-vector interaction.Here we do it for dark vector-like couplings as well 3,1 → 4,1 → g ′2 Λ 2 sin θ sin ϕ.
The g ′ and Λ are introduced to make the connection with the high-scale model that produces the effective operators of Eq. ( 4), g ′ is the gauge coupling and Λ is the energy scale of the new model.Since all processes in this work are assumed to be well below Λ, we treat g ′ /Λ as a combined parameter that sets the overall strength of interactions.The angles θ and ϕ generically parameterize the relative strengths of the three interactions.By adopting this parametrisation we are explicitly keeping interactions SM gauge-invariant and maintaining some connection with the full extension of the SM that explains DM.When one performs the mapping of these coefficients to the NREFT coefficients in Eq. ( 1), they obtain contributions to the NREFT operators O 1 , O 7 , O 9 .Typically because the coherent nuclear response to O 1 is so large, one ignores the other contributions.Since this paper is concerned with cancellation effects, it is conceivable that there are areas of parameter space where O 7 and O 9 are important, perhaps limiting the isospin suppression.
In the top row of Figure 1 we show the relative size of the differential recoil rate Eq. ( 2) for a xenon and argon target respectively.The rate depicted by the colour map is normalized simply to the highest value in the map itself.We do so according to the parametrisation in Eq. ( 7) and take the E R → 0 limit.It can be observed that there are substantial regions where there is a large suppression < ∼ 10 −1 of the differential rate.Furthermore, when one compares the suppression plot for different target nuclei, i.e. xenon and argon, the pattern is similar.This is because where the suppression occurs is largely dominated by O 1 and is extremal when c n /c p ≈ −Z/(A − Z), where A and Z are the use atomic mass and proton numbers respectively.Note that the minimum of the colour map is not the true minimum, but is instead capped at 10 −3 for aesthetic reasons.
Similarly, for the dark axial vector interaction, the prescription above is followed but replace ( χγ µ χ) → ( χγ µ γ 5 χ), and following the nomenclature of Ref. [51].As before, and exactly as in Ref. [50], we can eliminate any redundancies by parameterizing the coefficients to these operators Performing the NREFT mapping, one ends up with non-zero contributions for O 4 , O 6 , O 8 , and O 9 .Similarly to the top row of Figure 1, we map out the suppression of the SD interactions within this parametrisation in the  7) for 100 GeV dark matter mass.The top row reflects results for the spin-independent vector-like couplings model (where the contribution of the O1, O7, and O9 operators is considered) while the bottom row represents the situation in the spin-dependent axial vector-like couplings model (with operators O4, O6, O8, and O9).In the spin-independent case the contour of maximal suppression corresponding to a fixed neutron-to-proton coupling ratio is highlighted by a dashed line.The star and circle markers correspond to specific model points studied and detailed in Section IV.
bottom row.Here the situation is substantially different, the suppression pattern for xenon and argon no longer resemble each other.For xenon, the condition of gauge invariant interactions provokes other unavoidable nuclear responses that diminish the suppression in many cases, as shown in Ref. [50].This is because the event rate associated with the various contributing operators are similar in size, unlike in the SI case.For an argon target, the pattern we see is much simpler because in the argon target only the 40 Ar isotope is present in large quantities and it does not yield any response to the O 4 , O 6 , or O 9 operators.The corresponding patterns then show a high degree of complementarity between the two targets for the SD model.Since the plots of Figure 1 are normalised differential rates, the actual complementarity will depend on how sensitive xenon and argon target experiments will be to the unnormalised differential rate.The main point of showing Figure 1 is to exhibit where suppression occurs in the gauge invariant parametrisation, in a way that is as detector independent as possible.In later sections, we show how these parameter choices affect interpretations of current and future experimental limits.
Furthermore, in order to minimize the experimental input of Figure 1 we have chosen to show the differential recoil rate in the limit E R → 0. This means that q-dependent operators such as O 9 do not affect the patterns shown.It should be noted that other works that investigated isospin violation simply took the suppression in the E R → 0 limit and use such suppressions to recast calculated limits [8,10,12,13,15,22].Of course the presence of q-dependent operators means that this limit can no longer be taken.Additionally, doing this overlooks some interesting features of isospin-violating signals.To our knowledge, the most careful treatment to Recoil spectra in xenon (left) and argon (right) for varying interaction types, dark matter masses and neutron-to-proton coupling ratios, illustrating how isospin-violating interactions can induce large spectral changes in addition to overall rate.The grey-shaded region highlights regions typically inaccessible in the experimental setups considered here.
date is found in Refs.[54,55] where the authors do consider q and v-dependent operators and explore isospin violation in a model independent way.Ref. [54] specifically looks at the c n /c p values that maximally suppress direct detection signal and how it varies with m χ , while Ref. [55] employs a statistical approach to derive limits that account for isospin violation.This article builds on the existing literature by showing how the extent of the suppression varies with m χ as well as providing useful reference values for anyone with an isospin-violating model, our work also explicitly shows the important role argon-based detectors will have in this regard.
Figure 2 shows multiple spectral shapes for different NREFT responses, for a xenon target (left panel) and argon target (right panel).Because of the different strengths of nuclear response to certain operators, we show the differential rate normalised such that the integral over dE R is 1, as this better illustrates the spectral variations.We also show roughly where the expected experimental thresholds will be, E th , by plotting grey shaded regions for E R ≤ E th .Let us first examine why moving beyond the q → 0 limit, even for momentumindependent operators such O 1 (blue lines), may be important.We have chosen to show the isospin conserving case (solid) along with a result that is close to the maximum suppression for xenon, c n /c p = −0.7 (dashed).What we want to emphasise is that the isospin-violating spectrum goes to zero at significantly higher E R , even for the same values of m χ , and does not exhibit the characteristic exponential shape that would be expected for a typical O 1 interaction.Also notice how the rate for the c n /c p = −0.7 case in xenon is larger below E th , under-lining the potential importance that detector effects may have on the expected level of suppression.
The fact that the c n /c p = −0.7 spectrum for O 1 in the signal region is not monotonically decreasing is notable because this behaviour is not expected for q-independent recoils.We exhibit the behaviour of the momentumdependent operator O 11 with c n /c p = 1 for comparison.We have deliberately chosen a value of m χ such that the spectra start to resemble each other.Exemplifying that this can also occur in velocity-dependent operators, we include the O 8 spectrum, which much like O 1 , is monotonically decreasing in the isospin conserving case.However, when c n /c p = −0.7 we see that the spectrum matches very well with O 11 , even when the values for m χ are really very different.We have additionally included the c n /c p = −0.74 with O 11 case because we believe that this spectral shape is quite unique.Potentially there are inelastic dark matter models that would replicate this feature.
Other than underlining the importance of doing spectral analysis when determining how much suppression one expects, Figure 2 shows how, in the event of signal detection, specific isospin violating configurations may lead to signal misidentification.The right panel in Figure 2 shows the exact same benchmarks, but now the recoil spectra is with an argon target.We see that the effect on the signal shape is rather small.This is because the suppression for argon peaks at a different value of c n /c p .Post-discovery, once the overall rate of an observation is also taken into account, correct operator identification should be possible with both a xenon and argon target experiment with comparable and requisite sensitivities.This stresses the importance of complementarity between xenon and argon detectors that is present even in the SI particle model.
The non-standard signal shapes in Figure 2 come directly from the nuclear form factors and how they interfere at different recoil energies, see Eq. (3).Therefore, the consequences of isospin-violating models are highly dependent on nuclear responses.In order to check that our results will be in some way robust to any uncertainty in these nuclear form factors, we would require a systematic way of incorporating them into our analysis.However, as far as we are aware, this has not been made available, at least for the NREFT basis, see Ref. [56] for a discussion of SD structure functions.We leave a systematic study to future work, and instead present here illustrative examples of how, in the regions of parameter space where rate suppression from isospin violation is important, the choice of nuclear form factor has a large impact on the expected signal.
In Figure 3, we compare the O 1 form factors from Refs.[26,44], which we name F Fitz. and F Hof. after the first authors of each reference respectively.We are unaware of any reference where such a large difference has been pointed out.Recent works that use the F Fitz. form factors are [42,46,48,[57][58][59] whereas works like Ref. [49] use F Hof. .Other works such as Refs.[60,61] and popular codes such as DDcalc [62], WIMPy [63] and WimPyDD [64] use Ref. [27], which we find much more closely resembles that of F Hof. than F Fitz. .Figure 3 shows the differential recoil rates for the O 1 interaction for xenon (left) and argon (right) targets.We see that in the zero momentum limit, the form factors are in agreement.At higher E R , they start to disagree, underlining the importance of both having a reliable understanding of nuclear responses and determining the isospin suppression within the signal region and not just in the E R → 0 limit We can see that the differential rates for a xenon target are quite sensitive to the choice in form factor.The shaded region, highlighting the difference between F Fitz. and F Hof. , can be very large, even when the couplings are not so close to the maximal suppression of c n /c p = −0.9.Notice as well that the F Hof. form factors tend to be more suppressed over the whole range than F Fitz. for xenon.For argon, the simpler nucleus, the variation in rate is negligible in the isospin-conserving case.Additionally, when the variations are large as with the c n /c p = −0.84 and c n /c p = −0.9,we see that the effects tend to start at higher E R .In Appendix 1 we show how isospin suppression is affected by the choice of form factors and the effect can be substantial.In Appendix 1, we observe that the largest difference between form factors occurs for O 1 when isospin suppression is most pronounced.This is because the behaviour in the limit as E R → 0 no longer dominates over the signal.Moreover, since isospin suppression involves a fairly fine cancellation between form factors, and over a range of E R , it makes sense that the regions of greatest suppression are particularly sensitive to changes in the way form factors are parameterised.
It is beyond the scope of this paper to speculate on the origin of these differences, but given the size of the effect, we wanted to point it out.We are aware that considerable work is being undertaken currently to independently verify the accuracy of the nuclear form factors and to determine any sources of theoretical uncertainty [44,45,65].Furthermore, we see this as further motivation for the state-of-the-art techniques that accurately model nuclear responses as well as encouraging a further investigation of the uncertainties associated with such calculations.
In terms of accuracy, we believe that Refs.[44,45] are the current gold standard for SI interactions, and we will continue with our study using these.We take this position because Ref. [44] states that they have made use of improvements in nuclear modelling, including multinucleon interaction effects, to provide more precise form factors. Additionally, authors of Ref. [44] go to much greater lengths to validate their nuclear structure calculations by comparing nuclear spectra of their shell model to experimental results.For the SD responses, we use the form factors in Ref. [27].

III. EXPERIMENTAL SETUP AND LIMIT EVALUATION
In this section we outline our derivation of bounds and projections for the direct detection experiments we consider.As mentioned above, we focus on xenon and argon targets, which are the basis for what are or will be the most sensitive experiments for m χ > ∼ 10 GeV when one considers S1 (scintillation) and S2 (ionisation) analyses.Currently, the most recent argon experiments, DarkSide-50 [37] and DEAP-3600 [36] are substantially less sensitive than XENON1T [66] due to either differences in exposure (0.046 tn yr for DarkSide-50 versus 1 tn yr for XENON1T), or due to efficiency in the case of DEAP-3600.However, the argon detector that will be taking data in this decade, DarkSide-20k [3], with its projected 200 tn yr exposure, will be much more competitive with LZ [31] and XENON-nT [4] that are expected to have 15.33 and 20 tn yr exposures respectively.It is this set of experiments in which we are basing our projections for future-Xe and future-Ar experiments.We do not attempt to project sensitivities for next-generation experiments similar to that of ARGO [67] or DARWIN/XLZD [5] because they are at an earlier stage of their development and many experimental parameters and characteristics are yet to be determined or finalised.
We perform our analysis in terms of recoil energy and do not simulate detector responses in terms of S1 and S2 signals.We, therefore, do not replicate the full S1+S2 two-dimensional likelihood calculation performed by the experimental collaborations.To obtain the number of expected recoil events, N , from the differential rate, Eq. ( 2), one needs to include detector specific information, such as exposure, ε, the energy-dependent acceptance ϵ(E R ) and detector resolution.In particular, the resolution is modeled as a Gaussian and is accounted for via a convolution, Additionally, one might need to include an acceptance cut, which typically is some multiplicative factor and is sometimes already incorporated in the efficiency.
For XENON1T we make use of the public code associated with Ref. [68] to carry out an exclusion limit calculation in terms of the reconstructed energy, which approximates the XENON1T likelihood very well.This allows us to input a recoil spectra in energy units and we do not need to apply any additional acceptance, efficiencies or resolution in addition to what is being calculated in the package.For the other experiments that we consider in this work, the efficiencies are readily available in the literature.For our future xenon detector, we use the efficiency curve of LZ from Ref. [69] which already has the acceptance incorporated.For DarkSide-50, the curve from Figure 10 of Ref. [37] is used.For DEAP-3600 we use the information from Figure 21 of Ref. [36].For our future argon sensitivity projections, we use the DarkSide-20k NR acceptance, taken from Figure 92 of Ref. [3].
In order to model resolution effects we take the detector-dependent width σ(E ′ ), with E ′ being the "true" recoil energy.For XENON1T we use the function in Ref. [30], whereas for our future xenon experiment we use the function given in Ref. [70].For DarkSide-50, DEAP-3600, and our future argon detector we use the resolution function given in Ref. [71].The resolutions are given as functions of electron equivalent recoil energy, E ee .In order to translate from nuclear recoil E R , we make use of modified Lindhard factors [72,73].For xenon detectors we use the numerical fit for the parameter k in the Lindhard model determined in Ref. [74], whereas for argon, we make use of the effective Lindhard factor experimentally determined in Ref. [75].We have assumed that the electronic recoil receives negligible amounts of quenching allowing us to interpret E ′ in as the electron equivalent energy E ee .We note that this assumption may not be entirely true and for different experiments, the precise behaviour of ionised electrons and scintillation photons and how their energy is divided may alter how one models the resolution.The method implemented in this work is sufficient for this study and a more sophisticated analysis, like those performed in experimental collaborations will not change the results of this paper.
In order to calculate an exclusion limit or projected sensitivity of an experiment, one needs to determine the expected background signal, N b as well as the observed number of events, N obs .We use the Poissonian probability, where N th = N DM + N b is the sum of the expected background plus dark matter signal given by Eq. ( 9).For both DarkSide-50 and DEAP-3600, no background events were expected and no recoils were observed.
For our projections for DarkSide-20k and LZ, we use binned spectra for the expected background and signal rather than a total counts.The Likelihood from k independent energy bins is simply the product of the Poissonian probability of each bin, and we build a test statistic (TS) from the Likelihood, such that, using the Poisson probability, the TS can be written as where the sum over k is a sum over the bins.This test statistic can be used to determine a two-sided confidence interval.In this case in order to find the 90% exclusion limit or projected sensitivity, the TS is approximated as a χ 2 distribution.
With our future xenon and argon experiment projections being inspired by LZ and DarkSide-20k respectively, we consult the backgrounds reported by the relevant collaboration.For the DarkSide-20k inspired argon detector, the background discrimination between ER and NR events using pulse shape discrimination is expected to be extremely effective, resulting in a near backgroundfree experiment (to a level of < 0.1 events in the entire exposure [3]).After ER backgrounds have been discriminated the only background remaining is the CEνNS background from solar neutrinos [76,77], expected to contribute ∼ 1.6 background events within a 100 tonne-year exposure.We use background and signal spectra in a binned distribution with 1 keV width bins between 20 and 200 keV in terms of the nuclear recoil energy.
In order to attempt to replicate the LZ first search results with an exposure of 60 live days, we use the background spectrum from Figure 6 from Ref. [69] in units of the electron equivalent energy with 3 bins per keV ee between 0 and 17 keV.The dark matter spectra are calculating using the same binning, using Lindhard factor to convert between E R and E ee .To reproduce the effect of the S1-S2 background discrimination employed in the LZ collaboration's first data analysis we make use of Figure 4 from Ref. [69] and scale down the backgrounds to total 11 events (compared to 333 originally in the background distribution).As this region corresponds to the NR 90% quantile region, we also apply a 0.9 flat acceptance to the dark matter signal.This allows us to replicate the expected limit well for the entire mass range, and the observed limit at higher dark matter masses.However, we cannot perfectly match the observed limit in the region of maximum sensitivity around 20 − 30 GeV dark matter mass, which is due to an under-fluctuation of observed events in the S1-S2 search space, due to lack of public information on the exact distribution of expected backgrounds.The official LZ exclusion limit at around 20 − 30 GeV is around three times more stringent than that of our calculation.Our ability to replicate the LZ results is consistent in performance to the updated treatment of the LZ limits as implemented in DDcalc [78,79] by the GAMBIT working group [80].For the XENON1T limit, using the approximate likelihood code gives good agreement with the published limit.
For the argon experiments, we see a very good match with DarkSide-50 and DEAP-3600 for all dark matter candidate masses.For the DarkSide-20k projection, we can match the projection from the 2017 TDR very well by using a 1-bin method with 1.6 neutrino background events per 100 tonne-years of exposure.However, in this paper we use a binned spectrum which gives us an additional 20-50% sensitivity for most masses.Using binned spectra rather than a total number of counts allows us to determine the effect of the spectral shape differences in the isospin-violating models as well as the difference in total rate.Our analysis is able to replicate the isoscalar (c n /c p = 1), isovector (c n /c p = −1) and xenophobic (c n /c p = −0.7)working point limits for five different NREFT operators presented by DEAP-3600 in Ref. [47].
For our projections for the whole 1000 live days exposure of LZ we follow the same method as described above used for the observed limit with 60 live days, but with backgrounds scaled up for the exposure and assume the observed spectra equals the expected backgrounds.
We note that some of the above methods described are applied on data sets with low expected backgrounds.In these cases Wilks' theorem may no longer hold and the test statistic no longer converges to χ 2 , meaning the drawing of limits may be affected [81].We have checked that isospin suppression affects both of the experimental limits calculated with the methods described above and using toyMCs in the same way.Since our simplified limits replicate better the results from the full experimental analyses, we present our results using these.

IV. RESULTS
In this section we show the results for isospin-violating scenarios, using the detector setup described above.First, we present the work for the two gauge invariant cases described in Section II.Second, we present results considering only one NREFT operator at a time.We do this for an incomplete set of operators so as to not overload the main text but results for other operators can be found in Appendix 2.

A. Spin-(in)dependent model
As described in Section II, a gauge invariant formulation of dark matter interactions leads to the requirement of additional operator responses beyond the standard O 1 or O 4 formalism.For the SI model, this is of little consequence because of the large sensitivity difference between O 1 compared to O 7 and O 9 .Nonetheless, our calculations of the leading xenon and argon-based direct detec-tion experiments as well as the future projections serve as a timely update to the literature.This is even more true given the form factor dependence shown in Figure 3.The top row of Figure 4 shows the current and future sensitivities for the SI model.The panels show isospin configurations that are phenomenologically interesting, from left to right we have xenon-favoured, argon-favoured, and isospin-conserving scenarios.We choose the Xe/Arfavoured parameter values by taking the ratio of total events in the LZ exposure and the DarkSide-20k exposure in the ϕ and θ parameter space, and finding the minimum and maximum points.For the spin-independent model this method results in an argon-favoured point of θ/π = 0.62 and ϕ/2π = 0.98, however, there is a range of values with approximately the same relative argon-toxenon suppression ratio that we could have chosen, these follow the c n /c p = −0.72 contour shown in Figure 1.The equivalent maximal xenon-favoured point is found at values of θ/π = 0.45 and ϕ/2π = 0.41.These Ar-favoured and Xe-favoured points have been highlighted on Figure 1 using circle and star markers, respectively.
In the SD case, the situation is more interesting because gauge invariance requires that operator O 4 is accompanied by O 6 , O 8 , and O 9 .Here the different operators provoke nuclear responses that have roughly the same normalisations.From Figure 1, we saw that there should be a high degree of complementarity between argon and xenon target searches.The bottom row of Figure 4 shows the current and future sensitivities for the SD model.Once again we show three panels showing the xenon-favoured, argon-favoured and isospin-conserving scenarios.In the SD model, the relative dependence of the argon and xenon suppression rates with θ and ϕ is more complex due to the contributions of various xenon-sensitive operators.This results in the argon-favoured configuration not coinciding with the region of maximum xenon suppression (which can be seen in the bottom left of Figure 1).This is because where maximal xenon suppression occurs, the argon rate is similarly suppressed.Instead we find the argon favoured region occurs at a value of θ/π = 0.46 and ϕ/2π = 0.07.The xenonfavoured configuration does however align with the region of maximum argon suppression (albeit not the region of minimum xenon suppression) at θ/π = 0.33 and ϕ/2π = 0.72, again demonstrated by the star marker in the bottom panel of Figure 1.For the isospin conserving point we use a value of ϕ = π/4, for which the proton coupling equals the neutron coupling for all operators involved, for both models.At this ϕ value, isospin conservation holds for any value of θ and we show the limits for θ = 0.
Figure 4 shows the experimental exclusion limits which were calculated using existing public data from DarkSide-50, DEAP-3600, XENON1T, and LZ (with 60 days livetime), and derive projected sensitivities for DarkSide-20k and LZ (1,000 days livetime).For both the SI and SD models, results are presented as an exclusion limit at the 90% confidence level on (g ′ /Λ) 4 , as this quantity is proportional to the dark matter event rate and therefore scales with exposure in the same way as familiar zeromomentum limit cross-section constraints.
We note that the effects of isospin violation can have dramatic impacts on the apparent parameter space of dark matter currently excluded and expected to be probed in the future.For the SI case, Figure 4 shows a suppression at m χ ∼ 10 3 GeV of ∼ 10 6 and ∼ 10 5 for argon and xenon respectively.We see that currently, due to relatively small target exposures, argon limits are unable to lead in sensitivity even in the argon-favoured scenario, but this situation will change dramatically in the future.Future detectors will be of similar sensitivity, where at high m χ , the isospin-symmetric ∼ 10 difference is simply in line with the expectations due to the relative exposures.We see that this difference can be increased or decreased by ∼ 10 2 or ∼ 10 3 depending on the parameters chosen.Results of this kind have been widely reported in the literature before, but we feel having a dedicated up-to-date result will be instructive for the direct detection community.Furthermore, by including the effects of operators O 7 and O 9 we have been able to verify that even in manifestly gauge invariant models, the suppression scales seen by just considering O 1 , are not greatly affected.
Of greater novelty are the consequences we report for argon detectors in the gauge invariant SD model.The bottom row of Figure 4 shows how argon detectors are currently sensitive to SD dark matter models and will be able to beat xenon detectors in the future.We do see across the three panels, that the xenon sensitivities are fairly stable, whereas for argon the sensitivity in the high m χ region varies substantially, by a factor of ∼ 10 5 .This is thanks to the additional operators required by gauge invariance which substantially limits the points in parameter space where xenon suppression maximally occurs, see Fig. 1.Additionally, the regions where xenon is most suppressed are not aligned with areas where argon is least suppressed.If instead of showing the xenon and argon favoured points in Fig. 5, we showed minimally and maximally suppressed xenon (argon) we would see a factor ∼ 10 3 (∼ 10 6 ) difference in the respective limits.
We see that even in the isospin-conserving scenario, the future argon detector is more sensitive at high masses.This is quite the coup for a target that is insensitive to the canonical SD nuclear response, O 4 , and is simply a consequence of insisting on SM gauge invariance of the dark matter model under consideration.As shown in Ref. [50], SM gauge invariance requires that the O 8 operator is nonzero, provoking a nuclear response from the argon atom.Since our Majorana fermion model is simply motivated by the dark axial vector current χγ µ γ 5 χ, the standard source of spin-dependent interactions, we believe this result to be fairly general.
In this paper we only explicitly consider data and experimental analyses from LXe/LAr TPCs using both the S1 and S2 ionisation/scintillation signals targeting > ∼ 10 GeV mass dark matter candidates.However, we note that the XENON1T and DarkSide-50 experiments have established [38,39,82,83] that an ionisation-only analysis accounting for the impact of the Migdal effect and bremsstrahlung on electronic recoils allows for probes of dark matter to masses as low as 85 MeV and 40 MeV, respectively.Considering the published sensitivities of XENON1T and DarkSide-50 for dark matter masses below 10 GeV and the suppression results for the next generation experiments derived in this paper, we estimate that the relative dominance of future xenon and argon experiments illustrated in Figure 4 at intermediate-tohigh masses would extend down to the O(10) MeV-scale if data from ionisation-only analyses were also considered.An exception to this that we highlight is that it may be possible for DarkSide-20k to ultimately be more sensitive that LZ in the low mass region in the SI model even in the Xe-favoured regions of parameter space.

B. NREFT results
In our view, the two models considered above are wellmotivated and allow us to exhibit the main points of this paper.We are aware that there are numerous other possible candidates for dark matter and many of which are as theoretically attractive.In lieu of producing a comprehensive study of all the possible isospin-violating models, we do the next best thing and perform a study on NREFT operators in isolation.Those with their own dark matter models will be able to use the results presented here to estimate the implications isospin violation may have on the expected direct detection signal in current and future experiments.
The way we present our results is quite different from the previous section.This is because we want to maximise the utility for interested parties and present results in a way that comprehensively covers the parameter space.Therefore, we want to convey how expected counts is affected over a range of c n /c p .Additionally, by performing our analysis in a more realistic way, i.e. not simply observing suppressions in the E R → 0 limit, we are now sensitive to variations of the signal shape coming from the nuclear response functions, as outlined in Section II.This introduces a mass dependence on the suppression predicted because lower values of m χ have their spectral shape dominated by the tail of the dark matter velocity distribution, and the light incident dark matter is unable to impart enough momentum transfer to probe the structure function.For larger m χ , the signal is shallower, meaning that the isospin cancellation has to persist over a greater energy region.We quantify the impact of a change in coupling ratio on the dark matter -nucleon coupling rate (or Wilson coefficient) constraints via a suppression factor, S, which is determined as a ratio of the extracted limit at a given c n /c p to the limit at the extracted 90% C.L. limit evaluated in the isospin-conserving case: This suppression factor is determined by explicitly rederiving the experimental constraints with a new dark matter signal model, accounting for the energy recoil thresholds, efficiencies, and background contributions of the experiments under consideration as described in Sec.III.In Figure 5, the suppression factor for the DarkSide-20k-like argon limit and the LZ-like xenon experimental limit is plotted for neutron-to-proton coupling ratios between −1 and +1 for momentum-independent operators O 1 and O 8 and for the momentum-dependent operator O 11 .We present the inverse of the S suppression factor, meaning lower values correspond to a stronger suppression of the limits.The suppression behaviour is strongly dependent on the value of the neutron-to-proton coupling ratio, with target and operator-dependent values where experimental sensitivity to dark matter would be maximally suppressed.For the momentum-independent O 1 operator in a xenon target, this is around c n /c p = −0.7,compared to an argon target where this peak occurs around c n /c p = −0.8.The dependence of the suppression factor on the coupling ratio near the region of peak suppression is also affected by the form factors used in the calculation, as could be inferred by Figure 3, we explicitly show the effects in Figure 8 in the Appendix for the O 1 operator.
For O 1 , the dark matter mass dependence of the suppression is minimal.Still, results for the O 8 and O 11 operators illustrate significant mass dependence of the coupling ratio value corresponding to maximum suppression as well as the amount of suppression.This is particularly true for xenon and much less so for argon.Given this mass-dependent effect, we also present in the top row of Figure 6 a two-dimensional version of these plots (here only O 8 , for illustration), again redetermining the projected experimental constraints for each experiment at each mass and coupling ratio point, to illustrate how the limit suppression relative to the isospin-symmetric scenario varies as a function of dark matter candidate mass and coupling ratio.Again we present the inverse of the suppression factor in the top row 1/S, where large negative values (lighter colours) represent a stronger suppression.We also display the ratio of the suppression in argon (for a DarkSide-20k-like detector) and xenon (for an LZ-like detector), shown in the bottom left pane of Figure 6.This metric provides information about how the experimental limits for different targets will be affected in comparison to each other.For these plots, green regions denote points where the xenon LZ-like limit is suppressed more strongly, and in the purple areas, the argon DarkSide-20k-like experimental limit is suppressed to a greater degree.Finally, we also present in the bottom right panel the ratio of the exclusion limits, to demonstrate the mass and coupling ratio combinations for which DarkSide-20k or LZ, with their full exposures, would be projected to be the most sensitive.Rather than being a measure of the relative suppressions in Ar compared to Xe regardless of the nominal sensitivity, this gives a sense of which target is more sensitive for the specific exposures and setup of DarkSide-20k and LZ.We highlight that incorporating the efficiency (and the model used) affects the signal spectra and projected experimental limits, therefore the maximal suppression rate and c n /c p value are also modified.We have explored the effect of efficiencies by switching on and off certain experimental inputs and find that the effect on the suppression is negligible and thus the results presented here are insensitive to changes in the expected performance of these detectors.
The equivalent to Figure 6 for four other operators (O 1 , O 3 , O 5 , and O 11 ) can be found in the Appendix (Figure 9).We note that the result for O 1 is fairly massindependent in terms of the c n /c p value corresponding to maximum suppression, although the maximum amount of suppression reduces with increasing dark matter mass.The equivalent of the bottom right plot of Figure 6 for the other four operators, where we present the ratio of the exclusion limits, is presented in Figure 11  Additionally, we evaluate the experimental limits for the operators separately for the current and future detectors as we did for the models in the previous section.In Figure 7 we present the experimental limits for operator O 8 evaluated as an exclusion limit on the operator coefficient c 8 .The xenon-favoured and argonfavoured c n /c p values were chosen as the coupling ratio that maximised/minimised the suppression ratio shown in Figure 6 at a benchmark mass of 100 GeV, roughly the mass at which the experiments are most sensitive.The equivalent to Figure 7 for other operators can be found in the Appendix (Figures 12 and 13).
Our aim here is that the results presented will be useful for those who want to estimate the effects of specific isospin-violating models.For example, say one had a model that predicted a O 8 nuclear response with a neutron-to-proton coupling factor of c n /c p = −0.5, one could take our projected sensitivities for the isospinsymmetric case on the right-most panel of Figure 7, and use the top two panels of Figure 6 to determine the value of S. For c n /c p = −0.5, and the O 8 interaction as an example, a 100 GeV dark matter candidate has suppression factors of S ≈ 35 and ≈ 30 for xenon and argon respectively.Using the respective limits and the value for S as in Eq. ( 13), one can determine how future LZ and DarkSide-20k experiments will constrain the model in question.The isospin-symmetric limits at 100 GeV are c 2 8 ⟨v⟩ 4 Higgs = 8.0 × 10 −4 and 2.0 × 10 −4 for LZ and DarkSide-20k in the isospin symmetric case, so for the spin-dependent O 8 operator and c n /c p = −0.5 case, the limits would be 2.8 × 10 −2 and 6.0 × 10 −3 for LZ and DarkSide-20k, respectively.Given the insensitivity of our calculated suppression ratios to experimental features such as efficiency effects and exposure, these values could also be used to estimate the sensitivity of future xenon and argon experiments to isospin-violating dark matter.

V. CONCLUSION
In this article, we have revisited the case of isospinviolation in dark matter interactions with the SM.This occurs in a substantial number of models and has the potential to significantly suppress the recoil rate in direct detection and thus weaken the apparent current and future projected constraints by many orders of magnitude.We pay particular attention to the crucial complementarity of xenon and argon-target noble liquid detectors which, in the coming decade, will be the most sensitive for m χ ≥ 10 GeV for S1+S2 analyses and m χ ≥ 50 MeV for S2-only analyses.In this work, we have focused on the S1+S2 analysis strategy, but note that many of the qualitative results would extend down to 50 MeV when an ionisation analysis is incorporated.Additionally we dedicated some effort to highlight the importance of performing a spectral analysis for isospin-violating models.We showed that the magnitude of suppression can be altered by the dark matter mass, experimental configuration and nuclear form factors.We have tried to account for variations in the former two throughout this work, and have simply exhibited the importance of the latter, hoping to motivate more dedicated work.
To explore these points in a concrete way we have primarily focused on two gauge-invariant effective models of dark matter, which are inspired by Dirac and Majorana dark matter candidates that couple to the SM via vector currents.These two examples we refer to as SI and SD models, due to their dominant operator responses.We have performed new analyses of the best available public data (from the XENON1T, LZ, DarkSide-50, and DEAP-3600 experiments) to set new constraints on isospinviolating dark matter and determined the effects in the case where one target element is particularly favoured.We do the same for future sensitivity projections in LZlike and DarkSide-20k-like experimental setups.Interestingly we find that, due to gauge invariance, the SD model, which is canonically considered to be invisible to argon target experiments, is not.In fact, even in the isospin-conserving case, our projections show DarkSide-20k will be more sensitive than xenon for SD interactions for m χ > 100 GeV.Furthermore, we have performed the same analysis on the more model-independent and comprehensive NREFT operators developed for direct detection, in the hope that those in the model-building community can apply our projections to their isospin-violating models.In Section IV we only presented one operator fully, O 8 , in order to describe our presentation and how one may use our results.We have provided our results for operators O 1−11 in the Appendix.

DATA AVAILABILITY STATEMENT
Data available from the main results of this paper is available at Ref. [84].DarkSide-20k, , F Fitz. FIG. 8. Suppression of the experimental constraints on the dark matter -nucleon interaction strength with neutron-to-proton coupling ratio compared to the constraints for the isospin-conserving scenario at a representative dark matter mass of 20 GeV, for the O1 interaction for the two form factors considered in this paper.Efficiency effects and energy window constraints for the LZ and DarkSide-20k detectors have been incorporated.

Mass and coupling dependence of suppression of experimental constraints
Figure 9 and 10 show two dimensional distributions illustrating how experimental constraints on dark matternucleon interaction strength are affected by both m χ and c n /c p .We show this dependence for multiple NREFT operators that are not included in the main body of the text (Figure 9) and for the operators only present for xenon targets and not for argon (Figure 10).Additionally, we show the ratio of the projected exclusion limit in DarkSide-20k and LZ in the same parameter space, for four operators in Figure 11.

FIG. 1 .
FIG.1.Relative dR/dER for xenon (left) and argon (right) following the parametrisation in Eq. (7) for 100 GeV dark matter mass.The top row reflects results for the spin-independent vector-like couplings model (where the contribution of the O1, O7, and O9 operators is considered) while the bottom row represents the situation in the spin-dependent axial vector-like couplings model (with operators O4, O6, O8, and O9).In the spin-independent case the contour of maximal suppression corresponding to a fixed neutron-to-proton coupling ratio is highlighted by a dashed line.The star and circle markers correspond to specific model points studied and detailed in Section IV.

− 1 Xe, O 1 − 1 Ar, O 1 9 FIG. 3 .
FIG.3.Difference in the recoil spectra in xenon (left) and argon (right) for dark matter mass of 100 GeV between the two sets of form factors (dashed versus solid lines), for three different coupling ratios.Note we have set c p 1 ⟨v⟩Higgs = 10 −5 and c n 1 = cn/cp × c p 1 where cn/cp is the isospin ratio as signified in the caption.

11 FIG. 5 .
FIG.5.Inverse of the suppression factor for the experimental constraints on the dark matter -nucleon interaction strength with neutron-to-proton coupling ratio, for a LZ-like and DarkSide-20k-like experiment.Shown here is a comparison of the limit suppression factors for EFT operators O1, O8, and O11 for dark matter masses 20 GeV and 2 TeV.

FIG. 6 .
FIG. 6. Top row: Inverse the suppression factor for an LZ-like experiment (top left) and a DarkSide-20k-like experiment (top right) as a function of neutron-to-proton coupling ratio and dark matter candidate mass, for the O8 operator.Large negative values indicate a strong suppression of the experimental limit.Bottom left: Ratio of suppression in LZ-like to DarkSide-20k-like experiment.Negative z-values (green regions) identify regions where LZ would suffer from less suppression than Darkside-20k, with positive values (purple regions) identifying the opposite.Bottom right: : Ratio of the projected 90% exclusion limit in DarkSide-20k to LZ, where here negative values (green regions) indicate where the full exposure of an LZ-like detector would be the more sensitive experiment for a particular mass and coupling ratio while positive values (purple) indicates where DarkSide-20k would be the most sensitive.
of the Appendix.For completeness, two-dimensional suppression factors for operators for which xenon is sensitive, but argon is not due to its nuclear spin being zero (O 4 , O 6 , O 7 , O 9 , and O 10 ) are also presented in the Appendix (Figure 10) and are shown to be relatively mass independent with maximum suppression at c n /c p close to zero.For O 4 the maximum suppression occurs at c n /c p = −0.03.

FIG. 7 .
FIG. 7. Current LXe and dark matter exclusion regions [shaded] and future projections for two near-future liquid noble direct detection experiment sensitivities (a 200 tonne-year argon experiment and a 15 tonne-year xenon experiment) [dashed] as a function of nucleon -dark matter EFT interaction coefficient c8 for the O8 operator a xenon-favoured model scenario (left), an argon-favoured model scenario (middle), and for the standard isospin-conserving case (right).

APPENDIX 1 .
Limit suppression for the O1 interaction with neutron-to-proton coupling ratio and form factor variations

Figure 8
Figure8illustrates how the experimental constraints on dark matter -nucleon interaction strength are suppressed as a function of the coupling ratio c n /c p to the O 1 interaction for different targets at a dark matter candidate mass of 20 GeV, and how these are substantially affected by the choice of form factors F Hof. and F Fitz. .

10 (
FIG.9.Inverse suppression factor of experimental constraints on dark matter -nucleon interaction strength as a function of neutron-to-proton coupling ratio and dark matter mass for an LZ-like experiment (left), DarkSide-20k-like experiment (middle) and the ratio of DarkSide-20k-like to LZ-like (right) for the O1, O3, O5, and O11 operators.

FIG. 11 .
FIG.11.Ratio of the projected experimental constraints on dark matter -nucleon interaction strength as a function of neutron-to-proton coupling ratio and dark matter mass for an LZ-like experiment versus a DarkSide-20k-like experiment with their full projected exposures for the O1, O3, O5, and O11 operators.

TABLE I
[28]st of the NREFT operators constructed to obey Galilean invariance, introduced by[26].More operators have been introduced recently[28], but we list only operators we consider in this work for brevity.