Dark biportals at direct detection

We present a study of the constraining power of low-threshold, high-resolution direct-detection experiments for a biportal simplified dark matter model. In the scenario we consider here, dark matter and Standard Model particles interact through two dark vector mediators $-$ one light and one heavy with respect to the momentum transfer in the experiment. Interference effects at the level of scattering amplitudes can lead to novel marked shape features in the differential recoil spectra, which are best exploited by high-resolution, low-threshold experiments. We identify the region in parameters space for our model where such effects are dominant and show that composite-target experiments with large atomic mass differences are ideal to explore these scenarios. We develop a profile likelihood approach to analyze presently available and future data. Using published results by the CRESST-III experiment and projections of future sensitivities for the COSINUS experiment, we constrain the parameter space in our model, thereby showing the potential of such an analysis on a class of dark matter models which exhibit non-standard features in the recoil spectra.


Introduction
An important avenue in the search for particle dark matter (DM) is direct detection, which aims to reveal interactions between DM and Standard Model (SM) fields in ultra-sensitive, Earth-bound detectors. Traditionally, these experiments have focused on measuring elastic scattering of nonrelativistic DM particles off atomic nuclei, with the dominant background typically arising from electron recoils. The key observable to identify a possible DM signal in these analyses is the differential nuclear recoil spectrum, whose shape is dictated by underlying particle physics models.
The interactions underlying scattering processes involved in these direct DM searches can be well approximated by effective contact operators if the mediator fields are sufficiently heavy compared to the typical momentum transfer [1][2][3][4], which is of the order of 1-100 keV. If this is not the case, the expected changes in the shape of the recoil spectrum make low-threshold experiments respond better to the interaction mechanism [5] (see also Refs. [6,7]).
In the last years, enormous technological progress has been made in lowering the nuclear recoil threshold of direct DM detection. A recent compilation of current results may be found, for instance, in Ref. [8]. In order to better understand the reach of these experiments, it is thus interesting and timely to analyze DM models that predict nonstandard features of the shape of the recoil spectrum which these low-threshold experiments might be sensitive to. Such analyses not only help demonstrate the importance of these experiments but may also point to previously unexplored signatures, which could be investigated with improved detector design. Motivated by these considerations, in this paper we consider the case study of a biportal simplified model of DM-quark interactions and focus on two cryogenic experiments with composite targets, namely, CRESST [9][10][11][12], whose detector material is calcium tungstate (CaWO 4 ) and which features a low threshold, and the upcoming COSINUS experiment [13], which uses sodium iodide (NaI). The choice of COS-INUS is motivated by the large mass difference between the target elements, which (as we demonstrate later) is important for probing our scenarios.
The very concept of DM-quark interaction via multiple mediators is not new and has been previously studied in the literature; see, e.g., Ref. [14] and references therein. For example, in supersymmetric gauge theories, DM-nucleon scattering can receive contributions from both Higgs and squark exchanges. Models with both a scalar and a Z portal to the dark sector naturally arise in extensions of the SM gauge group involving a new U (1) [15,16]. Models with two vector mediators were studied in Ref. [17] in the context of DM indirect detection and in Ref. [18]. For the case of two scalar mediators, interference effects on LHC searches were studied in Ref. [19], while direct-detection and associated flavor-physics constraints were discussed in Refs. [20,21]. A Higgs-leptoquark portal to dark matter was the subject of Ref. [22]. Two-scalar-mediator models were also explored in experimental searches at the LHC, as in Ref. [23]. The phenomenology of the two-Higgsdoublet plus pseudoscalar model was studied in Refs. [24,25], and more general scenarios of multi-pseudoscalar-mediated DM were recently analyzed in Ref. [26], with the goal of overcoming inherent limitations of the approach based on simplified models with a single mediator. Multiple DM-quark interactions have also been studied in the framework of an effective field theory (EFT) taking into account interference effects [27].
Our focus here is on the direct-detection phenomenology of a simplified model with two dark vector mediators, V l and V h , coupling to both DM and SM quarks. Crucially, a large hierarchy between the masses of these mediators is assumed, namely, m 2 V l ≤ q 2 and m 2 V h q 2 , where q 2 denotes the squared momentum transfer. The resulting interference effects at the level of scattering amplitudes lead to interesting novel features in this model like sharp dips in the recoil spectra. Our analysis demonstrates the importance of lowthreshold, high-resolution experiments in possibly detecting these features. Novel aspects of our study include the application of a profile likelihood approach to determine exclusion limits on the model parameters. For this purpose, we make use of the public CRESST-III data release, as well as the future sensitivity of the COSINUS experiment and study the target dependence of the recoil spectra in our model setup.
The paper is organized as follows. In Sec. 2 we introduce our vector biportal model, derive the expression for the differential recoil spectrum from the tree-level DM-nucleus cross section, discuss its characteristic features due to constructive and destructive interference and assess the effects of detector-specific quantities such as threshold and resolution for benchmark values of the model parameters. In Sec. 3 we discuss the likelihood formalism we applied to use published real data from the CRESST-III experiment and the projected sensitivity of the COSINUS experiment. Exclusions limits for the model parameters are presented and discussed in Sec. 4 and Sec. 6 contains our conclusions.

Vector biportal model at direct-detection experiments
The simplified model setup we consider consists of a Dirac DM particle χ interacting in a flavor-independent way with the SM quarks via a light and a heavy vector mediator (V l and V h , respectively), with different couplings. The corresponding interaction terms in the relativistic Lagrangian amount to where the SM quarks are denoted by Q in order to avoid confusion with the momentum transfer q. The dimensionless constants g Ql , g χl denote the couplings between V l and the quarks, and V l and the DM particle, respectively. Analogous definitions hold for g Qh , g χh in the case of V h . The Lagrangian above gives a general, effective, low-energy description of the interactions relevant for our analysis and is representative of a whole class of UV-complete theories featuring two vector mediators with a large mass hierarchy and flavor-independent vector couplings. This theoretical scenario, as we are going to show, is characterized by the fact that the resulting nuclear recoil spectra show nonstandard features like marked dips, due to interference effects. We point out that in the future it will also be interesting to investigate the modifications to these features of the recoil spectra corresponding to a generalization of our biportal simplified model to include further coupling structures as well as flavor-dependent mediator couplings.
We first derive the differential DM-nucleus cross section in terms of couplings and masses in our model. This leads to the expression for the recoil rate, which our numerical analysis is based on. The total unpolarized squared scattering amplitude from the diagrams shown in Fig. 1 is given in the nonrelativistic limit by . (2.2) In Fig. 1, the contact interaction is generated by integrating out the heavy mediator field V h . The magnitude of the momentum transfer amounts to q ≡ q 2 = √ 2m N E R . Here A denotes the mass number of the target nucleus and m N is its mass. The structure of the nucleus at nonzero momentum transfer is approximated by the standard form factor F 2 (q). 1 In our numerical analysis, we use the functional form for F (q) first proposed by Helm [30], with the nuclear radius parametrization by Lewin and Smith [31].
An essential feature of the model is the interference term arising in the squared scattering amplitude. This plays an important role in the resulting phenomenology, as we are going to discuss in detail. It is also noteworthy that such an interference term would also be present in the case of the analogous exchange of two scalar mediators [19][20][21]32].
For the light mediator we define an overall dimensionless coupling g l ≡ g χl g Ql since the tree-level squared amplitude depends on this product and not on the separate values of g χl and g Ql . Also for the heavy mediator, since we are only sensitive to the ratio between the product of the couplings and the mediator mass squared, we introduce an effective The dependencies on the couplings are relevant when thinking about the complementarity of the constraints on the biportal models from collider experiments and relic density mechanisms.
From Eq. (2.2) the tree-level differential cross section for nonrelativistic kinematics is given by where the third term in brackets is the contribution from interference. For g l g h > 0 we have constructive interference, while for g l g h < 0 there is destructive interference. For a given target, the free parameters in this equation are g l , g h and m V l . The observable of interest in direct-detection experiments is the differential nuclear recoil rate. As we will see later on, experiments with high resolution and corresponding low threshold are needed in order to investigate the full phenomenology of the differential recoil spectrum in our model. This is given as a function of the differential cross section above by where ρ 0 is the local DM density, N is the (fiducial) mass of the detector, | v| is the DM speed in the lab frame, v min (E R ) is the minimum DM speed required to generate a nuclear recoil with energy E R for an elastic collision and v esc is the Galactic escape speed. The DM velocity distribution, denoted by f ( v), is in principle modulated in time due to the Earth's motion around the Sun, but we average out this effect here as we consider measurement times of 1 year or longer. We assume a normalized isotropic Maxwell-Boltzmann distribution in the Galactic rest frame 2 v ∞ is the root-mean-square velocity, which is related to the asymptotic value of the rotational velocity v ∞ , and z 2 = 3v 2 esc /(2σ 2 v ) [33]. The distribution f ( v) in the Earth's rest frame is obtained by applying a Galilean transformation. In our numerical analysis we take ρ 0 = 0.3 GeV c −2 cm −3 [34], v ∞ = 220 km s −1 , and set the Earth's average speed equal to 1.05 v ∞ [33]. We also take v esc = 544 km s −1 , which is in agreement with the bounds set in Refs. [35,36]. In the following we will neglect uncertainties in the determination of the local DM density (see, e.g., Ref. [37] for a review) as well as in the DM velocity distribution since a treatment of these effects is outside the scope of this work. Finally, for detectors with composite targets, a weighted sum over all targets in Eq. (2.4) needs to be performed.

Features of the recoil spectra in the biportal model
We now illustrate and discuss the peculiar features of the differential recoil spectra in our biportal model as compared to the single-mediator case. For this purpose, in this section we fix the four parameters that are left free in Eq. (2.4 to the benchmark values listed in Table 1 and explore the effects of the variations of the parameters around these values. In the benchmark model the mass of the DM particle and the light mediator V l are set to 5 GeV and 10 MeV, respectively. The values of g l and the effective coupling g h are chosen such that differential recoil rates are compatible with current limits on the spin-independent DM-nucleon cross section [38]. For a DM mass of 5 GeV with σ SI < 10 −42 cm 2 [38], we estimate g l = O(10 −12 ) and g h = O(10 −14 MeV −2 ) for an m V h of about 100 MeV. The benchmark values of the couplings for the light and heavy mediator are chosen to satisfy this constraint while aiming to maximize interference effects.
We briefly comment on the relevance of collider constraints on our model parameters before proceeding further with our discussion. 2   Since what is relevant in our analysis is the product of DM-V l and SM-V l couplings, we can probe the parameter space where DM-V l is larger and thus can escape the collider constraints. Similar arguments apply to heavier mediators as well. The relevant searches for heavier mediators arise from fixed-target experiments as well as colliders. The relevant constraints depend on the mass of the mediator as well as the mediator branching ratio to specific final states. For example, if the heavy mediators have masses of the order of 100 GeV, LHC constraints via resonance searches are applicable. In this case the mediator couplings to leptons are more strongly constrained than mediator couplings to quarks or DM. For mediators lighter than 100 GeV, other experimental constraints apply and depending on the mediator branching ratio, the exact limits differ. However, at direct detection experiments mediators with masses larger than about 10 MeV are considered heavy 3 and hence the exact constraints become model dependent.
We now discuss the features of the differential nuclear recoil spectra in the benchmark vector biportal model. The cases of a single-element target, namely germanium (Ge), and the composite target sodium iodide (NaI) are shown in Figs. 2a and 2b, respectively, for both constructive (solid turquoise line) and destructive interference (dashed turquoise line).
To compare the biportal model with single-mediator models, we also show the spectra for the heavy (blue line) and light mediator component (grey line).
As expected, the constructive biportal model predicts an overall higher differential recoil rate, irrespective of the target. The shape is a nonlinear combination of the two individual single-mediator rates due to the presence of a non-negligible interference term. In the case of destructive interference, the situation becomes even more interesting.
As shown in Fig. 2a, an overall lower rate is observed, but more interestingly a strong interference effect is present at around E R = 0.75 keV. This dip appears precisely where the rates for the single heavy and single light mediators overlap and, due to the negative interference term, lead to a strong suppression. We stress that this dip is unique to the biportal model, as similar spectral features are not present in the case of interactions via one mediator. While the presence of this kink is a universal feature of biportal models, the exact location and depth depend on the target material as well as all of the free model parameters.
The impact of the interference term becomes even more apparent for composite targets. As shown in Fig. 2b, the biportal model with constructive interference (solid turquoise line) inherits a damped version of the kink present for the single heavy mediator case around a recoil energy of 2 keV. For NaI, this kink in the single heavy mediator spectrum marks the change of the dominant target which is iodine at low recoil energies and sodium at higher recoil energies. The recoil spectrum for the biportal model with destructive interference (dashed turquoise line) in Fig. 2b now shows two interference effects. The first feature appears at rather low recoil energies, which is the iodine-dominated regime of the spectrum, and the other feature is present in the sodium-dominated regime. In the case of destructive interference, the biportal model can thus lead to multiple pronounced features in the expected recoil spectra of composite targets. However, these features are maximally enhanced when the various components of the target material vary strongly in their atomic mass/weight (e.g., Al 2 O 3 only exhibits one spectral feature despite being a multielement target 4 ). In order to probe the biportal model, especially the destructive case, composite targets with a pronounced mass difference between their elements are thus advantageous.
In the next step, we investigate the impact of the choice of couplings in the biportal  Figure 3: Differential event rates for the biportal mediator model on a NaI target for the constructive (a) and destructive cases (b). The light mediator mass is fixed at 10 MeV, the DM mass is fixed at 5 GeV, and the effective coupling of the heavy mediator is fixed at 1 × 10 −14 MeV −2 . Each line shows the spectrum for a different value of the light mediator coupling g l . Neither the detector threshold nor the energy resolution are considered.
model. In order to do this in a systematic way, the effective heavy mediator coupling is fixed while we vary g l . This analysis allows us to study the limiting cases of coupling hierarchies corresponding to either single heavy or single light mediator scenarios and, more importantly, to identify interesting ranges of coupling hierarchies which lead to distinct biportal model phenomenology. In Fig. 3, the left panel shows the constructive case on the NaI target. For the highest value of the light coupling (light turquoise line), the light mediator dominates the interaction and the spectrum thus looks like that for a single light mediator. On the other hand, for the smallest value of g l the heavy mediator is the dominant portal and the spectrum thus resembles the single heavy mediator case, including the characteristic kink for NaI at a recoil energy of around 2 keV. In the destructive case shown in the right panel of Fig. 3 one can observe the same limiting cases of the spectra for the highest and lowest absolute values of g l . However, for all other values the interference characteristics are visible. In order for a dip to appear in the biportal spectrum, the individual spectra of the single light and single heavy mediator components need to overlap or lie very close to each other, corresponding to complete or almost complete effacement of the total rate. In the case of smaller contributions to the interaction from the light mediator, such an overlap or approach can only take place at low recoil energies, where the light mediator component shows an enhanced rate. For lower absolute values of g l , dips thus only appear in the iodine-dominated part of the recoil spectrum close to the threshold. If g l is increased, dips also start to show in the sodium-dominated area (dark blue line). For higher absolute values of g l , a possible overlap of rates is only achievable at higher recoil energies outside the sodium-dominated regime. This dip gets more pronounced with further increased coupling of the light mediator and moves towards higher energies. At some point, the contribution to the differential recoil rate from the light mediator is much larger than that from the heavy mediator for all recoil energies. In addition to differential recoil spectra, the total number of expected events also offers a tool to compare the phenomenological predictions of various DM-SM interaction models with the experimental data. The total rate per unit exposure is given by where E thr is the energy threshold of the experiment and E max is the maximum energy of the recoil spectrum. In the biportal model, understanding the total event rate as a function of the light mediator coupling for a fixed value of g h is especially instructive. Setting for the moment E thr = 0, the predictions for the total rate R(|g l |) are shown in Fig. 4 for the biportal model (turquoise; constructive as solid, destructive as dashed lines), as well as for a single light (grey) and a single heavy (blue) mediator model. Showing the total number of events as a function of g l results in a constant rate for the single heavy mediator case and a rate scaling linearly with the squared coupling for the single light mediator case, respectively. As previously seen for the differential event rate, for low values of g l the total event rate predicted by the biportal model coincides with the single heavy mediator case.
On the other hand, for high values of g l the light mediator dominates the interaction. While the constructive case is characterized by a smooth transition between the two limiting cases in Fig. 4, the destructive case predicts a significantly reduced number of total events for certain values of g l . This suppression takes place for precisely those values of the light mediator coupling for which the destructive interference effects are enhanced for the given choice of g h . Moreover, in the destructive case the total event rate as a function of the light mediator coupling R(|g l |) is not injective. In other words, a certain number of events can be expected for different combinations of the two couplings. This is an important feature of the destructive biportal model and will play a role later in the discussion of exclusion limits.

Impact of detector-specific quantities
Detector-specific parameters, such as threshold or energy resolution, were not considered in the previous subsection. However, these are important features of direct-detection experiments and have to be taken into account whenever a statement is made about the sensitivity of DM searches to the model parameter space. In the following we investigate the impact of the detector energy resolution on differential recoil spectra off NaI, as well as on total event rates. We model the finite energy resolution of a detector by a Gaussian distribution of width σ. The differential event rate is then convolved with this Gaussian distribution.
In Fig. 5 we demonstrate the impact of this convolution with the resolution on the shape of the differential event rate for the destructive interference case in the biportal model on a NaI target. The masses are set equal to their benchmark values, the effective coupling of the heavy mediator to 10 −14 MeV −2 and g l = −4 × 10 −12 . The various lines correspond to different values of the resolution σ, while the turquoise dashed line shows the spectrum without finite resolution. Poor detector resolutions, i.e. high values of σ, can wash out the destructive interference features significantly. However, for high-resolution  Event rates are shown for single light mediator (grey), single heavy mediator (blue), and biportal models (turquoise) for both constructive (solid lines) and destructive interference (dashed lines). The DM particle mass is fixed at 5 GeV, the light mediator mass is fixed at 10 MeV, and the effective coupling for the heavy mediator is set to 1 × 10 −14 MeV −2 . Neither the detector threshold nor the energy resolution effect are considered.
experiments with values of σ 0.05 keV, such as CRESST-III, 5 measurement of the interference dips is in principle possible. Moreover, one has to consider that resolution and threshold are not two independent quantities, but in fact are closely related. The threshold of an experiment is usually around 4-6 times the value of the resolution σ [40]. In addition to the flattening of peaks, the corresponding threshold thus might cause an interference feature to lie completely out of the sensitive region of a detector. While this might be the case for the first feature in the spectrum shown in Fig. 5 for the COSINUS experiment with a projected resolution of 0.2 keV, the second one is still within reach. A change of couplings might also reposition the first destructive feature in an energy region accessible to COSINUS.
The impact of the detector resolution on the expected total event rate as a function of g l is shown in Fig. 6 for both constructive (solid lines) and destructive (dashed lines) scenarios. Three different cases of detector resolution corresponding to the ideal case (turquoise), 0.02 keV (dark blue), and 0.2 keV (light blue) are considered when computing the total event rate. Our results show that the constructive case is less affected by the detector resolution than the destructive case. For a poor detector resolution, the total event rate for the destructive biportal model can be significantly reduced. As mentioned above, there is a close connection between threshold and resolution. Therefore, as the resolution worsens, the threshold increases, effectively discarding the exponential rise at low recoil energies. Accounting for both of these effects in Fig. 6 leads to a further reduction of the overall rate for poor detector resolution. Although the integration is carried out from the threshold over the whole available spectrum, one should note that due to the overall exponentially falling spectrum, the number of total events is dominated by the shape of the spectrum within the first few keV of recoil energy. A higher threshold can thus significantly lower the number of total events. The vertical distance between the various lines in Fig. 6 grows slightly with increasing |g l | as the impact of the threshold is stronger for interactions that are dominated by light mediators. This can be explained by the fact that the differential recoil rates for light mediator interactions are enhanced at low recoil energies, and those low-energy events are then cut off by a high threshold. Before proceeding further, a word about the shape of the event rate is in order. The interference contributions play a significant role only away from the limiting cases where either the light or the heavy mediator are dominant. While constructive interference generally increases the rate, in the case of destructive interference, it is possible to obtain a drastic reduction in the total rate. Such a reduction will play an important role in the upcoming discussion.

Profile likelihood analysis setup
In the next step we do not set the free parameters of the biportal model equal to their benchmark values but instead constrain them using results of current experiments and projections for future ones. As pointed out in the previous section, low-threshold experiments with composite targets are are especially sensitive to the features of the the biportal model with mass hierarchy. We thus perform our analysis using data published by the CRESST Collaboration for their best-performing CaWO 4 detector in the first run of the CRESST-III experiment [9,41], as well as for the future COSINUS experiment [13]. The latter is a logical choice due to its target material NaI, which has proven especially suitable for our study due to the large difference in mass numbers for sodium and iodine. The experiments we consider and the corresponding parameters are summarized in Table 2.  The DM mass is fixed at 5 GeV, the light mediator mass is fixed at 10 MeV, and the effective coupling for the heavy mediator is set to 1 × 10 −14 MeV −2 . Different colored lines correspond to different values of the detector resolution. For each resolution an appropriate threshold of E thr = 5 × σ was set.
We use a profile likelihood approach (similarly to e.g., Ref. [42]) to calculate the exclusion limits, where the test statistic, as a function of the free parameters of interest µ ⊂ {m χ , m Z , g l , g h }, is given by Here l b denotes the background level, which is set to the value that maximizes the likelihood function L with respect to the given value of µ in the numerator. The expressionl b is called the conditional maximum likelihood estimator of l b . In the denominator, L is maximized both in µ and l b , and so bothμ andl b are their unconditional maximum likelihood estimators. The background level is the only nuisance parameter over which we profile in this simplified approach. Hence, astrophysical quantities (like the galactic escape velocity) and detector specific quantities (like the threshold) are assumed to be precisely known.
A test statistic defined as in Eq. (3.1) follows a χ 2 distribution of degree n = dim(µ) [43]. One can thus exclude a specific set of parameters for the model at confidence level 1 − p if where CDF −1 χ 2 (n, x) is the inverse of the cumulative distribution function (quantile) for a χ 2 distribution with n degrees of freedom.
CRESST and COSINUS are being (will be) operated at the Gran Sasso underground laboratory LNGS. In the case of CRESST, these low background levels, together with a gross exposure of 5.594 kg days (net exposure of 3.46 kg days), lead to a lower number of resulting events, allowing us to work with an unbinned likelihood function. This is CRESST Table 2: Summary of experimental parameters used throughout this work. The constants in the CRESST-III background function are derived using data.
advantageous, as computational artifacts due to the choice of binning can be avoided. We thus use the following extended likelihood function, to account for the fact that the expected number of events N sample in a measured energy sample {E 1 , E 2 , . . . , E N sample } is itself a random variable: In the above expression f (E i |µ, l b ) denotes the non-normalized energy distribution of a certain DM model as measured by a certain detector, and is the total number of expected events that one achieves when integrating the distribution from the energy threshold E thr to the maximum energy E max (m χ ) = 2µ 2 N v 2 /m N . For the projection for COSINUS, the exposure and thus the number of expected events are significantly higher, calling for a binned likelihood approach to reduce the computation time. Thus the data have to be arranged in N bins bins with (n 1 , . . . , n N bins ) entries, where the bin width w bin has to be carefully attuned to the threshold and resolution of the experiment. The likelihood function is then given by When finding an expression for f (E i |µ, l b ), the finite energy resolution of the detector is taken into account by convolving the energy distribution as predicted by some DM model, with a Gaussian of width σ, the energy resolution of the detector. In general the detector energy resolution is an energy-dependent quantity. However, we assume that it is constant, as this is a valid simplification especially at the low recoil energies considered in this analysis. The final expression for the energy distribution thus amounts to where dR χ / dE and dR b / dE denote the differential recoil spectra for the DM model and the background model of the detector, respectively, while ε gives the exposure. Note here that the background component is not convolved with the resolution. The background for each of the experiments is modeled directly from the data and thus any detector-related properties are already considered. In the case of CRESST-III, the background is fitted to the data under the assumption of zero DM signal. The fit was done using a maximum likelihood approach and it was possible to achieve lowest values of the negative log-likelihood with a background containing a flat plus an exponentially decaying component (see also For the COSINUS projection the background distribution is fully determined by the flat model, which is used to generate the mock background sample. Within the analysis of real data released by the CRESST Collaboration, two more variables are needed in order to describe the expected energy distribution of the detector. The first variable is the signal survival probability ξ χ , which accounts for the loss of signal events from applied data selection criteria (cuts). Since the data from the CRESST-III experiment is taken from a so-called acceptance region in the light yield versus energy plane, 7 the probability of a signal event lying within the acceptance region, ξ χ,acc (E), has to be considered as well. Hence, for the CRESST-III data set the energy distribution amounts to We validate our approach by rederiving the official CRESST-III exclusion limits. The comparison between our derived limit and the official limit is shown in Fig. 7. In general, a good agreement between the two is seen. The small differences in the limit originate from the fact that the limits published by CRESST were calculated with Yellin's optimum interval method [44] and we are using a likelihood approach. A main difference between the two is the treatment of the background. The Yellin method makes no assumption on the background and solely exploits differences in the shape of the background and potential signal to set a limit on the cross section. In the likelihood approach we assume a flat background and an exponential for the excess at low energies. Their scaling parameters

DM-nucleon cross section [pb]
Limits CRESST-III likelihood approach Limits CRESST-III public data release Projection COSINUS likelihood approach are the only nuisance parameters of the likelihood. In this way, considering an exponentially decaying background component at low energies leads to stricter limits at low DM masses. Another factor that may lead to some differences between the results is the method with which the expected DM spectrum is determined. In this work this is done by means of the first term in Eq. (3.7), while CRESST uses a more complex procedure in their analysis described in Ref. [9]. The differences between these two methods were shown in Ref. [41]. Given the differences between the two approaches and a good agreement with the official limits, we consider our approach to be validated and use our implementation to derive limits for the parameters of the biportal model. In addition, we show the expected sensitivity for COSINUS in Fig. 7.

Results
We now present our results for the exclusion limits of the four free parameters in the biportal model to which direct-detection event rates are sensitive. In the following, we will assume that the DM relic density constraint is satisfied. 8 There are multiple planes in which the exclusion limits can be drawn, and each of them are instructive. We begin with the analysis in the m χ -g l plane for the real data from the CRESST-III experiment in Fig. 8. While the mass of the light mediator is fixed at 10 MeV in this figure, each panel corresponds to a different value of the heavy mediator coupling. We show the limits for the constructive (turquoise) and destructive (blue) case of the biportal model in comparison with the limit for a single light mediator (i.e. g h = 0) (dashed line). The excluded parameter space is shaded in all cases. Fig. 8 is best explained by starting with the bottom right panel, where g h has its smallest value of 1 × 10 −14 MeV −2 and the light mediator thus dominates the interaction  process. As expected, the biportal limits in both the constructive and destructive case agree with those for a single light mediator. For a slightly larger heavy-mediator coupling g h = 1 × 10 −13 MeV −2 , displayed in the bottom left panel, the biportal limits start to diverge from the single mediator case. While constructive interference leads to stronger limits -due to a higher number of predicted events -the destructive case gives weaker limits. In the upper two panels, the biportal model shows its differences in the allowed parameter space compared to the single-mediator case. For constructive interference the exclusion limits behave as expected: for increasing g h only lower values of g l are compatible with the background-only hypothesis and the limits thus get stronger. However, for the destructive case, the nonexcluded region is now not only limited from above, but for certain masses also limited from below.
This behavior can be understood as follows. Concentrating on the top two panels of Fig. 8, when g l is sufficiently large, the total differential cross section is dominated by the light mediator contribution. This together with the threshold of the detector leads to similar limits for the light DM mass around 0.2 GeV irrespective of the value for g h . As g l decreases, the total rate decreases, as explained also in Sec. 2.2. For intermediate values of g l the destructive interference plays the most important role, thus drastically reducing the total number of events obtained, and hence no limits are obtained. Reducing g l further renders the contribution of the light mediator negligible; however, sufficient discrimination  power is still present in the recoil spectrum due to the heavy mediator contribution g h . This leads to a somewhat artificial upper bound on g l . This is further confirmed by comparing the top two panels of Fig. 8, where increasing the heavy mediator coupling g h leads to a stronger upper bound on g l .
The parameter space of a single light mediator, biportal constructive and biportal destructive interference models are thus complementary to each other. While the constructive interference biportal model generically leads to stronger limits compared to the single light mediator models, limits from destructive interference biportal scenarios show a different constraining behavior. Generally this leads to two-sided limits and opens up some parameter space which is otherwise disfavored by single light mediator models. We stress that the biportal model studied in this work is capable of testing regions of parameter space that are otherwise not accessible to single-mediator models.
Turning our attention to the COSINUS experiment, in Fig. 9 we show the projected exclusion limits in the m χ − g l plane for fixed light mediator mass and heavy mediator coupling. Comparing this with the limits obtained for CaWO 4 , the qualitative similarity can be immediately seen. There are however crucial differences which arise from the experimental setup. The first and foremost difference is in the reach of the experiments for low DM masses. Given the low threshold, CRESST probes much smaller DM masses compared to COSINUS. Second, the expected COSINUS exposure is about 14 times larger compared to the CRESST data set used in this work, 9 which would make COSINUS probe much 9 The lower exposure of CRESST is not a general limitation of CRESST, but only true for the CRESST-    smaller values of g l . Finally, an interesting difference between the two can be observed in the limits for the destructive interference case. In the COSINUS experiment, the allowed parameter space corresponding to the strip between the lower and upper bounds on g l is III data set used in this work which is from a single detector only [9,41]. As is shown in Ref. [11], exposures on the scale of COSINUS or even higher would be within reach also for the small (24 g), low-threshold CRESST-III detectors.  In principle, tungsten would be a very good element, leading to a large mass gap in the composition of the CRESST target, but the recoil spectrum from scattering off this element is mostly below threshold and hence tungsten does not contribute for most of the parameter space, except potentially for very light DM masses.
In Fig. 10 exclusion limits (at 90% C.L.) using CRESST real data sample are again shown in the m χ − g l plane. Here the effective coupling of the heavy mediator is fixed, while each panel shows limits in the biportal model for different values of m V l . Excluded regions are again shaded. Overall, larger mediator masses lead to higher allowed values for g l . This can be explained by the inverse proportionality of the differential recoil rate (and also the total event rate) to the mediator mass. For lower m V l a light-mediator interaction model predicts a higher number of events and the limit on the coupling thus has to be lower. This also has a strong effect on the biportal model. Especially for m V l = 1 MeV in the upper left panel, the model is dominated by the light mediator and there are no noticeable visible differences between the constructive and destructive cases. Moreover, since the impact of the interference term is very low compared to the contribution from the light mediator, there is no lower limit on g l in the destructive case. In the next panel with m V l = 5 MeV, interference effects start to appear again. For higher DM masses the limits are very similar to those in the next panel for m V l = 10 MeV, which is identical to the already discussed top-right panel of Fig. 8. For lower DM masses the limits are however flatter, due to the lower mediator mass.
In Fig. 11, we show the limits derived using CRESST-III real data in the g l − m V l plane for a fixed dark matter mass of 5 GeV. The bounds for the constructive (cyan) and destructive (blue) cases of the biportal model are displayed together with the limit resulting from DM-nucleus interaction via a single light mediator (gray dashed). Higher values for the light mediator coupling are allowed for higher mediator masses. This behavior can be explained by the inverse mass proportionality in all light mediator terms contributing to the matrix element, which describes the DM-nucleus interaction. When going from higher to lower mediator masses, the limits get stronger as the predicted number of events increases. For mediator masses under 10 MeV the limits start to flatten. This is the mediator mass regime for which m V l < q 2 = 2m N E R and the DM-nucleus scattering can be described by a long-range interaction. The momentum transfer is then the governing quantity in the denominator of the intermediate vector-boson propagator and the mediator mass is expected to only have an impact at very low recoil energies for which m V l ≥ 2 m N E R . Also, the threshold of the experiment plays an important role here, as contributions from the strong exponential behavior of long-range interactions to the total rate at low recoil energies get cut off. As in the g l − m χ plane above, the constructive case of the biportal model gives stronger limits than the single light mediator case. The behavior of the destructive limit is also consistent with Fig. 8, showing an additional bound on the mediator coupling from below for higher mediator masses, where the rate is further reduced and the heavy mediator dominates the process. This results again in a band-like structure of the allowed parameter space.
In Fig. 12, we turn again to the COSINUS experiment, and discuss limits in the g l −m V l plane for a fixed DM mass of 5 GeV. This figure shows bounds for various values of the heavy mediator coupling to allow for comparison with the panels in Fig. 9. In the lower two panels one can again observe the single light mediator limiting case, while in the upper panels the overall behavior is comparable to the CRESST case.

Relic density
Finally, we briefly provide a discussion about relic density generation mechanisms and comment on the possibility of generating relic density within the parameter space we consider. We should remind ourselves that there are in principle seven free parameters in our theory: m V h , m V l , g Qh , g Ql , g χ h , g χ l , m χ . 10 In order to compute the relic density, we need to set all seven parameters.
Qualitatively speaking, there are three primary annihilation channels that contribute to the relic density, specifically, when m V h m χ as we consider here. The first two are the s-channel-mediated χχ → ff final states. The third is the t-channel diagram χχ → V l V l , which is proportional to g 2 χ l . For the s-channel diagrams there is a contribution from the interference term between the light and heavy mediators. Quantitatively, we implement the model in Micromegas 5.0.2 [45] to perform the calculations for the relic density. For completeness, we also include vector couplings to leptons, which we set to g Ql , g Qh for the respective interactions.
For concreteness and ease of comparison, we choose the same benchmarks used previously. In particular, we use m V l = 20 MeV. We set m V h = 500 GeV, g χ h = 1, g Qh = 0.02, such that the resulting effective heavy coupling g h = 10 −14 MeV −2 . We then vary m χ , g Ql .
In the left panel of Fig. 13, we show the resulting relic density contours in the g χ lm χ plane. When g χ l is large, the t-channel diagram contributes and the resulting relic density is driven by the χχ → V l V l process. Due to the large annihilation cross section in this region, the relic density is very small. As g χ l decreases, the contribution from the t-channel process decreases. For smaller values of g χ l the relic density is driven by two s-channel diagrams, χχ → V h /V l → ff . This cross section is directly proportional to g 2 χ l , and decreases as g χ l decreases. Correspondingly, the relic density increases. We obtain the required value of the relic density for g χ l ∼ 10 −4 -10 −5 . Decreasing g χ l further leads to overabundant dark matter.
We also note that there are strong constraints of such light mediators coupling to the SM particles [46]. For 10 MeV dark photons, the strongest constraints limit the dark photon-SM photon mixing parameter ≤ 10 −8 . Correspondingly, in the right panel of Fig. 13 (right panel) we show relic density contours for g Ql = g ll = 10 −8 with m V l = 10 MeV. We observe that for such suppressed couplings to the SM particles, the observed relic density is satisfied for g χl ∼ 10 −2 . Smaller values of g χl typically correspond to overabundant relic density due to suppressed cross sections.
Finally, we point out that there are several ways in which our results can be made compatible with observations without affecting the conclusions of our previous direct-detection analysis. For example, couplings to leptons could be increased while keeping the ones to quarks unchanged. Furthermore, in a complete theory, there could be additional particles to which dark matter can annihilate, therefore depleting the relic density. For example, introducing axial-vector couplings along with vector couplings would increase the χχ → ff cross section and thus reduce the relic density. Finally, the relic density can be diluted in the early Universe by additional entropy production.

Conclusions
Establishing the nature of particle DM and its interactions with SM fields is one of the outstanding open issues in particle physics and cosmology. In the context of direct DM detection, tremendous progress has been recently made on the front of low-threshold, highresolution experiments, which open up a new testing ground for DM models. Here we presented the case study of a simplified model with one heavy and one light dark vector mediator where interference effects at the level of DM-nucleus scattering amplitudes lead to sharp features in the differential recoil rates. Here, light and heavy are defined with respect to the momentum transfer at direct detection, which depends on the target and the DM mass. Our analysis showed that this kind of model is best explored at experiments with composite targets characterized by a large gap in the atomic mass of their components. Furthermore, we also demonstrated the importance of low-threshold, high-resolution experiments which maximally exploit the information contained in the shape of the recoil spectrum.
We applied a profile likelihood approach to analyze both public real data from the CRESST-III experiment and mock data based on the sensitivity of the future COSINUS experiment and derived limits on the free parameters of our biportal model with vector couplings. We found that CRESST-III already provides strong bounds, which will be further improved by COSINUS projections. To this end, we examined a large region of biportal model parameter space. We showed that CRESST-III public data are in general capable of probing DM masses as low as 0.2 GeV, and constrain light mediator couplings to be lower than 10 −10 for a light mediator mass of 10 MeV. The exact limits on light mediator couplings depend on the strength of the heavy mediator coupling as well. With its expected sensitivity the COSINUS experiment can set stronger limits on mediator couplings, albeit looses sensitivity to DM masses less than 1 GeV. This is due to a larger expected exposure, combination of target materials and higher threshold. We note that future CRESST exposures might also be larger, which we have not considered here. Once again, it is important to remember that the exact values of the allowed light mediator couplings depend on the heavy mediator couplings as well. We verified that constraints on the (overall) coupling weaken for light mediators with masses ¿ 50 MeV. In the future, it will be interesting to study further biportal simplified models and UV-complete theories characterized by sizable interference effects e.g., involving axial-vector couplings, and investigate the power of direct-detection experiments in determining viable regions in parameter space. Further complementary studies can involve constraints provided by other DM searches and DM relic density.
Our study provides a specific example of the importance that low-threshold, highresolution direct-detection experiments can have in guiding model building in particle DM by opening up the possibility to test novel features of recoil spectra and unexplored signatures.