#  Spin-singlet and spin-triplet pairing correlations in antiferromagnetically
coupled Kondo systems

Haoyu Hu  Department of Physics and Astronomy, Rice Center for Quantum
Materials, Rice University, Houston, Texas, 77005, USA  Ang Cai  Department of
Physics and Astronomy, Rice Center for Quantum Materials, Rice University,
Houston, Texas, 77005, USA  Lei Chen  Department of Physics and Astronomy,
Rice Center for Quantum Materials, Rice University, Houston, Texas, 77005, USA
Qimiao Si  Department of Physics and Astronomy, Rice Center for Quantum
Materials, Rice University, Houston, Texas, 77005, USA

###### Abstract

Recent experiments in quantum critical heavy fermion metals have pushed to the
fore the question about whether antiferromagnetic fluctuations can promote
both spin-singlet and spin-triplet superconductivity. Here we address the
issue through non-perturbative calculations in antiferromagnetically
correlated Kondo systems. We identify Kondo-destruction quantum critical
points in both the SU(2) symmetric and Ising-anisotropic cluster Bose-Fermi
Anderson models. The spin-singlet pairing correlations are significantly
enhanced near the quantum critical point in the SU(2) case; however, with
adequate but still realistic degree of Ising anisotropy, the spin-triplet
pairing correlations are competitive. Our results demonstrate that spin-flip
processes strengthen the spin-singlet pairing at the Kondo-destruction quantum
critical points, and point towards a way for antiferromagnetic correlations to
drive spin-triplet pairing. Further implications of our findings in the
broader context of strongly correlated superconductivity are discussed.

Quantum criticality has the potential to be a unifying theme in strongly
correlated metals Coleman and Schofield (2005); Paschen and Si (2021). A
quantum critical point (QCP) arises at the point of continuous transition
between two different types of ground states. Various types of novel order,
including unconventional superconductivity, may emerge in its vicinity Broun
(2008). The materials in which this physics may be relevant include the
cuprates Lee  _et al._ (2006), organic superconductors Powell and McKenzie
(2006), heavy fermion metals Kirchner  _et al._ (2020); Steglich and Wirth
(2016), and iron-based superconductors Si  _et al._ (2016).

Heavy fermion metals represent a particularly important case study for the
connection between quantum criticality and unconventional superconductivity.
There are a large number – about 50 – of them which superconduct. In many
cases, superconductivity develops out of a strange metal normal state, which
is often associated with a QCP at the border of an antiferromagnetic (AF)
order. Prototype examples are CeCu2Si2 Smidman  _et al._ (2018), and CeRhIn5
Park  _et al._ (2006); Knebel  _et al._ (2008); D. Thompson and Fisk (2012).
Recently, superconductivity has been discovered at ultra-low temperatures in
the canonical quantum critical heavy fermion metal, YbRh2Si2, which exhibits
AF order and a field-induced QCP Smidman  _et al._ (2018). In addition to the
superconducting phase at zero field, which has been interpreted as arising
from the intrinsic electronic quantum criticality unmasked by the hyperfine
coupling to nuclear spins Schuberth  _et al._ (2016), a new superconducting
phase occurs near the critical magnetic field Nguyen  _et al._ (2021). The
latter appears to have spin-triplet pairing, raising a profound question of
how antiferromagnetic correlations can drive spin-triplet pairing in quantum
critical Kondo systems.

To make progress, it is necessary to start from a proper description of the
quantum criticality in the normal state. It has been recognized that heavy
fermion quantum criticality can go beyond the Landau description based on the
slow fluctuations of a spin-density-wave (SDW) order Hertz (1976); Millis
(1993); Moriya (2012). Instead, a Kondo-destruction type of QCP Si  _et al._
(2001); Coleman  _et al._ (2001); Senthil  _et al._ (2004), which involves a
partial-Mott (delocalization-localization) transition and a small-to-large
Fermi surface jump, plays an important role Paschen  _et al._ (2004);
Gegenwart  _et al._ (2007); Friedemann  _et al._ (2010); Shishido  _et al._
(2005); Park  _et al._ (2006); Knebel  _et al._ (2008); Schröder  _et al._
(2000). Essentially, the Kondo effect is destroyed at the QCP by the dynamical
magnetic fluctuations. As such, the 4$f$ electrons that manifest as the Kondo-
driven delocalized composite fermions in the paramagnetic phase become
localized in the AF-ordered phase, leading to a sudden jump of the Fermi
surface from large to small.

The questions are, then, two-fold. First, does the Kondo destruction QCP
necessarily drive superconductivity? Second, under what conditions does the
quantum criticality promote spin-singlet or spin-triplet Cooper pairing?
Theoretically, this problem is challenging because in the quantum critical
regime there is no well-defined quasi-particle anywhere on the Fermi surface.
It is thus necessary to search for microscopic models which possess this kind
of QCP so that pairing correlations can be studied in a concrete setting.

Both of these questions can be addressed within cluster Bose-Fermi
Anderson/Kondo models (BFAM/BFKM). These models involve multiple local moments
coupled to each other through a Ruderman-Kittel-Kasuya-Yosida (RKKY)
interaction and simultaneously to both fermionic and bosonic baths. The
bosonic bath captures the dynamical magnetic fluctuations of a Kondo lattice
model, and formally develops from the later through a cluster extended
dynamical mean field approach (C-EDMFT) Pixley  _et al._ (2015a). In the Ising
limit (with infinite Ising anisotropy), the cluster BFAM has been studied by
the continuous-time quantum Monte Carlo (CT-QMC) method for its quantum
critical behavior and associated pairing correlations Pixley  _et al._
(2015b). Going beyond the pure Ising case, one can expect that the spin-flip
processes associated with the transverse component of the RKKY interaction
will significantly influence the nature of the pairing that develops. However,
this is a challenging problem, because the introduction of transverse exchange
couplings requires a triple expansion in the CT-QMC approach. Recently, we
have laid the groundwork by introducing an SU(2) CT-QMC method for a single-
impurity BFAM that involves a double expansion in the hybridization and
transverse Bose-Kondo coupling Cai and Si (2019) such that it accesses an
extended dynamical range in the quantum critical regime (see also, Ref.
Otsuki, 2013.)

In this Letter, we carry out the first study of an SU(2) symmetric cluster
BFAM as well as its spin-anisotropic counterpart with finite Ising anisotropy.
To do so, we develop the aforementioned triple-expansion scheme within the CT-
QMC approach. We identify a Kondo-destruction QCP in the SU(2) symmetric
model. Near the QCP, we find that spin-singlet pairing correlations are
significantly enhanced, more strongly than its Ising-anisotropic counterpart.
At the same time, the spin-triplet pairing correlations are suppressed, which
is to be contrasted with what happens in the realistically Ising-anisotropic
case. We understand both features from the physics that spin-flip processes
strengthen the spin-singlet pairing at the Kondo-destruction QCPs.

_Model and solution methods:_ The two-impurity SU(2) BFAM Hamiltonian is

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{i,\sigma}\epsilon_{d}d_{i\sigma}^{\dagger}d_{i\sigma}+\sum_{i}Un_{i\uparrow}n_{i\downarrow}+I\sum_{\alpha}S^{\alpha}_{1}S^{\alpha}_{2}$
(1) $\displaystyle+$ $\displaystyle\sum_{{\bf k},\sigma}\epsilon_{\bf
k}c_{{\bf k}\sigma}^{\dagger}c_{{\bf k}\sigma}+\sum_{i,{\bf
k},\sigma}\left(Ve^{i{\bf k\cdot r_{i}}}d_{i\sigma}^{\dagger}c_{{\bf
k}\sigma}+h.c.\right)$ $\displaystyle+$
$\displaystyle\sum_{p,\alpha}\omega_{p}{\phi^{\alpha}_{p}}^{\dagger}\phi^{\alpha}_{p}+g\sum_{p,\alpha}(S_{1}^{\alpha}-S_{2}^{\alpha})({\phi_{p}^{\alpha}}^{\dagger}+\phi_{p}^{\alpha})\,.$

The first line describes the interactions among the local degrees of freedom:
$d_{i\sigma}^{\dagger}$ creates an electron on site $i=1,2$ with spin
$\sigma=\uparrow,\downarrow$, $\epsilon_{d}$ is the onsite energy, $U$ is the
coulomb repulsion and $I$ is the RKKY interaction between the two impurities.
$S_{i}^{\alpha}=\frac{1}{2}\sum_{\mu,\nu}d^{\dagger}_{i\mu}\sigma^{\alpha}_{\mu\nu}d_{i\nu}$
is the spin operator and $\sigma^{\alpha}_{\mu\nu}$ are the Pauli matrices.
The second line describes the fermionic bath and its coupling with the
impurities: $c^{\dagger}_{{\bf k}\sigma}$ creates a conduction electron with
wave vector ${\bf k}$, spin $\sigma$ and energy $\epsilon_{\bf k}$, which
hybridizes with the local $d$ electron through the matrix element $V$. The
third line describes the 3-component vector bosonic bath and its coupling to
the anti-ferro component of the two local moments through coupling constant
$g$. $\alpha$ is the spin index and it is summed over $x,y$ and $z$ such that
the Hamiltonian is SU(2) invariant. We will also consider the Ising-
anisotropic version of the model (see below).

For the fermionic bath, we consider a flat density of states
$\rho_{F}(\epsilon)=\sum_{\bf k}\delta(\epsilon-\epsilon_{\bf
k})=\rho_{0}\Theta(|D-\epsilon|)$, with a hybridization function
$\Gamma(\epsilon)=\Gamma_{0}\Theta(|D-\epsilon|)$ and
$\Gamma_{0}=\pi\rho_{0}V^{2}$. In addition, the impurity location ${\bf
r_{1}}$ and ${\bf r_{2}}$ are set to be infinitely separated to avoid double
counting of the RKKY interactions. The bosonic bath takes a sub-ohmic density
of states specified by
$\rho_{B}(\omega)=\sum_{p}\delta(\omega-\omega_{p})=K_{0}\omega^{s}e^{-\omega/\Lambda}\Theta(\omega)$.
The prefactor $\rho_{0}$ and $K_{0}$ are determined by the normalization
conditions: $\int_{-D}^{D}\rho_{F}(\epsilon)d\epsilon=1$ and
$\int_{0}^{\infty}\rho_{B}(\omega)d\omega=1$. The BFAM can be solved by a
recently developed sign-problem-free CT-QMC method. The triple expansion
method builds on prior work that has dealt with various double expansions
Steiner  _et al._ (2015); Otsuki (2013); Cai and Si (2019). Further details of
the method are described in the Supplementary Material. In the calculation, we
set $D=1$, $\Lambda=1$, and preserve particle-hole symmetry by choosing
$U=-2\epsilon_{d}=0.1$. In addition, we take a sub-ohmic bosonic bath with
exponent $s=0.6$ and fix the hybridization parameter $\Gamma_{0}$ at
$\Gamma_{0}=0.5$. The RKKY interactions $I$ and bosonic coupling $g$ are
tuning parameters to explore phase digram and search for QCPs.

Figure 1: Phase diagram for the two-impurity SU(2) Bose-Fermi Anderson model
determined via CT-QMC in the $I$ vs $g^{2}$ plane at fixed hybridization
parameter $\Gamma_{0}=0.5$. Blue circles represent QCPs from Kondo screened
phase to impurity singlet (IS) phase governed by the critical point at $g=0$.
Orange circles represent QCPs from IS phase to antiferromagnetic local moment
(LM) phase. Red circles represent QCPs from the Kondo screened phase to the LM
phase, whose critical behavior is governed by the Kondo-destruction fixed
point at $I=0$. Dashes line represent the path
$(g^{2},I)=(2.22\lambda,\lambda)$ taken to study the pairing correlation,
which mimics the trajectory in self-consistent C-EDMFT calculations.

_Quantum critical properties:_ Fig. 1 shows the $T=0$ phase diagram of SU(2)
symmetric BFAM in the $g^{2}-I$ parameter space obtained from CT-QMC. In the
absence of the bosonic bath, the model reduces to the two impurity Anderson
model Jones and Varma (1987), which has a single QCP from Kondo screened phase
at small RKKY interaction to an impurity singlet (IS) phase at large RKKY
interaction. Turning on the bosonic coupling $g$, we find that the transition
persists with the same critical behavior but at a reduced critical value of
$I$. At large bosonic coupling, there will be a second transition from the IS
phase to an antiferromagnetic local moment (LM) phase. For $I$ smaller than
$0.3$, the IS phase disappears and we have a direct transition between a Kondo
screened phase and the LM phase. This turns out to be a Kondo destruction QCP
that has not been studied before.

Given its “beyond-Landau” nature, we characterize the Kondo-destruction QCP
with fidelity susceptibility that directly evaluates the changes of
wavefunctions and detects quantum phase transition without introducing any
order parameters Zanardi and Paunković (2006). The fidelity susceptibility
$\chi_{F}$ is defined as the second derivative of fidelity’s logarithm, where
the fidelity measures the distance between the ground-state wavefunctions at
two different values of parameters $V,V+\delta V$. This definition can be
extended to finite temperature: $\chi_{F}=\int_{0}^{\beta/2}\left(\langle{\cal
T_{\tau}}H_{V}(\tau)H_{V}\rangle-\langle H_{V}\rangle^{2}\right)\tau d\tau$,
where $H_{V}$ includes all the terms in the Hamiltonian that are proportional
to the hybridization strength $V$ Albuquerque  _et al._ (2010); Wang  _et al._
(2015). Due to the sharp changes of ground-state wavefunctions near the QCP,
we observe a singular behaviors of $\chi_{F}$:

$\displaystyle\chi_{F}(\beta,g)/\beta$ $\displaystyle=$
$\displaystyle\beta^{2/\nu}\tilde{\chi}\left(\beta^{1/\nu}(g-g_{c})/g_{c}+A/\beta^{\phi/\nu}\right)\,,$
(2)

with $g_{c}=0.66(2)$ and $\nu^{-1}=0.31(6)$. The critical value and critical
exponent are obtained by fitting the numerical data as shown in Fig. 2 (a),
(b). The scaling behavior is consistent with a continuous transition and,
moreover, suggests that $V$ is a relevant perturbation at QCP. Thus, the
transition is indeed a Kondo destruction QCP. Besides, we study the Binder
cumulants and correlation length of the spin correlation function
(Supplementary Material). The results are consistent with the conclusions
drawn above from the fidelity susceptibility.

Figure 2: (a) Fidelity susceptibility $\chi_{F}$ vs. bosonic coupling $g$,
showing divergence near the critical coupling $g_{c}=0.665(5)$. (b) Finite
size scaling of $\chi_{F}$ according to Eq.(2), with $\nu^{-1}=0.31$,
$g_{c}=0.66$. (c) Staggered spin susceptibility $\chi_{s}(T)$ in the Kondo
screened phase (green squares), at the critical point (black circles), and in
the LM phase (blue diamonds). The red and grey dashed lines are fits to
$\chi_{s}(T)$ at the critical point $g=g_{c}$ and in the LM phase,
respectively.

To further characterize the QCP, we analyze the static spin susceptibility,
defined as $\chi_{s}(T)=\int_{0}^{\beta}\chi_{s}(\tau)d\tau$. As seen in Fig.
2 (c), $\chi_{s}(T)$ diverges as $\chi_{s}(T)\sim T^{-x}$ with $x=0.60(1)$ at
the QCP. The result is consistent with the relation $x=s$, suggesting that the
critical behavior in the spin channel is completely determined by the
properties of the bosonic bath. Note that this is one of the necessary
conditions for the Kondo-destruction QCP to be realized in the lattice case Si
_et al._ (2001).

Figure 3: (a) Static singlet pairing susceptibility $\chi^{(s,0)}(T)$ along
the cut $\lambda$: $(g^{2},I)=(2.22\lambda,\lambda)$ as shown in Fig. 1. The
Kondo destruction QCP is located at $\lambda_{c}=0.2$. (b) Imaginary time
pairing correlation function $\chi^{(s,0)}(\tau)$ at the QCP, consistent with
a $1/\tau^{2-\eta}$ decay with $\eta=0.6$. The dashed orange line is a fit
according to $\chi^{(s,0)}(\tau)\sim\left(\pi T/sin(\pi\tau T)\right)^{1.4}$.

_Pairing susceptibilities:_ Within our 2-site construction we can study four
different types of inter-site pairing, specified by the following pairing
operators: $\Delta_{(t,\pm
1)}^{\dagger}=\Delta^{\dagger}_{\uparrow\uparrow/\downarrow\downarrow}$,
$\Delta_{(t,0)}^{\dagger}=\left(\Delta^{\dagger}_{\uparrow\downarrow}+\Delta^{\dagger}_{\downarrow\uparrow}\right)/\sqrt{2}$,
and
$\Delta_{(s,0)}^{\dagger}=\left(\Delta^{\dagger}_{\uparrow\downarrow}-\Delta^{\dagger}_{\downarrow\uparrow}\right)/\sqrt{2}$,
with
$\Delta_{\sigma\sigma^{\prime}}^{\dagger}=d_{1\sigma}^{\dagger}d_{2\sigma^{\prime}}^{\dagger}$,
$\\{\sigma,\sigma^{\prime}\\}\in\\{\uparrow,\downarrow\\}$. They are,
respectively, the pairing operators in the triplet channel with
$m^{z}_{tot}=\pm 1$, $m^{z}_{tot}=0$ and the singlet channel with
$m^{z}_{tot}=0$, where $m^{z}_{tot}$ is the $z$ component of the total spin of
the cooper pair. Accordingly we can define the dynamic and static pairing
correlation functions respectively: $\chi^{(\alpha,m)}(\tau)=\langle{\cal
T_{\tau}}\Delta_{(\alpha,m)}^{\dagger}(\tau)\Delta_{(\alpha,m)}\rangle$,
$\chi^{(\alpha,m)}(T)=\int_{0}^{\beta}\chi^{(\alpha,m)}(\tau)d\tau$. The
details of the calculational procedure is given in the Supplementary Material.

We study the behavior of the pairing susceptibility following a particular cut
$(g^{2},I)=(2.22\lambda,\lambda)$ illustrated in Fig. 1. This cut mimics the
trajectory of the C-EDMFT self-consistent calculation, and will meet the line
of Kondo destruction QCPs at $\lambda=\lambda_{c}=0.2$. As shown in Fig. 3(a),
the pairing susceptibility in the singlet channel $\chi^{(s,0)}(T)$ increases
as we move from the Kondo screened phase to the QCP, and is peaked slightly
before the transition. There is also a large enhancement near the QCP as we
lower the temperature: $\chi^{(s,0)}(\tau)$ acquires a large anomalous
dimension at the critical fixed point. In Fig. 3(b) we plot
$\chi^{(s,0)}(\tau)$ against the conformal form
$\chi^{(s,0)}(\tau)\sim\left(\pi T/\sin(\pi\tau T)\right)^{2-\eta}$. We
observe a collapse across different temperatures, which suggests that
$\chi^{(s,0)}(\tau)\sim 1/\tau^{2-\eta}$, with an anomalous exponent
$\eta\simeq 0.6$. At the same time, we find that the pairing susceptibility in
the triplet channel is strongly suppressed and does not acquire any anomalous
dimension.

_SU(2) vs. Ising-anisotropic model_ For the SU(2) model, Fig. 4(a) displays
the evolution of the pairing susceptibilities across the QCP. It clearly shows
the enhancement of the spin-singlet pairing correlation near the QCP, and the
concomitant suppression of the spin-triplet pairing correlation in all the
three (degenerate) channels.

We now consider an Ising-anisotropic model. The extreme Ising-anisotropic
model of Ref. Pixley  _et al._ , 2015b corresponds to having only the Ising
component for both the bosonic-Kondo coupling and RKKY interaction. More
realistically, we consider the model whose effective RKKY interaction is
written in an anisotropic form:

$\displaystyle
J_{z}S^{z}_{1}S^{z}_{2}+J_{p}(S^{x}_{1}S^{x}_{2}+S^{y}_{1}S^{y}_{2})\,.$ (3)

Here $J_{z}$ and $J_{p}$ are respectively the longitudinal and spin-flip RKKY
exchange interactions. The details of the model are described in the
Supplementary Material.

We consider the parameters such that the Ising anisotropy near the QCP has a
realistic $J_{p}/J_{z}\approx 0.2$. Fig. 4(b) shows that the pairing
correlation in the $S_{z}=0$ triplet channel is enhanced near the QCP, along
with the spin-singlet pairing correlation. (By contrast, the pairing
correlations in the $S_{z}=\pm 1$ triplet channels are suppressed.) In the
case of extreme Ising anisotropy, with $J_{p}=0$, the two dominant pairing
susceptibilities are exactly equal to each other. Here, the pairing
susceptibility in the $S_{z}=0$ triplet channel is almost equal to its
counterpart in the spin-singlet channel. An analysis of the Landau free-energy
functional shows that, a Zeeman coupling will further enhance the $S_{z}=0$
spin-triplet channel over the spin-singlet channel, thereby making the spin-
triplet pairing dominate Hu et al. (2021).

Figure 4: Static pairing susceptibilities of the SU(2) symmetric model (a) and
the Ising anistropic model (b) at $\beta=400$. (a) SU(2) symmetric model,
where the trajectory is the same as in Fig. 1, with $\lambda_{c}=0.2$. The
degenerate pairing susceptibilities in three triplet channels,
$\chi^{(t,0)},\chi^{(t,\pm 1)}$, are suppressed near the QCP. (b) The Ising
anisotropic model, with anisotropy $J_{p}/J_{z}=0.2$. We take a cut along
$\lambda:(g,I)=(\lambda,2.688\lambda)$, with $\lambda_{c}\sim 0.45$. The
pairing susceptibility in the triplet $S^{z}=0$ channel, $\chi^{(t,0)}$, is
similarly enhanced near the QCP as for the singlet channel $\chi^{(s,0)}$.

_Discussion and Summary:_ Our findings shed light on the recent experiment in
YbRh2Si2 near its field-driven QCP Nguyen  _et al._ (2021). YbRh2Si2 has
antiferromagnetic order at the ambient conditions. It has a strong
$xy$-anisotropy in the spin space, which is turned to Ising-anisotropic by the
applied magnetic field. Our results for the Ising-anisotropic cluster BFAM
with an AF RKKY interaction suggest that the $S_{z}=0$ spin-triplet pairing
channel will be competitive against the spin-singlet pairing near the Kondo-
destruction QCP, and can be further turned into the dominant pairing channel
by the Zeeman coupling of the critical magnetic field. This provides a natural
understanding of the surprising experimental results.

By contrast, when the spin anisotropy is relatively weak, our results derived
from the SU(2) symmetric cluster BFAM suggests that the spin-singlet pairing
dominates. This result lays the groundwork for the understanding of the
superconductivity that has been experimentally observed in a host of
antiferromagnetically correlated heavy fermion metals, such as CeRhIn5 and
CeCu2Si2.

Our results also shed light on other classes of strongly correlated metals.
For example, UTe2 has a strong Ising anisotropy in its spin response Ran  _et
al._ (2019). Recent inelastic neutron scattering experiments Duan  _et al._
(2021) show that both the normal and superconducting states are
antiferromagnetically correlated; this raises the exciting possibility that,
here too, antiferromagnetic correlations are primarily responsible for the
potential spin-triplet pairing.

Going beyond the heavy-fermion context, the electronic delocalization-
localization effect has been implicated by experiments in both the hole-doped
cuprates Badoux  _et al._ (2016) and organic charge-transfer salts Oike  _et
al._ (2015). Recently, based on analyses of related dynamical equations, it
has been suggested that beyond-Landau quantum criticality that bears
similarities to the Kondo-destruction QCP operates in doped Mott insulators
Chowdhury  _et al._ (2021). These calculations have stayed at the spin-
rotational invariant level. Our results suggest that in these models too, it
would be instructive to analyze the role of spin anisotropy that may arise in
doping Mott insulators such as Sr2IrO4.

To summarize, we have studied a cluster Bose-Fermi Anderson model as an
effective Kondo system for beyond-Landau quantum criticality. In both the
SU(2) symmetric and Ising anisotropic models, we identify Kondo-destruction
quantum critical points. The spin-singlet pairing correlations are more
strongly enhanced near the QCP in the SU(2) symmetric case than its Ising-
anisotropic counterpart. This reveals that spin-flip processes strengthen the
spin-singlet pairing. In accordance with this insight, we find the $S_{z}=0$
spin-triplet pairing channel to be competitive in the realistically Ising-
anisotropic version of the model. The results provide a natural understanding
of the recent experiments on the superconductivity near the field-driven
quantum critical point of the prominent heavy fermion system YbRh2Si2, in a
way that incorporates the salient quantum critical characteristics in its
normal state. Finally, our findings bring out new insights pertinent to the
unconventional superconductivity in a broad range of strongly correlated
metals.

_Acknowledgements._ We thank P. C. Dai, K. Ingersent, S. Paschen, J. H.
Pixley, F. Steglich, J. D. Thompson and J.-X. Zhu for useful discussions. The
work was in part supported by the National Science Foundation under Grant No.
DMR-1920740 (H.H. and L.C.), the Air Force Office of Scientific Research under
Grant No. FA9550-21-1-0356 (Q.S.), the Robert A. Welch Foundation Grant No.
C-1411 (A.C.), the Data Analysis and Visualization Cyberinfrastructure funded
by NSF under grant OCI-0959097 and an IBM Shared University Research (SUR)
Award at Rice University, and the Extreme Science and Engineering Discovery
Environment (XSEDE) by NSF under Grant No. DMR170109. Q.S. acknowledges the
hospitality of the Aspen Center for Physics (NSF under Grant No. PHY-1607611).

## References

  * Coleman and Schofield (2005) Piers Coleman and Andrew J. Schofield, “Quantum criticality,” Nature 433, 226–229 (2005).
  * Paschen and Si (2021) Silke Paschen and Qimiao Si, “Quantum phases driven by strong correlations,” Nat. Rev. Phys. 3, 9 (2021).
  * Broun (2008) D. M. Broun, “What lies beneath the dome?” Nature Physics 4, 170–172 (2008).
  * Lee  _et al._ (2006) Patrick A. Lee, Naoto Nagaosa, and Xiao-Gang Wen, “Doping a Mott insulator: Physics of high-temperature superconductivity,” Rev. Mod. Phys. 78, 17–85 (2006).
  * Powell and McKenzie (2006) B. J. Powell and Ross H. McKenzie, “Strong electronic correlations in superconducting organic charge transfer salts,” J. Phys. Condens. Matter 18, R827–R866 (2006).
  * Kirchner  _et al._ (2020) Stefan Kirchner, Silke Paschen, Qiuyun Chen, Steffen Wirth, Donglai Feng, Joe D. Thompson, and Qimiao Si, “Colloquium: Heavy-electron quantum criticality and single-particle spectroscopy,” Rev. Mod. Phys. 92, 011002 (2020).
  * Steglich and Wirth (2016) Frank Steglich and Steffen Wirth, “Foundations of heavy-fermion superconductivity: lattice Kondo effect and Mott physics,” Rep. Prog. Phys 79, 084502 (2016).
  * Si  _et al._ (2016) Qimiao Si, Rong Yu, and Elihu Abrahams, “High-temperature superconductivity in iron pnictides and chalcogenides,” Nat. Rev. Mater 1, 1–15 (2016).
  * Smidman  _et al._ (2018) M. Smidman, O. Stockert, J. Arndt, G. M. Pang, L. Jiao, H. Q. Yuan, H. A. Vieyra, S. Kitagawa, K. Ishida, K. Fujiwara, T. C. Kobayashi, E. Schuberth, M. Tippmann, L. Steinke, S. Lausberg, A. Steppke, M. Brando, H. Pfau, U. Stockert, P. Sun, S. Friedemann, S. Wirth, C. Krellner, S. Kirchner, E. M. Nica, R. Yu, Q. Si, and F. Steglich, “Interplay between unconventional superconductivity and heavy-fermion quantum criticality: CeCu2Si2 versus YbRh2Si2,” Philos. Mag. 98, 2930–2963 (2018).
  * Park  _et al._ (2006) Tuson Park, F. Ronning, H. Q. Yuan, M. B. Salamon, R. Movshovich, J. L. Sarrao, and J. D. Thompson, “Hidden magnetism and quantum criticality in the heavy fermion superconductor CeRhIn5,” Nature 440, 65–68 (2006).
  * Knebel  _et al._ (2008) G. Knebel, D. Aoki, J.-P. Brison, and J. Flouquet, “The quantum critical point in $\mathrm{CeRhIn_{5}}$: a resistivity study,” J. Phys. Soc. Jpn. 77, 114704 (2008).
  * D. Thompson and Fisk (2012) Joe D. Thompson and Zachary Fisk, “Progress in heavy-fermion superconductivity: $\mathrm{Ce}$115 and related materials,” J. Phys. Soc. Jpn. 81, 011002 (2012).
  * Schuberth  _et al._ (2016) Erwin Schuberth, Marc Tippmann, Lucia Steinke, Stefan Lausberg, Alexander Steppke, Manuel Brando, Cornelius Krellner, Christoph Geibel, Rong Yu, Qimiao Si, and Frank Steglich, “Emergence of superconductivity in the canonical heavy-electron metal YbRh2Si2,” Science 351, 485–488 (2016).
  * Nguyen  _et al._ (2021) D. H. Nguyen, A. Sidorenko, M. Taupin, G. Knebel, G. Lapertot, E. Schuberth, and S. Paschen, “Superconductivity in an extreme strange metal,” Nat. Commun 12, 4341 (2021).
  * Hertz (1976) John A. Hertz, “Quantum critical phenomena,” Phys. Rev. B 14, 1165–1184 (1976).
  * Millis (1993) A. J. Millis, “Effect of a nonzero temperature on quantum critical points in itinerant fermion systems,” Phys. Rev. B 48, 7183–7196 (1993).
  * Moriya (2012) Toru Moriya,  _Spin fluctuations in itinerant electron magnetism_ , Vol. 56 (Springer Science & Business Media, 2012).
  * Si  _et al._ (2001) Qimiao Si, Silvio Rabello, Kevin Ingersent, and J Lleweilun Smith, “Locally critical quantum phase transitions in strongly correlated metals,” Nature 413, 804 (2001).
  * Coleman  _et al._ (2001) P. Coleman, C. Pépin, Qimiao Si, and R. Ramazashvili, “How do Fermi liquids get heavy and die?” J. Phys. Condens. Matter 13, R723 (2001).
  * Senthil  _et al._ (2004) T. Senthil, M. Vojta, and S. Sachdev, “Weak magnetism and non-Fermi liquids near heavy-fermion critical points,” Phys. Rev. B 69, 035111 (2004).
  * Paschen  _et al._ (2004) S. Paschen, T. Lühmann, S. Wirth, P. Gegenwart, O. Trovarelli, C. Geibel, F. Steglich, P. Coleman, and Q. Si, “Hall-effect evolution across a heavy-fermion quantum critical point,” Nature 432, 881–885 (2004).
  * Gegenwart  _et al._ (2007) P. Gegenwart, T. Westerkamp, C. Krellner, Y. Tokiwa, S. Paschen, C. Geibel, F. Steglich, E. Abrahams, and Q. Si, “Multiple energy scales at a quantum critical point,” Science 315, 969–971 (2007).
  * Friedemann  _et al._ (2010) Sven Friedemann, Niels Oeschler, Steffen Wirth, Cornelius Krellner, Christoph Geibel, Frank Steglich, Silke Paschen, Stefan Kirchner, and Qimiao Si, “Fermi-surface collapse and dynamical scaling near a quantum-critical point,” Proc. Natl. Acad. Sci. U.S.A 107, 14547–14551 (2010).
  * Shishido  _et al._ (2005) Hiroaki Shishido, Rikio Settai, Hisatomo Harima, and Yoshichika Ōnuki, “A drastic change of the Fermi surface at a critical pressure in CeRhIn5: dHvA study under pressure,” J. Phys. Soc. Jpn 74, 1103–1106 (2005).
  * Schröder  _et al._ (2000) A. Schröder, G. Aeppli, R. Coldea, M. Adams, O. Stockert, H. v Löhneysen, E. Bucher, R. Ramazashvili, and P. Coleman, “Onset of antiferromagnetism in heavy-fermion metals,” Nature 407, 351 (2000).
  * Pixley  _et al._ (2015a) J. H. Pixley, Ang Cai, and Qimiao Si, “Cluster extended dynamical mean-field approach and unconventional superconductivity,” Phys. Rev. B 91, 125127 (2015a).
  * Pixley  _et al._ (2015b) J. H. Pixley, Lili Deng, Kevin Ingersent, and Qimiao Si, “Pairing correlations near a Kondo-destruction quantum critical point,” Phys. Rev. B 91, 201109 (2015b).
  * Cai and Si (2019) Ang Cai and Qimiao Si, “Bose-Fermi Anderson model with SU(2) symmetry: Continuous-time quantum Monte Carlo study,” arXiv preprint arXiv:1902.10094 (2019).
  * Otsuki (2013) Junya Otsuki, “Spin-boson coupling in continuous-time quantum Monte Carlo,” Phys. Rev. B 87, 125102 (2013).
  * Steiner  _et al._ (2015) Karim Steiner, Yusuke Nomura, and Philipp Werner, “Double-expansion impurity solver for multiorbital models with dynamically screened $U$ and $J$,” Phys. Rev. B 92, 115123 (2015).
  * Jones and Varma (1987) B. A. Jones and C. M. Varma, “Study of two magnetic impurities in a Fermi gas,” Phys. Rev. Lett. 58, 843–846 (1987).
  * Zanardi and Paunković (2006) Paolo Zanardi and Nikola Paunković, “Ground state overlap and quantum phase transitions,” Phys. Rev. E 74, 031123 (2006).
  * Albuquerque  _et al._ (2010) A. Fabricio Albuquerque, Fabien Alet, Clément Sire, and Sylvain Capponi, “Quantum critical scaling of fidelity susceptibility,” Phys. Rev. B 81, 064418 (2010).
  * Wang  _et al._ (2015) Lei Wang, Ye-Hua Liu, Jakub Imriška, Ping Nang Ma, and Matthias Troyer, “Fidelity susceptibility made simple: A unified quantum Monte Carlo approach,” Phys. Rev. X 5, 031007 (2015).
  * Hu et al. (2021) Haoyu Hu et al., “unpublsihed,” (2021).
  * Ran  _et al._ (2019) Sheng Ran, Chris Eckberg, Qing-Ping Ding, Yuji Furukawa, Tristin Metz, Shanta R. Saha, I-Lin Liu, Mark Zic, Hyunsoo Kim, Johnpierre Paglione, and Nicholas P. Butch, “Nearly ferromagnetic spin-triplet superconductivity,” Science 365, 684–687 (2019).
  * Duan  _et al._ (2021) Chunruo Duan, R. E. Baumbach, Andrey Podlesnyak, Yuhang Deng, Camilla Moir, Alexander J. Breindel, M. Brian Maple, E. M. Nica, Qimiao Si, and Pengcheng Dai, “Resonance from antiferromagnetic spin fluctuations for superconductivity in UTe2,” Nature (2021).
  * Badoux  _et al._ (2016) S. Badoux, W. Tabis, F. Laliberté, G. Grissonnanche, B. Vignolle, D. Vignolles, J. Béard, D. A. Bonn, W. N. Hardy, R. Liang, N. Doiron-Leyraud, Louis Taillefer, and Cyril Proust, “Change of carrier density at the pseudogap critical point of a cuprate superconductor,” Nature 531, 210 (2016).
  * Oike  _et al._ (2015) H. Oike, K. Miyagawa, H. Taniguchi, and K. Kanoda, “Pressure-induced Mott transition in an organic superconductor with a finite doping level,” Phys. Rev. Lett. 114, 067002 (2015).
  * Chowdhury  _et al._ (2021) Debanjan Chowdhury, Antoine Georges, Olivier Parcollet, and Subir Sachdev, “Sachdev-Ye-Kitaev models and beyond: A window into non-Fermi liquids,” arXiv:2109.05037 (2021).
  * Werner and Millis (2007) Philipp Werner and Andrew J. Millis, “Efficient dynamical mean field simulation of the Holstein-Hubbard model,” Phys. Rev. Lett. 99, 146404 (2007).
  * Lang and Firsov (1962) I. J. Lang and Yu A. Firsov, “Kinetic theory of small mobility semiconductors,” Sov. Phys. JETP 16, 1301–1317 (1962).
  * Metropolis  _et al._ (1953) Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller, “Equation of state calculations by fast computing machines,” J. Chem. Phys 21, 1087–1092 (1953).
  * Hastings (1970) W. K. Hastings, “Monte Carlo sampling methods using Markov chains and their applications,” Biometrika 57, 97–109 (1970).
  * Binder (1981) Kurt Binder, “Finite size scaling analysis of Ising model block distribution functions,” Z. Phys., B Condens. matter 43, 119–140 (1981).
  * Sandvik (2010) Anders W. Sandvik, “Computational studies of quantum spin systems,” in  _AIP Conf. Proc_ , Vol. 1297 (AIP, 2010) pp. 135–338.
  * Gull  _et al._ (2007) Emanuel Gull, Philipp Werner, Andrew Millis, and Matthias Troyer, “Performance analysis of continuous-time solvers for quantum impurity models,” Phys. Rev. B 76, 235123 (2007).
  * Hoshino and Werner (2016) Shintaro Hoshino and Philipp Werner, “Electronic orders in multiorbital Hubbard models with lifted orbital degeneracy,” Phys. Rev. B 93, 155161 (2016).
  * Zhu and Zhu (2011) Lijun Zhu and Jian-Xin Zhu, “Coherence scale of coupled Anderson impurities,” Phys. Rev. B 83, 195103 (2011).

##  I Continuous-time Quantum Monte Carlo Method

In the continuous-time quantum Monte Carlo method, we perform a unitary
transformation Werner and Millis (2007); Lang and Firsov (1962) to remove the
$z$ component spin-boson coupling, using the generator ${\cal
S}=g(S_{1}^{z}-S_{2}^{z})\sum_{p}\frac{1}{\omega_{p}}(\phi^{z\dagger}_{p}+\phi^{z}_{p})$.
The new Hamiltonian $\tilde{H}=e^{\cal S}He^{-\cal S}$ can be separated into
two parts $\tilde{H}_{0}+\tilde{H}_{1}$, where
$\tilde{H}_{0}=\tilde{\epsilon}_{d}\tilde{d}_{i\sigma}^{\dagger}\tilde{d}_{i\sigma}+\sum_{i}\tilde{U}n_{i\uparrow}n_{i\downarrow}+\tilde{I}S_{1}^{z}S_{2}^{z}+\sum_{{\bf
k},\sigma}\epsilon_{\bf k}c_{{\bf k}\sigma}^{\dagger}c_{{\bf
k}\sigma}+\sum_{p,\alpha}\omega_{p}{\phi^{\alpha}_{p}}^{\dagger}\phi^{\alpha}_{p}$,
$\tilde{H}_{1}=\sum_{i}\left(H_{c,i}+(-1)^{i+1}H_{b,i}\right)+H_{I}$ \+ h.c.,
with $H_{c,i}=\sum_{\bf{k},\sigma}Ve^{i{\bf k\cdot
r_{i}}}\tilde{d}_{i\sigma}^{\dagger}c_{{\bf k},\sigma}$,
$H_{b,i}=g\tilde{S}_{i}^{+}\Phi^{-}$, and
$H_{I}=(I/2)\tilde{S}_{1}^{+}\tilde{S}_{2}^{-}$. Here
$\tilde{d}_{i,\sigma}=e^{\cal S}{d}_{i,\sigma}e^{-\cal
S}=d_{i,\sigma}e^{-\frac{g}{2}s_{\sigma}p_{i}\sum_{p}\frac{1}{\omega_{p}}(\phi_{p}^{z\dagger}-\phi_{p}^{z})}$
is the dressed fermionic operator, with $s_{\sigma}=\pm 1$ for
$\sigma=\uparrow/\downarrow$ and $p_{i}=\pm 1$ for $i=1,2$.
$\tilde{S}^{+}_{i}=\tilde{d}_{i\uparrow}^{\dagger}\tilde{d}_{i\downarrow}$,
$\Phi^{-}=1/\sqrt{2}\sum_{p}\left((\phi_{p}^{x\dagger}+\phi_{p}^{x})-i(\phi_{p}^{y\dagger}+\phi_{p}^{y})\right)$.
$\tilde{U}$, $\tilde{\epsilon_{d}}$ and $\tilde{I}$ are the renormalized
coupling constant, given by $\tilde{U}=U+(g^{2}/2)\sum_{p}1/\omega_{p}^{2}$,
$\tilde{\epsilon_{d}}=\epsilon_{d}-(g^{2}/4)\sum_{p}1/\omega_{p}^{2}\textbf{
}$ and $\tilde{I}=I+2g^{2}\sum_{p}1/\omega_{p}^{2}$.

Now, we are able to express the partition function ${\cal Z}$ in terms of an
expansion of $\tilde{H}_{1}$ under the interaction representation of
$\tilde{H}_{0}$,

$\displaystyle{\cal Z}$ $\displaystyle=$
$\displaystyle\sum_{n^{c}_{i},n^{b}_{i},\bar{n}^{c}_{i},\bar{n}^{b}_{i},n^{I},\bar{n}^{I}}(-1)^{p}\times\text{Tr}[e^{-\beta\tilde{H}_{0}}{\cal
T_{\tau}}$ (4) $\displaystyle\times$
$\displaystyle\frac{1}{n^{I}!}\left(\int_{0}^{\beta}H_{I}(\tau)d\tau\right)^{n^{I}}\
\frac{1}{\bar{n}^{I}!}\left(\int_{0}^{\beta}H^{\dagger}_{I}(\tau)d\tau\right)^{\bar{n}^{I}}$
$\displaystyle\times$
$\displaystyle\prod_{i=1,2}\frac{1}{n^{c}_{i}!}\left(\int_{0}^{\beta}H_{c,i}(\tau)d\tau\right)^{n_{i}^{c}}\frac{1}{\bar{n}^{c}_{i}!}\left(\int_{0}^{\beta}H^{\dagger}_{c,i}(\tau)d\tau\right)^{\bar{n}_{i}^{c}}$
$\displaystyle\times$
$\displaystyle\frac{1}{n^{b}_{i}!}\left(\int_{0}^{\beta}H_{b,i}(\tau)d\tau\right)^{n_{i}^{b}}\frac{1}{\bar{n}^{b}_{i}!}\left(\int_{0}^{\beta}H^{\dagger}_{b,i}(\tau)d\tau\right)^{\bar{n}_{i}^{b}}]\,,$

where
$p=n^{c}_{1}+\bar{n}^{c}_{1}+n^{c}_{2}+\bar{n}^{c}_{2}+n_{1}^{b}+\bar{n}_{1}^{b}+n^{I}+\bar{n}^{I}$.
Then each term in the expansion is sampled by the Metropolis-Hastings
algorithm Metropolis  _et al._ (1953); Hastings (1970).

The $(-1)^{p}$ factor in the expression could potentially cause a sign
problem. However, it is always positive in this model. Notice that the number
of conduction electrons at each site is a conserved quantity, meaning that
$n_{i}^{c}=\bar{n}_{i}^{c}$. Besides, $S_{1}^{z}$, $S_{2}^{z}$, and number of
bosons are all conserved quantity, such that
$n_{1}^{b}+n^{I}=\bar{n}_{1}^{b}+\bar{n}^{I}$,
$n_{2}^{b}+\bar{n}^{I}=\bar{n}_{2}^{b}+{n}^{I}$, and
$n_{1}^{b}+n_{2}^{b}=\bar{n}_{1}^{b}+\bar{n}_{2}^{b}$. Under these conditions,
$p$ is always even and thus $(-1)^{p}=1$.

##  II Quantum phase transition

Figure 5: Left: Binder cumulant $U_{2}$ (a) and the correlation length divided
by inverse temperature $\xi/\beta$ (b) vs. bosonic coupling $g$ with $I=0.2$
at the labeled inverse temperature $1/\beta$. Both quantities gives a crossing
at critical bosonic coupling $g_{c}=0.665(5)$. Right: Scaling collapse of
$U_{2}$ (c) and $\xi$ (d) with a correlation length exponent $\nu^{-1}=0.34$,
critical coupling $g_{c}=0.66$ (from $U_{2}$) and $\nu^{-1}=0.36$,
$g_{c}=0.67$ (from $\xi/\beta$).

To further establish the presence of a QCP between the Kondo screened phase
and the LM phase, we calculate the Binder cumulant for a 3 components order
parameter,
$U_{2}=5/2-{3{\langle(\bf{m}\cdot\bf{m})^{2}\rangle}}/{2{\langle\bf{m}\cdot\bf{m}\rangle^{2}}}$,
and correlation length $\xi$ of staggered spin correlation function
$\chi_{s}^{\alpha}(\tau)$ Binder (1981); Sandvik (2010). Here ${\bf m}$ is the
staggered magnetization with
$m^{\alpha}=\langle\frac{1}{\beta}\int_{0}^{\beta}\left(S_{1}^{\alpha}(\tau)-S_{2}^{\alpha}(\tau)\right)d\tau\rangle$
and $\chi^{\alpha}_{s}(\tau)=\langle
T_{\tau}(S_{1}^{\alpha}(\tau)-S_{2}^{\alpha}(\tau))(S_{1}^{\alpha}-S_{2}^{\alpha})\rangle$
Cai and Si (2019).

At the transition from Kondo screened to LM phase, both $U_{2}$ and
$\xi/\beta$ should display a crossing at the critical coupling $g=g_{c}$, or
in other words, become independent of $\beta$, the system size. As shown in
Fig. 5 (a)(b), this is indeed what we have observed, where we plot $U_{2}$ and
$\xi/\beta$ versus tuning parameter $g$ for various choices of $\beta$. In
fact, the finite-size scaling hypothesis at a second-order QCP predicts that
they should be described by the following scaling form near the QCP,

$\displaystyle U_{2}(g,\beta)$ $\displaystyle=$
$\displaystyle\tilde{U}_{2}\left(\beta^{1/\nu}(g-g_{c})/g_{c}+A/\beta^{\phi/\nu}\right)\,,$
(5) $\displaystyle\xi(g,\beta)$ $\displaystyle=$
$\displaystyle\beta\tilde{\xi}\left(\beta^{1/\nu}(g-g_{c})/g_{c}+A/\beta^{\phi/\nu}\right)\,.$
(6)

$\tilde{U}_{2}$ and $\tilde{\xi}$ are universal functions, $g_{c}$ the
critical coupling, $\nu$ the correlation length exponent. The term
$A/\beta^{\phi/\nu}$ captures the correction to scaling from the leading
irrelevant operator. By performing scaling collapse analysis, we find the
optimal choices for $g_{c}$ and $\nu$ are $g_{c}=0.66(1)$, $\nu^{-1}=0.34(4)$
(from $U_{2}$) and $g_{c}=0.67(1)$, $\nu^{-1}=0.36(4)$ (from $\xi$). From this
we estimate the actual critical value to be $g_{c}=0.665(5)$ and the
correlation length exponent of the Kondo destruction QCP to be
$\nu^{-1}=0.35(5)$, which are consistent with the results derived from
fidelity susceptibility $\chi_{F}$. The discrepancy is due to the larger
scaling corrections of $\chi_{F}$.

##  III The procedure to calculate the pairing susceptibilities

To measure the pairing susceptibility, we use the formula of the four-point
correlation function shown in Ref. Gull  _et al._ (2007). However, due to the
infinite separation in the SU(2) case, $\langle
T_{\tau}\Delta^{\dagger}_{\uparrow\downarrow}(\tau)\Delta_{\downarrow\uparrow}\rangle$
and $\langle
T_{\tau}\Delta^{\dagger}_{\downarrow\uparrow}(\tau)\Delta_{\uparrow\downarrow}\rangle$
can not be measured via this formula, since they do not conserve the
occupation number of $d_{i\sigma}$. We can then utilize the SU(2) symmetry and
obtain these pairing susceptibilities via: $\langle{\cal
T_{\tau}}\Delta^{\dagger}_{\uparrow\downarrow}(\tau)\Delta_{\downarrow\uparrow}\rangle=\langle{\cal
T_{\tau}}\Delta^{\dagger}_{\uparrow\uparrow}(\tau)\Delta_{\uparrow\uparrow}\rangle-\langle{\cal
T_{\tau}}\Delta^{\dagger}_{\uparrow\downarrow}(\tau)\Delta_{\uparrow\downarrow}\rangle,\langle{\cal
T_{\tau}}\Delta^{\dagger}_{\downarrow\uparrow}(\tau)\Delta_{\uparrow\downarrow}\rangle=\langle{\cal
T_{\tau}}\Delta^{\dagger}_{\downarrow\downarrow}(\tau)\Delta_{\downarrow\downarrow}\rangle-\langle{\cal
T_{\tau}}\Delta^{\dagger}_{\downarrow\uparrow}(\tau)\Delta_{\downarrow\uparrow}\rangle$
Hoshino and Werner (2016).

##  IV Ising anisotropic model

The Hamiltonian of Ising anisotropic BFAM is

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{i,\sigma}\epsilon_{d}d_{i\sigma}^{\dagger}d_{i\sigma}+\sum_{i}Un_{i\uparrow}n_{i\downarrow}+IS^{z}_{1}S^{z}_{2}$
(7) $\displaystyle+$ $\displaystyle\sum_{{\bf k},\sigma}\epsilon_{\bf
k}c_{{\bf k}\sigma}^{\dagger}c_{{\bf k}\sigma}+\sum_{i,{\bf
k},\sigma}\left(Ve^{i{\bf k\cdot r_{i}}}d_{i\sigma}^{\dagger}c_{{\bf
k}\sigma}+h.c.\right)$ $\displaystyle+$
$\displaystyle\sum_{p}\omega_{p}{\phi^{z}_{p}}^{\dagger}\phi^{z}_{p}+g\sum_{p}(S_{1}^{z}-S_{2}^{z})({\phi_{p}^{z}}^{\dagger}+\phi_{p}^{z})\,.$

, where only $z$-component of RKKY interaction and spin-boson coupling are
kept. The effect of conduction electrons is described by two hybridization
functions

$\displaystyle\Gamma(\epsilon)=\pi V^{2}\sum_{\bf
k}\delta(\epsilon-\epsilon_{\bf k})\approx\Gamma_{0}\Theta(D-|\epsilon|)$
$\displaystyle\Gamma_{a}(\epsilon)=\pi V^{2}\sum_{\bf
k}\delta(\epsilon-\epsilon_{\bf k})\cos({\bf k}\cdot({\bf r_{1}}-{\bf
r_{2}}))\approx\gamma\epsilon\Gamma_{0}\Theta(D-|\epsilon|)\,,$ (8)

where $\Gamma(\epsilon)$($\Gamma_{a}(\epsilon)$) describes the effective
hybridization from fermionic bath between two $d$ electrons of the same
(different) sites Zhu and Zhu (2011). By tuning parameter $\gamma$, we are
able to change the separation between two $d$ electrons: $\gamma=0$ denotes
infinite separation $|{\bf r_{1}}-{\bf r_{2}}|=\infty$ and $\gamma=1$
represents $|{\bf r_{1}}-{\bf r_{2}}|=a$ with $a$ the lattice spacing. The
finite separation generates an effective SU(2)-symmetric Heisenberg
interactions of strength $2\gamma^{2}\Gamma^{2}/(\pi^{2}U)$ between two $d$
electrons. Thus, we can estimate the effective spin-spin interactions
$J_{z}S^{z}_{1}S^{z}_{2}+J_{p}(S^{x}_{1}S^{x}_{2}+S^{y}_{1}S^{y}_{2})$ from
fermionic bath, bosonic bath and RKKY interactions:
$J_{z}=I_{z}+2g^{2}+2\gamma^{2}\Gamma^{2}/(\pi^{2}U)$ and
$J_{p}=2\gamma^{2}\Gamma^{2}/(\pi^{2}U)$. In the calculation, we take
$\gamma=0.2$ and quantify the Ising anisotropy near QCP $J_{p}/J_{z}\approx
0.2$. Other parameters are taken to be
$U=-2\epsilon_{d}=0.001,\Gamma=0.25,s=0.8,D=1,\Lambda=1$


# Displaying dark matter constraints from colliders with varying simplified
model parameters

Andreas Albert  Antonio Boveia  Oleg Brandt  Eric Corrigan  Zeynep Demiragli
Caterina Doglioni  Etienne Dreyer  Boyu Gao  Josh Greaves  Ulrich Haisch
Philip Harris  Greg Landsberg  Alexander Moreno  Katherine Pachal  Priscilla
Pani  Federica Piazza  Tim M. P. Tait  David Yu  Felix Yu  Lian-Tao Wang

###### Abstract

The search for dark matter is one of the main science drivers of the particle
and astroparticle physics communities. Determining the nature of dark matter
will require a broad approach, with a range of experiments pursuing different
experimental hypotheses. Within this search program, collider experiments
provide insights on dark matter which are complementary to direct/indirect
detection experiments and to astrophysical evidence. To compare results from a
wide variety of experiments, a common theoretical framework is required. The
ATLAS and CMS experiments have adopted a set of simplified models which
introduce two new particles, a dark matter particle and a mediator, and whose
interaction strengths are set by the couplings of the mediator.

So far, the presentation of LHC and future hadron collider results has focused
on four benchmark scenarios with specific coupling values within these
simplified models. In this work, we describe ways to extend those four
benchmark scenarios to arbitrary couplings, and release the corresponding code
for use in further studies. This will allow for more straightforward
comparison of collider searches to accelerator experiments that are sensitive
to smaller couplings, such as those for the US Community Study on the Future
of Particle Physics (Snowmass 2021) [1], and will give a more complete picture
of the coupling dependence of dark matter collider searches when compared to
direct and indirect detection searches. By using semi-analytical methods to
rescale collider limits, we drastically reduce the computing resources needed
relative to traditional approaches based on the generation of additional
simulated signal samples.

##  1 Introduction

The nature of dark matter is one of the most compelling open questions
confronting physics today. Across the particle and astroparticle physics
communities, a wide range of experiments are attempting to observe and
characterise dark matter (DM) and understand what its relationship to the
Standard Model (SM) may be. The range of possibilities for the nature of DM
and its interactions is vast (see [2] for a review) and as a result, the
experimental scope is similarly enormous and the approaches taken highly
varied. One important experimental direction is the search for DM production
at colliders [3, 4].

To understand how all the different experiments across different areas of
physics complement one another in their collective search for DM, a common set
of benchmarks for comparing results is required. The complementarity between
experimental frontiers is both highly important and difficult to capture.
Indirect detection, direct detection, and collider experiments each have
unique strengths in particular areas of the space of possible DM masses and
models. To fully explore this space, understand what is covered
experimentally, and identify gaps in the broader experimental program, common
frameworks have to be defined in which exclusion limits from current and
future experiments can be visualised together. This contextualisation requires
that the models in which the various experiments interpret their results can
be meaningfully compared.

The ATLAS and CMS collaborations produce a vast array of searches which can be
interpreted in DM contexts, including final states that correspond to
supersymmetric, 2HDM+a, and complicated dark sector models. These specific
interpretations make powerful statements about the theories in question but
are difficult to compare to one another and to results from other fields. To
provide a common framework that would allow for broader interpretation of
ATLAS and CMS DM results, the ATLAS/CMS Dark Matter Forum established a
specific set of simplified models which the two collaborations have used as
the basis for many of their DM searches throughout the LHC Run 2 [5]. In
addition to the simplified models themselves, a smaller set of specific
benchmarks with fixed coupling values was established, allowing for clear
comparisons between and within the LHC experiments [6, 7, 8, 9].

There are contexts, however, in which these benchmark scenarios are not as
flexible as an optimal comparison would require. For example, the HL-LHC and
future colliders should be able to probe smaller coupling values than the LHC,
low-mass experiments often present results for dark photon models that don’t
align clearly with the current LHC DM benchmarks, and LHC limits reframed in
terms of DM-nucleon interaction cross section can vary with coupling choices
in a way which is not fully captured by the existing benchmarks alone.

Historically, reinterpreting ATLAS and CMS analysis limits for a new set of
couplings within the simplified model framework required generating new Monte
Carlo signal simulations and re-running the statistical analysis. This made
testing more than a handful of scenarios prohibitive. But for two signals with
the same physical signature in the detector and differing only in cross
section, the limit set on one signal can be reinterpreted as a limit on the
other by direct rescaling, as long as the ratio between the cross sections of
the two signals is known. Having reliable analytical approximations for the
cross sections of the simplified models relative to each other will then
enable analyses to begin from a single limit in one set of couplings and
rescale to any other, as long as the original result includes sufficient
information.

This paper presents analytical and semi-analytical methods enabling rescaling
for two of the standard DM simplified models used by ATLAS and CMS with
reasonable accuracy and without resorting to generating events. The rescaling
formulas are valid for any hadron collider and are based on leading order
cross sections for the production of the most important final states for these
models. The Python package developed to streamline the application of these
techniques can be found in GitHub here. The present version of the code and
this whitepaper is intended to support Snowmass 2021 [1]. We will document the
full code as part of a future version of this whitepaper that we plan to
publicly release within the LHC Dark Matter Working Group.

To briefly summarise what is not included in this framework: rescaling certain
of the LHC simplified models is not currently supported (see Section 2), limit
rescaling for lepton colliders is not currently supported, and both analysis
acceptances and $k$-factors are treated as constant by the rescaling (see
Section 3). The first two elements can be added in the future by extensions to
the current work if there is sufficient interest from the community. The
third, treatment of acceptances and $k$-factors, must be considered on an
individual basis by users of the framework.

Additionally, this paper uses a similar analytical approach to determine the
minimum allowed coupling values within the simplified models for which the
resulting scenario is compatible with the relic abundance of DM. The results
provided make it possible to compare the non-excluded points in a mass-mass
plane for any set of couplings to the minimum allowed coupling compatible with
the relic density in those areas and draw a conclusion about the relevance of
those regions. They enable a physically realistic interpretation of the
simplified models and illuminate their most interesting regions.

##  2 Models considered

The simplified models considered in this study are detailed extensively in
Refs. [5, 7]. Each simplified model adds two new particles beyond the Standard
Model: a Dirac fermion dark matter particle $\chi$ and an $s$-channel mediator
particle. This mediator can be spin-1, in which case its couplings can be
either vector or axial-vector in nature, or spin-0, in which case its
couplings can be scalar or pseudo-scalar. In this paper, we will refer to
these four different mediator types as different models.

Each model introduces two types of new vertex: a vertex between the mediator
and a pair of DM particles, and a vertex between the mediator and a pair of SM
fermions. There is no interaction between the mediator and SM bosons, and no
interaction between the DM particle and any particle other than the mediator.
The strength of the interactions at each vertex is governed by a coupling $g$.
ATLAS and CMS treat the mediator-SM interactions as having a single universal
coupling strength $g_{q}$ for all mediator-quark interactions and a single
universal coupling strength $g_{l}$ for all mediator-lepton interactions,
including neutrinos. Mediators with smaller masses will be kinematically
forbidden from decaying into heavy quarks or leptons, but the coupling is
still taken to be constant for all fermion generations. For vector and axial-
vector mediators, this leads to universal quark and lepton couplings. For
scalar and pseudoscalar mediators, the fermion couplings are proportional to
the SM Higgs couplings, where the overall relative strengths of the quark and
lepton coupling types are scaled by $g_{q}$ and $g_{l}$. No Higgs mixing is
included in the spin-0 simplified models, and no mixing with the Z-boson is
included in the spin-1 models. These choices mean that the four simplified
models are in fact fully defined by five free parameters: the mediator mass
$M_{\rm{med}}$, the DM particle mass $M_{\chi}$, the coupling strength between
the mediator and DM particles $g_{\chi}$, and the two mediator-SM universal
coupling strengths $g_{q}$ and $g_{l}$. In this paper, we will refer to
different sets of coupling values within the same model as different
scenarios.

The rescaling methods presented here are computed using leading-order
approximations for the cross sections of the key processes. For vector and
axial-vector mediators, the relevant processes are resonant SM particle pair
production, with or without accompanying initial state radiation (ISR), and
invisible decay of a mediator in the presence of ISR leading to a final state
of a single SM particle and significant missing energy. Examples of these
processes, a dijet and a monojet scenario, are illustrated in Fig. 1. Decays
of the mediator to leptons, or the presence of an ISR particle (jet, photon,
etc), also provide visible signatures with significant exclusion power; mono-
photon or other mono-boson analyses are also relevant invisible signatures.
However, in all such cases, the change in cross section due to a variation in
the DM model or scenario studied can be estimated from the core
$qq\rightarrow\text{med}\rightarrow f\bar{f}$ or
$qq\rightarrow\text{med}\rightarrow\chi\bar{\chi}$ process; the effect of an
ISR particle cancels out in the cross section ratio used for rescaling.

av-dijet (100,70)thick i1,i2 o1,o2
wiggly,tension=0.6,label=$V,,A\text{-}V$v1,v2 fermion,tension=0.6o2,v2,o1
fermion,tension=0.6i2,v1,i1 v1,v2 ${g_{q}}$v1 ${g_{q}}$v2 ${\bar{q}}$i1
${q}$i2 ${\bar{q}}$o1 ${q}$o2

av-monojet (100,70)thick i1,i2 o1,o2 isr pisr
wiggly,tension=0.6,label=$V,,A\text{-}V$v1,v2 fermion,tension=0.6o2,v2,o1
fermioni2,visr,v1 plainv1,pvisr,i1 fermion,tension=0v1,i1 v1,v2 ${g_{q}}$v1
${g_{DM}}$v2 ${\bar{q}}$i1 ${q}$i2 ${\bar{\chi}}$o1 ${\chi}$o2
gluon,tension=0visr,isr phantom,tension=0pvisr,pisr ${g}$isr

Figure 1: Representative Feynman diagrams for leading-order processes with a
vector or axial-vector $s$-channel mediator. The visible decay mode (left)
uses a dijet final state as an example. The invisible decay mode (right) uses
a monojet final state as an example.

For scalar and pseudoscalar mediators, the leading order DM or SM production
diagrams are significantly more complex and additional signatures like
$t\bar{t}+E_{T}^{\mathrm{miss}}$ become important. For this reason, although
similar rescaling formulas could be computed for the scalar and pseudoscalar
models, they are beyond the scope of the present paper. Rescaling options
presented here will therefore focus only on the vector and axial-vector
models. Any readers interested in using or helping to develop
scalar/pseudoscalar rescaling are invited to contact the authors. Relic
density calculations are available for all four models.

##  3 Assumptions and caveats

Although the rescaling formulae are based on the leading order cross section
calculations for the visible and invisible final states, recent experimental
results are typically reported at next-to-leading order (NLO). NLO accuracy
can be reasonably approximated in many cases by multiplying an original
exclusion limit calculated at NLO by the LO/LO cross section scale factor
computed using the methods here. The assumption involved is that the
$k$-factors are unchanged between the start and end points of the rescaling.
For example, for the same ($M_{\rm{med}}$, $M_{\chi}$) point in an axial model
with couplings $g_{a,b,c}$ being rescaled to a vector model with couplings
$g_{e,f,g}$, as long as the $k$-factors at these points are the same in the
two models and scenarios, a rescaled NLO limit at ($M_{\rm{med}}$, $M_{\chi}$)
will be NLO also. Similarly, if a visible final state limit defined as a
function of $M_{\rm{med}}$ alone for fixed $M_{\chi}$ is used to extract
limits with various $M_{\chi}^{\prime}$, it is assumed that the $k$-factor at
($M_{\rm{med}}$, $M_{\chi}$) is the same as the $k$-factor at ($M_{\rm{med}}$,
$M_{\chi}^{\prime}$) within the model and scenario being studied. This
assumption should be checked by the user at a few representative points using
an NLO Monte Carlo simulation. If the $k$-factors are found to be
inconsistent, a correction should be derived and applied or the bounds of the
rescaling adjusted to remove the inconsistent points.

A second assumption of the rescaling method is that the effects of a signal on
an analysis result are constant across the start and end points of a
rescaling, less the effect of the cross section. Therefore a limit calculated
for a signal point with a particular acceptance should not be rescaled to act
as a limit on a signal point with very different acceptance, and a resonance
search which is sensitive to the peak shape of a signal should reinterpret
points only to signals with a similar peak shape and location. For
$E_{T}^{\mathrm{miss}}+X$ analyses, the acceptance is not overly sensitive to
model parameters, but much like the $k$-factor assumption, the
$E_{T}^{\mathrm{miss}}$ for a few representative points towards the extremes
of the rescaling should be checked to ensure the acceptance will not change
meaningfully. For resonant final states, the relevant parameter is the
mediator’s intrinsic width. For intrinsic width to mass ratios significantly
smaller than the analysis resolution, varying the couplings will have no
effect on the mediator width so long as the values remain small, but when the
analysis resolution is quite good or the signals become wider, the effects can
be non-negligible. To account for this, the maximum intrinsic width to which a
given limit is valid should be specified by the analyser, and the rescaling
will not be applied for points outside this value. To support analyses like
dilepton where the resolution is excellent, the rescaling code has been
developed to support an input range of observed resonant limits applicable for
different mediator widths, up to the maximum width limit provided. Details of
this can be found in Section 4.2.

##  4 Rescaling resonant final states

For resonant final states, as discussed in Ref. 2, the effects on the cross
section of varying model parameters can be explored by considering the process
$qq\rightarrow\text{med}\rightarrow f\bar{f}$ regardless of the presence of
ISR. The kinematics, acceptance, and amount of background observed in the
final state do depend on the presence of ISR, as well as strongly depending on
the mass of the mediator. These factors enforce that any reinterpretation of
limits must be within a single analysis final state, that is, must not alter
the presence or absence of additional radiation, and cannot rescale a limit
based on one mediator mass to constrain a different mediator mass. This last
point is very clear when one considers how the exponentially decreasing
background in a resonance search must result in a very different number of
observed events corresponding to two different invariant masses.

Resonant analysis limits can be reinterpreted across changing couplings, since
varying couplings can affect only the cross section (the target for the
rescaling) and the intrinsic width but not the invariant mass of an observed
signal. Variations in intrinsic width will affect the detector signature of a
signal point if that intrinsic width of the model is on the order of the
experimental mass resolution. This is the case with high-resolution final
states such as dilepton and can be handled by using analysis limits
parameterised by intrinsic width. This will be discussed further in Section
4.2.

One key point distinguishing visible from invisible final states is that
resonant (visible) limits can also be reinterpreted across varying DM masses
$M_{\chi}$. With no coupling to DM present in the leading order interaction
for these final states, changing $M_{\chi}$ will only affect the result
through variation in the cross section and intrinsic width, as in the case of
varying couplings. This means the analysis limits used as inputs to derive
exclusions in the $M_{\rm{med}}$-$M_{\chi}$ plane can be one-dimensional and
parameterised only as a function of $M_{\rm{med}}$.

The dijet and dilepton analyses typically publish a one-dimensional limit
corresponding to a single model and set of couplings, which can take one of
two standard forms. The first form is an observed 95% confidence level upper
limit on the cross section of new physics $\sigma_{\text{obs}}$ alongside a
theoretical prediction for the cross section with any appropriate acceptance
and efficiency factors $\sigma_{\text{th}}$, both as functions of
$M_{\rm{med}}$, and whose intersection point marks the upper limit on
$M_{\rm{med}}$ for the model and scenario of the theory curve. The second form
is a limit on $g_{q}$ or $g_{l}$ as a function of $M_{\rm{med}}$ with all
other couplings and $M_{\chi}$ fixed. These forms are in reality highly
interconnected: in a dijet final state with all other couplings and $M_{\chi}$
fixed, for a specific $M_{\rm{med}}$, the the theoretical prediction for the
cross section will scale proportionally to $g_{q}^{2}$ and the limit set on
the coupling $g_{q,\text{lim}}$ is just that value for which the theoretical
prediction matches the observed limit. The two limit forms can then be
converted directly from one to the other via

$g_{q,\text{lim}}=g_{q,\text{th}}\sqrt{\frac{\sigma_{\text{obs}}}{\sigma_{\text{th}}}}\,,$
(4.1)

where $g_{q,\text{th}}$ is the value of $g_{q}$ in calculating
$\sigma_{\text{th}}$. The equivalent for a dilepton final state, where the
cross section scales proportionally to $g_{q}g_{l}$ and $g_{q}$ is fixed, is
$g_{l,\text{lim}}=g_{l,\text{th}}(\sigma_{\text{obs}}/\sigma_{\text{th}})$.

Regardless of starting point, the goal is to produce an estimate of the
exclusion depth in a two-dimensional plane defined by $M_{\rm{med}}$ and
$M_{\chi}$. We define exclusion depth $d_{\text{ex}}$ to be the ratio of the
experimental cross section upper limit divided by the theoretical cross
section for any point in the plane. The observed limit in the mass-mass plane
is set by the contour $d_{\text{ex}}=1$. Points with $d_{\text{ex}}<1$ are
excluded while those with $d_{\text{ex}}>1$ are not. With this information,
rescaling to any model and scenario is as simple as multiplying the exclusion
depth by the ratio of cross sections in the source and target scenario.

This requires a good approximation of the cross section for the processes in
question. Fortunately, these resonant processes will be on-shell everywhere.
For any $M_{\rm{med}}$ of interest, in these visible final states, there will
always be quarks or leptons with $m_{f}<2\,M_{\rm{med}}$ into which the
mediator can decay. Therefore this simple on-shell approximation can be
applied for all points:

$\tilde{\sigma}\sim\frac{\Gamma_{i}\Gamma_{f}}{\Gamma_{\text{tot}}}\,,$ (4.2)

where $\Gamma_{i}$ and $\Gamma_{f}$ are the partial widths of the initial and
final states respectively and $\Gamma_{\text{tot}}$ is the total intrinsic
width of the mediator. This quantity is easy to compute since formulas for the
partial widths are well-known in these models [7].

The first step is to rescale the one-dimensional $M_{\chi}$ limit into a full
grid of exclusion depths in ($M_{\rm{med}}$, $M_{\chi}$). In the initial one-
dimensional limit on the cross sections, in which
$d_{\text{ex}}=\sigma_{\text{obs}}/\sigma_{\text{th}}$, $\sigma_{\text{th}}$
corresponds to a specific set of couplings $g_{q}$, $g_{\chi}$, $g_{l}$ as
well as a fixed $M_{\chi}$ and the $M_{\rm{med}}$ value of each point on the
x-axis. To get the exclusion depth in a target point with
$M^{\prime}_{\rm{med}}$ = $M_{\rm{med}}$ but varied $g_{q}^{\prime}$,
$g_{\chi}^{\prime}$, $g_{l}^{\prime}$ and $M^{\prime}_{\chi}$, the original
exclusion depth can simply be multiplied by the ratio of the cross section
approximations for the two scenarios:

$d_{\text{ex}}^{\prime}=d_{\text{ex}}(\frac{\tilde{\sigma}}{\tilde{\sigma}^{\prime}})=\frac{\sigma_{\text{obs}}}{\sigma_{\text{th}}}(\frac{\tilde{\sigma}}{\tilde{\sigma}^{\prime}})\,.$
(4.3)

Note that $\sigma_{\text{th}}$ and $\tilde{\sigma}$ are not the same:
$\sigma_{\text{th}}$ is the full cross section including any analysis-relevant
factors while $\tilde{\sigma}$ is an approximation missing various
multiplicative factors and which is only relevant when compared to the same
approximation at other model and scenario points.

When instead beginning from a limit on $g_{q}$ as a function of
$M_{\rm{med}}$, the exclusion depth in the initial plot is defined by the
relation in Eq. 4.1: $d_{\text{ex}}=(g_{q,\text{lim}}/g_{q,th})^{2}$.
Combining this with Eq. 4.3,

$d_{\text{ex}}^{\prime}=(\frac{g_{q,\text{lim}}}{g_{q,th}})^{2}\frac{\tilde{\sigma}}{\tilde{\sigma}^{\prime}}\,,$
(4.4)

where $g_{q,th}$ and $\tilde{\sigma}$ are defined for a theory curve which,
combined with a matching observed limit curve, would result in the input
$g_{q}$ limit. These are not defined as there is no input cross section curve,
but their effect in fact cancels out: $\tilde{\sigma}\propto g_{q}^{2}$ for
dijets and so selecting any arbitrary $g_{q}$ and evaluating $\tilde{\sigma}$
at the same will produce the same result. For convenience, we choose unity.
Then using a $g_{q}$ limit at a hadron collider:

$d_{\text{ex}}^{\prime}=g^{2}_{q,\text{lim}}(\frac{\tilde{\sigma}(g_{q}=1)}{\tilde{\sigma}^{\prime}})\,,$
(4.5)

all of which quantities can be calculated for a grid of ($M_{\rm{med}}$,
$M_{\chi}$) values from the one input curve. The equivalent exclusion depth
extracted from a $g_{l}$ limit at a hadron collider with fixed $g_{q}$ is:

$d_{\text{ex}}^{\prime}=g_{l,\text{lim}}(\frac{\tilde{\sigma}(g_{l}=1)}{\tilde{\sigma}^{\prime}})\,.$
(4.6)

Once one of the above formulas has been used to calculate
$d_{\text{ex}}^{\prime}$ at each point of interest in the mass-mass plane for
a single model and coupling scenario, the next step is to rescale to any other
scenario or model of interest. This is achieved through multiplying the
exclusion depth by a simple factor of
$\tilde{\sigma}/\tilde{\sigma}^{\prime}$, where $\tilde{\sigma}$ is computed
using Eq. 4.2 for the set of couplings of the scenario from which to begin and
$\tilde{\sigma}^{\prime}$ is in the target scenario. Since the expressions for
the partial widths include the factors that differ between vector and axial-
vector couplings, these differences are included in the calculation and the
simple $\tilde{\sigma}/\tilde{\sigma}^{\prime}$ rescaling can be used to
convert between the two models. The dijet and dilepton resonant final states
and their variants do not provide strong constraints on scalar and
pseudoscalar models, where other final states are dominant. Conversion between
vector and axial-vector as well as across couplings within each model
therefore covers the full space of interesting rescaling possibilities for
these two models.

The rescaling method discussed here for resonant final states is based on a
method used by CMS collaborators to present the dijet and dilepton limits in
the CMS summary plots since the beginning of Run 2. The technique has been
since been adopted by ATLAS for the presentation of their limit summaries.
This paper makes one small extension on the well-established rescaling
technique by enabling limits to be extracted from $g_{q}$ limits as a function
of mass regardless of whether $M_{\chi}$ is taken to be decoupled or light.

###  4.1 Dijet and dijet+X

The simplest example of resonant limit rescaling is for an analysis when the
intrinsic mediator width-to-mass ratio is smaller than or comparable to the
experimental mass resolution. This is common for dijet analyses where the mass
resolution is on the order of 4% at $M_{\rm{med}}=1$ TeV and 3% at
$M_{\rm{med}}=2$ TeV [10, 11]. The rescaling methods derived in Section 4 are
illustrated using the results of the CMS dijet analysis in 36 $\rm{fb}^{-1}$
of data [12]. Since CMS has been using these resonance search rescaling
methods for some time, the mass-mass exclusion limits presented in that paper
are already derived using the procedures given here and beginning from the
limit on $g_{q}$ as a function of $M_{\rm{med}}$ with no coupling to dark
matter. To validate the equivalency of rescaling from cross-section or $g_{q}$
limits with or without decoupled dark matter, Figure 2 shows that the same
mass-mass limits can be obtained beginning from the cross-section limit and
theory curve corresponding to $g_{q}=0.25$, $g_{\chi}=1.0$, $M_{\chi}$ $=1$
GeV instead. The observed limit on $\sigma\times B\times A$ is taken from the
public data while the theory curve and is visually extracted from the limit
plot.

(a)

(b)

Figure 2: Exclusion contours in the ($M_{\rm{med}}$, $M_{\chi}$) plane for
axial-vector mediators (2(a)) and vector mediators (2(b)) with couplings of
$g_{q}=0.25,g_{\chi}=1.0,g_{l}=0.0$ using the CMS dijet analysis results. The
shape of this exclusion surface is shown as a blue gradient. Each point at
which the exclusion depth was calculated is plotted as a dot: excluded points
are plotted in white, while points which are not excluded are plotted in red.
The solid red curves correspond to the published CMS results in these planes
computed using the rescaling methods discussed here and beginning from the
limit on $g_{q}$ as a function of mass for a mediator with no coupling to dark
matter. The solid white curves correspond to the result derived using the
rescaling methods discussed here and beginning from the observed limit on
$\sigma\times B\times A$ as a function of $M_{\rm{med}}$ and the corresponding
theoretical prediction for a mediator with coupling $g_{\chi}=1$ and a dark
matter particle with mass $M_{\chi}=1$ GeV. The level of agreement between the
two contours is a measure of the similarity in performance between the
rescaling procedures beginning from these two different initial limits.

Determining the largest intrinsic width for which the input observed limit is
applicable is the responsibility of each user of the provided software. The
software can then take this largest acceptable width as an input and will only
return rescaled limits when the point being tested falls within this region of
validity. This is intended to avoid accidental reporting of a limit in an
unjustified context. For cases where a single observed limit curve cannot
apply to all the cases of interest, a modification of this approach is
applied, as described in the following section.

###  4.2 Dilepton and other high-resolution resonances

In many resonance search situations, the experimental mass resolution is very
good and can be significantly smaller than the intrinsic width of these
mediators, leading to observed limits that would vary significantly for
different mediators with the same $M_{\rm{med}}$. One example is a di-electron
final state, in which the mass resolution is about 1% for all mediator masses
over $\sim 800$ GeV [13]. In this case, multiple observed limits can be
provided corresponding to different intrinsic mediator widths. During the
rescaling process, at each point to be studied, the observed limit used to
compute $d_{\text{ex}}$ will be found by interpolating between the input
observed limits at the two nearest widths to that of the studied point. The
implementation provided with this paper uses a linear interpolation between
each neighbouring pair of observed limit inputs. All points with intrinsic
widths smaller than that of the narrowest provided limit curve will take their
observed values from that narrowest limit curve with no interpolation or
extrapolation. Points with a larger intrinsic width than that of the widest
provided limit curve will be considered to fall outside the bounds of analysis
applicability and will not return rescaled limits.

This method can be applied to dijet final states as well, or to any resonant
final state where it is desirable to eliminate the assumption that the
observable signal shape and acceptance are unchanged as a function of DM mass
and coupling values. By providing a range of limits corresponding to various
intrinsic width-to-mass ratios for a single mediator mass point, the most
appropriate one can be chosen for each coupling and DM mass.

##  5 Rescaling $E_{T}^{\mathrm{miss}}+X$ final states

For $E_{T}^{\mathrm{miss}}+X$ signatures, the relevant final state involves
decay of a mediator to dark matter particles, and therefore the off-shell case
has to be handled appropriately. The approximations used for visible resonant
final states in Section 4 require a well-defined decay width to the final
state $\Gamma_{f}$, which goes to zero at $M_{\rm{med}}=2M_{\chi}$ for an
invisible final state. The approximation in fact loses validity before this
diagonal is fully reached, as the transition across the on-shell to off-shell
boundary is smooth. A method for re-scaling $E_{T}^{\mathrm{miss}}+X$
signatures has therefore been developed with the goal of handling this
transition smoothly such that it is applicable in all regimes.

Previously, ATLAS and CMS $E_{T}^{\mathrm{miss}}+X$ analyses have generated a
grid of signal points in one of the four benchmark scenarios at full
reconstruction level and used these to determine the region of
($M_{\rm{med}}$, $M_{\chi}$) space excluded in that scenario [14, 15]. Three
additional grids of signal points are then generated at LO and particle level
in order to obtain cross section estimates for each point in the additional
three benchmark scenarios. The ratio of cross sections between the points in
the original scenario and the target scenarios is used to scale the limits,
resulting in a new estimate of the excluded and non-excluded points within
each target scenario. The goal of these studies is to use the same starting
information, a full grid of results in the ($M_{\rm{med}}$, $M_{\chi}$) plane
for one coupling scenario, but rescale to another target scenario quickly and
mathematically instead of through signal generation. Note that unlike in the
resonant case where a one-dimensional limit in $M_{\rm{med}}$ is sufficient, a
two-dimensional grid is required here. Because dark matter particles appear in
the final state for $E_{T}^{\mathrm{miss}}+X$ signatures, both $M_{\chi}$ and
$M_{\rm{med}}$ affect the $E_{T}^{\mathrm{miss}}$ distribution and analysis
acceptance, and it is therefore not possible to rescale a limit across either
of these dimensions.

The semi-analytical rescaling method developed here for
$E_{T}^{\mathrm{miss}}+X$ signatures has two separate components. One can be
used to rescale an initial scenario to another set of couplings within the
same overall model, while the other should be used when a scenario in one
model is rescaled to a target in a different model (e.g. vector mediator to
axial-vector mediator).

In the ATLAS/CMS Dark Matter Forum whitepaper [5] defining the simplified
models studied here, it is specified that the cross section scaling can be
estimated using the integral of the Breit-Wigner propagator for the mediator.
Several approximations of this integral are given, corresponding to the
different regimes of on-shell mediators, off-shell mediators, and effective
field theories. In order to smoothly handle the on-shell to off-shell
transition, we instead use the full integral of the propagator over the
permitted phase space $s\geq 4M_{\chi}^{2}$:

$\displaystyle\mathcal{I}_{\text{prop}}=$ $\displaystyle\
g_{q}^{2}g_{\chi}^{2}\int_{4M_{\chi}^{2}}^{\infty}\frac{ds}{(s-M_{\rm{med}}^{2})^{2}+M_{\rm{med}}^{2}\Gamma^{2}}$
(5.1) $\displaystyle=$ $\displaystyle\ \frac{g_{q}^{2}g_{\chi}^{2}}{\Gamma
M_{\rm{med}}}\left(\frac{\pi}{2}+\arctan{\left(\frac{M_{\rm{med}}^{2}-4M_{\chi}^{2}}{\Gamma
M_{\rm{med}}}\right)}\right)\,.$ (5.2)

Coupling dependence comes not only from the explicit $g_{q}^{2}g_{\chi}^{2}$
factor but is also embedded in the value of $\Gamma$ (see width definitions in
Ref. [7]). $\mathcal{I}_{\text{prop}}$ is easily analytically calculable and
provides a robust estimate of the cross section when comparing two different
coupling scenarios, so long as both assume the same mediator type. A starting
grid of exclusion depths can then be easily rescaled to a target scenario of
the same model by multiplying each $d_{\text{ex}}$ value by the ratio of
$\mathcal{I}_{\text{prop}}$ in the initial and target scenarios at that point.

The $\mathcal{I}_{\text{prop}}$-based rescaling procedure is illustrated in
Figure 3 using the observed limits from the ATLAS
$E_{T}^{\mathrm{miss}}+\gamma$ analysis in 36 $\rm{fb}^{-1}$ of data [16].
Circles show the locations of the signal points used. Each signal point has a
$z$-axis value of theoretical cross section divided by excluded cross section;
white points are those with $z<1$ and are excluded while red points have $z>1$
and are not excluded. The red curves show the exclusion contours reported by
the published analysis in each plane. The blue shades are a linear
interpolation between the points used to illustrate the cross section
exclusion surface. In Figure 3(a) no reinterpretation is done: the values at
each point are taken directly from the paper and overlaid with the
corresponding exclusion curve. In Figure 3(b) the values at each point are
calculated using the $\mathcal{I}_{\text{prop}}$ scaling procedure starting
from the values in Figure 3(a). The white and red colour coding of points is
based on these rescaled values, as is the linear interpolation in blue. The
good agreement between the contour curve taken from the paper and the points
marked as excluded by the rescaling method serves as a validation of the
rescaling.

(a)

(b)

Figure 3: Original cross section limits with couplings
$g_{q}=0.25,~{}g_{\chi}=1.0,~{}g_{l}=0.0$ (3(a)) and rescaled cross section
limits with couplings $g_{q}=0.1,g_{\chi}=1.0,g_{l}=0.1$ (3(b)) for an axial-
vector mediator model using the ATLAS $E_{T}^{\mathrm{miss}}+\gamma$ analysis
results [16]. Limits are presented as a function of mediator mass and dark
matter mass for fixed coupling values. At each signal point, the value of
$d_{\text{ex}}=\sigma_{\text{theory}}/\sigma_{\text{excluded}}$ is calculated:
points where $d_{\text{ex}}<1$ are excluded and plotted in white, while points
where $d_{\text{ex}}>1$ are not excluded and plotted in red. A linear
interpolation between the $d_{\text{ex}}$ values at the points illustrates the
shape of this exclusion surface and is shown as a blue gradient. The solid red
curves are taken from the published $E_{T}^{\mathrm{miss}}+\gamma$ results and
show the exclusion contour $d_{\text{ex}}=1$. The level of agreement between
the excluded points and the published contour in (3(b)) is a measure of the
similarity in performance between the $\mathcal{I}_{\text{prop}}$ scaling
procedure and the Monte Carlo based method used in the ATLAS analysis.

Rescaling cross sections using a ratio based on $\mathcal{I}_{\text{prop}}$ is
no longer sufficient when the target scenario uses a different mediator model
than the initial scenario. Additional mass-dependent terms in the cross
section vary between models, so while they cancel within a single model and
can be ignored in favour of the simple propagator-based expression, they must
be correctly accounted for when rescaling from one model to another. We
therefore calculate a more complete cross section whose ratio can be used to
perform cross-model rescaling.

As discussed in Section 2, the LO signal contribution to each
$E_{T}^{\mathrm{miss}}+X$ analysis can be calculated from the
$qq\rightarrow\text{med}\rightarrow\chi\bar{\chi}$ process independent of
additional ISR radiation. Therefore, to calculate the LO scale factor
translating between two scenarios with different mediator models, the cross
section for that single diagram in both models is sufficient. The parton-level
cross sections as functions of the interaction scale $S$ are given by the
following equations, discounting scale factors which will cancel in the ratio:

$\sigma_{V}(S)\propto\frac{g_{q}^{2}g_{\chi}^{2}(S+2M_{\chi}^{2})\sqrt{S-4M_{\chi}^{2}}}{(\Gamma^{2}M_{\rm{med}}^{2}+(M_{\rm{med}}^{2}-S)^{2})\sqrt{S}}$
(5.3)

and

$\sigma_{AV}(S)\propto\frac{g_{q}^{2}g_{\chi}^{2}(S-4M_{\chi}^{2})^{3/2}}{(\Gamma^{2}M_{\rm{med}}^{2}+(M_{\rm{med}}^{2}-S)^{2})\sqrt{S}}\,,$
(5.4)

where $\Gamma$ is the total width of the mediator (non-zero even when decays
to DM are off-shell).

To calculate cross sections with sufficient accuracy for the rescaling
procedure, a hadron-level quantity must be determined accounting for parton
distribution functions (PDFs). The full cross section is thus calculated semi-
analytically by performing a two dimensional integral over the longitudinal
momentum fractions $x_{1},x_{2}$ of the interacting partons within the allowed
range. Where $S$ is the squared centre-of-mass energy (i.e. $\sqrt{S}=13$
TeV), $\hat{s}=x_{1}x_{2}S$ is the scale of the interaction, and $f^{a}$ and
$f^{b}$ are the PDFs of the incoming partons, the total cross section is

$\sigma_{V/AV}^{\text{tot}}=\int_{4M_{\chi}^{2}}^{\text{S}}\int_{4M_{\chi}^{2}}^{x_{2}S}f^{a}(x_{1},\hat{s})f^{b}(x_{2},\hat{s})\sigma_{V/AV}(\hat{s})dx_{1}dx_{2}\,.$
(5.5)

This integral is performed numerically in the two models and the ratio at each
point is taken as the scale factor to convert between them. The implementation
uses LHAPDF and so provides a wide range of PDF sets and flavour schemes,
although results are found to be fairly independent of the selection [17].

A demonstration of the $\sigma_{V/AV}^{\text{tot}}$-based scaling procedure is
given in Figure 4. The $z$-value at each point in Figure 4(b) is found using
the ratio of the vector and axial-vector cross sections
$\sigma_{V}^{\text{tot}}$ and $\sigma_{AV}^{\text{tot}}$ applied to the axial-
vector limits from Figure 3(a). The success of this rescaling method is again
illustrated by the agreement between the excluded points (white) determined by
the rescaling and the solid curve obtained from the analysis paper. Figure
4(b) is created by further applying a $\mathcal{I}_{\text{prop}}$ rescaling to
the results in Figure 4(a). Even after this iterative application of
$\sigma_{V/AV}^{\text{tot}}$ and $\mathcal{I}_{\text{prop}}$ rescaling, the
agreement between the points and the published results is fair.

(a)

(b)

Figure 4: Rescaled cross section limits with couplings
$g_{q}=0.25,g_{\chi}=1.0,g_{l}=0.0$ (4(a)) and
$g_{q}=0.1,g_{\chi}=1.0,g_{l}=0.01$ (4(b)) for a vector mediator model using
the ATLAS $E_{T}^{\mathrm{miss}}+\gamma$ analysis results. Limits are
presented as a function of mediator mass and dark matter mass for fixed
coupling values. At each signal point, the value of
$d_{\text{ex}}=\sigma_{\text{theory}}/\sigma_{\text{excluded}}$ is calculated:
points where $d_{\text{ex}}<1$ are excluded and plotted in white, while points
where $d_{\text{ex}}>1$ are not excluded and plotted in red. A linear
interpolation between the $d_{\text{ex}}$ values at the points illustrates the
shape of this exclusion surface and is shown as a blue gradient. The solid red
curves are taken from the published $E_{T}^{\mathrm{miss}}+\gamma$ results and
show the exclusion contour $d_{\text{ex}}=1$. The level of agreement between
the excluded points and the published contour is a measure of the performance
of the $\sigma_{V/AV}^{\text{tot}}$ scaling procedure in (4(a)) and of the
combination of both procedures in (4(b)).

The full recommended procedure for rescaling $E_{T}^{\mathrm{miss}}+X$
analysis results is therefore to use the $\sigma_{V/AV}^{\text{tot}}$-based
rescaling procedure to convert from the original model with one type of
mediator into a single scenario with the desired new mediator type (e.g.
vector mediator to axial-vector mediator). The
$\mathcal{I}_{\text{prop}}$-based rescaling procedure can then be used to
convert to other coupling scenarios with the same mediator type.

##  6 Combined examples of coupling-scaled exclusion plots

An example of the opportunities which rescaling across couplings presents is
shown in Figure 5. The limits from the CMS dijet and ATLAS monophoton analyses
studied above are combined into single exclusion contours for a fixed set of
couplings and a range of couplings are overlaid. For this figure, $g_{\chi}$
and $g_{l}$ are held constant while $g_{q}$ is gradually reduced. Both dijet
and monophoton exclusion limits become weaker as $g_{q}$ decreases. Such a
detailed study of the effect of couplings on exclusion limits is possible
because of the drastically reduced computing time involved in determining each
rescaling.

(a)

(b)

Figure 5: Exclusion limits of the dijet and monophoton analyses combined,
corresponding to 36 fb-1 of LHC data, when the coupling to quarks is varied.
The dark matter coupling is held fixed at $g_{\chi}=1.0$ and there is no
coupling between the mediators and leptons. The shaded region is excluded at
the corresponding $g_{q}$ value. Axial-vector (5(a)) and vector (5(b))
mediators show similar patterns, though with a more pronounced cutoff along
the $M_{\rm{med}}=2M_{\chi}$ line for axial-vector mediators. All exclusion
regions shrink as $g_{q}$ decreases.

Similar plots can be made varying $g_{\chi}$ or $g_{l}$ instead. Decreasing
$g_{l}$ greatly reduces the power of dilepton exclusions while slightly
increasing the power of dijet exclusions, since the dijet signatures then
takes a higher fraction of the branching ratio; it has a smaller but nonzero
impact on $E_{T}^{\mathrm{miss}}+X$ exclusion curves. Decreasing $g_{\chi}$
increases the power of visible final state exclusions and decreases that of
$E_{T}^{\mathrm{miss}}+X$ analyses.

##  7 Relic densities in DM simplified models

Calculation of the dark matter relic density involves computing the cross
section for dark matter annihilation to standard model particles. In this
scenario, two production modes tend to dominate, single mediator production
and double mediator production. Single mediator production occurs when the
dark matter annihilates with itself to create a single mediator that then
decays to standard model particles. This production mode is the inverse of the
production mode probed when searching for $E_{T}^{\mathrm{miss}}+X$ final
states, and is resonantly enhanced along the $2M_{\chi}=M_{\rm{med}}$ line.
The second production mode arises from double mediator production. Diagrams
for both of these production modes are shown in Figure 6 [18].

Figure 6: Feynman diagrams of DM $s$-channel annihilation to quarks (left) and
$t$-channel annihilation to a pair of mediators subsequently decaying to
quarks (right). The exchanged $\Phi$ particle(s) can be either (pseudo-)scalar
or (axial-)vector mediator(s).

When computing the relic density, we can approximate broad features of the
calculation by understanding the form of the leading order production cross
section. The leading order cross section can be written as

$\sigma_{{\rm ann},s}^{V}\cdot v=\sum_{q}\frac{N_{c}^{q}\hskip
0.35565ptg_{\chi}^{2}\hskip 0.7113ptg_{q}^{2}\hskip
1.42262pt\beta_{q}}{2\pi}\,\frac{2M_{\chi}^{2}+m_{q}^{2}}{\left(M_{\rm{med}}^{2}-4M_{\chi}^{2}\right)^{2}+M_{\rm{med}}^{2}\Gamma_{\rm
med}^{2}}\,,$ (7.1)

for the single mediator vector production, and it can be written as

$\sigma_{{\rm ann},t}^{V}\cdot v=\frac{g_{\chi}^{4}\hskip 1.42262pt\beta_{\rm
med}}{4\pi}\,\frac{M_{\chi}^{2}-M_{\rm{med}}^{2}}{\left(M_{\rm{med}}^{2}-2M_{\chi}^{2}\right)^{2}}\,,$
(7.2)

for the double mediator production.

When we consider single mediator production, if we take $\sigma_{{\rm
ann},t}^{V}\cdot v$ to be proportional to the thermalized cross section
obtained from the full relic density calculation, and we fix both dark matter
mass $M_{\chi}$ and mediator mass $M_{\rm{med}}$, we find that both equations
above can be written in the form of a quadratic equation. Namely, for some
constants $A,B,C,D,E,F$, we can write

$\displaystyle
Ag_{q}^{2}g_{\chi}^{2}+Bg_{q}^{2}+Cg_{\chi}^{2}+Dg_{q}^{4}+Eg_{\chi}^{4}+F=0$
(7.3)

This also follows from the fact that the widths for SM particles $\Gamma_{\rm
SM}\propto g_{q}^{2}$ while the width for dark matter $\Gamma_{\chi}\propto
g_{\chi}^{2}$. From this form, we observe that if $g_{\chi}$ is a constant,
then this equation immediately reduces to a quadratic equation in $g_{q}^{2}$,
and likewise for $g_{\chi}^{2}$ when $g_{q}$ is a constant. Furthermore, the
quartic terms $D$ and $E$ are also positive. As a consequence, when all
couplings are fixed with the exception of either $g_{q}$ or $g_{\chi}$ there
is a well defined coupling minimum. For couplings smaller than the minimum the
relic is larger, and the same is true for couplings larger than the maximum.
By scanning the couplings, we can find both the minimum and maximum allowed
coupling that would satisfy the dark matter relic density.

Figure 7 shows the result of a procedure whereby we scan the standard model
coupling, compute the relic density, and find the minimum allowed coupling
that would not overproduce dark matter. We present this plot for each of the
simplified models used within the LHC Dark Matter Working Group. We find that
for large regions of phase space, in particular light dark matter with a heavy
mediator, there is no solution that will produce an underabundance of dark
matter. Additionally, we find that couplings of order $g_{q}=1$ are required
to satisfy the dark matter relic density for spin-0 scalar and pseudo-scalar
mediators, and couplings of order $g_{q}=0.1$ or larger are needed to explain
spin-1 vector and axial-vector mediators.

The bounds on these models are clearly substantial, and help to serve as a
guide. We would like to stress that as models become more complicated these
bounds can change substantially. However, this provides a clear set of
benchmarks that can viewed as a target for both invisible, for light dark
matter, and visible, for heavy dark matter, searches.

(a)

(b)

(c)

(d)

Figure 7:  Minimum allowed coupling for the scalar (7(a)), pseudo-
scalar(7(b)), vector(7(c)), and axial-vector(7(d)) mediators. The minimum
coupling is computed on the z-axis. For the scalar and pseudo-scalar
mediators, this coupling is treated as a correction to the Yukawa coupling.

##  8 Conclusion

The simplified dark matter models defined by the LHC Dark Matter Working Group
provide interpretations covering a wide range of signatures and allowing
multiple experiments to frame results in a single set of benchmarks. However,
the LHC experiments have historically been restricted in the range of these
dark matter scenarios they can probe with resonant and
$E_{T}^{\mathrm{miss}}+X$ analyses due to the processing time required to
generate Monte Carlo signal points in each scenario of interest. The methods
described here for analytically rescaling between vector and axial-vector
mediators as well as between different sets of couplings in the framework of
the simplified DM models make this process significantly quicker by avoiding
the need for computationally-expensive simulation. The resonance analysis
rescaling technique has been used in published CMS results for some time, but
is publicly documented here for the first time by its developers. The
$E_{T}^{\mathrm{miss}}+X$ rescaling technique is new and, when combined with
the resonance rescaling, makes it practical and time-efficient to reinterpret
all major limits on the vector and axial-vector mediators, allowing
experiments to create full limit summary plots in any scenario desired.

This paper acts as a user’s guide to understanding and applying analytical
rescaling methods to their own analysis limits. Certain caveats relating to
the consistency of acceptances and $k$-factors across the rescaled scenarios
have been described and recommendations outlined. The present version of this
paper is intended to support the Snowmass 2021 process [1].In a planned second
version of this paper, a Python code package to handle the rescaling will be
included. A preliminary version of this code is already available on Git, or
through request to the authors, and has been passed to members of the ATLAS
and CMS collaborations. Although no rescaling formulae have so far been
developed for lepton colliders or for scalar and pseudoscalar mediators, these
are open points of investigation and may be added to the software in future if
there is evident user interest.

This paper also provides a guide to the relic density implications of the LHC
simplified models. The smallest possible couplings for which dark matter would
not be overproduced are presented for all four simplified models. While these
models are generally used more for their simplicity and ability to be directly
compared across analyses and experiments more than for any probability as
realistic dark matter models, their consequences for dark matter abundance are
valuable to understand and provide an important additional axis of information
for all users of the models.

## Acknowledgments

Research by C. Doglioni is part of projects that have received funding from
the European Research Council under the European Union’s Horizon 2020 research
and innovation program (grant agreement 679305 and 101002463) and from the
Swedish Research Council. K. Pachal’s research is supported by TRIUMF, which
receives federal funding via a contribution agreement with the National
Research Council (NRC) of Canada. Research for P. Harris is supported by a
Department of Energy Early Career Award. Research by A. Boveia is supported by
the US Department of Energy (Grant DE-SC0011726).

## References

  * [1] US Community Study on the Future of Particle Physics (Snowmass 2021), 2022\.  https://snowmass21.org.
  * [2] J. L. Feng, Dark Matter Candidates from Particle Physics and Methods of Detection, Annual Review of Astronomy and Astrophysics 48 (2010), no. 1 495–545, [arXiv:1003.0904].
  * [3] F. Kahlhoefer, Review of LHC Dark Matter Searches, Int. J. Mod. Phys. A32 (2017), no. 13 1730006, [arXiv:1702.02430].
  * [4] A. Boveia and C. Doglioni, Dark Matter Searches at Colliders, Annual Review of Nuclear and Particle Science 68 (2018), no. 1 429–459, [arXiv:1810.12238].
  * [5] D. Abercrombie et al., Dark Matter benchmark models for early LHC Run-2 Searches: Report of the ATLAS/CMS Dark Matter Forum, Physics of the Dark Universe 27 (2020) 100371, [arXiv:1507.00966].
  * [6] A. Boveia et al., Recommendations on presenting LHC searches for missing transverse energy signals using simplified s-channel models of dark matter, Physics of the Dark Universe 27 (2020) 100365, [arXiv:1603.04156].
  * [7] A. Albert et al., Recommendations of the LHC Dark Matter Working Group: Comparing LHC searches for dark matter mediators in visible and invisible decay channels and calculations of the thermal relic density, Physics of the Dark Universe 26 (2019) 100377, [arXiv:1703.05703].
  * [8] ATLAS Collaboration, Dark matter summary plots: update July 2020, Tech. Rep. ATL-PHYS-PUB-2020-021, CERN, Geneva, Jul, 2020.  http://cds.cern.ch/record/2725266.
  * [9] CMS Collaboration, “Dark Matter Searches from CMS.”  http://cms-results.web.cern.ch/cms-results/public-results/publications/EXO/DM.html.
  * [10] ATLAS Collaboration, Search for low-mass resonances decaying into two jets and produced in association with a photon using $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector, Phys. Lett. B 795 (2019) 56–75, [arXiv:1901.10917].
  * [11] ATLAS Collaboration, Search for new resonances in mass distributions of jet pairs using 139 fb-1 of $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector, JHEP 03 (2020) 145, [arXiv:1910.08447].
  * [12] CMS Collaboration, Search for narrow and broad dijet resonances in proton-proton collisions at $\sqrt{s}=13$ TeV and constraints on dark matter mediators and other new particles, JHEP 08 (2018) 130, [arXiv:1806.00843].
  * [13] ATLAS Collaboration, Search for high-mass dilepton resonances using 139 fb-1 of $pp$ collision data collected at $\sqrt{s}=$13 TeV with the ATLAS detector, Phys. Lett. B 796 (2019) 68–87, [arXiv:1903.06248].
  * [14] ATLAS Collaboration, Search for dark matter and other new phenomena in events with an energetic jet and large missing transverse momentum using the ATLAS detector, JHEP 01 (2018) 126, [arXiv:1711.03301].
  * [15] CMS Collaboration, Search for dark matter produced with an energetic jet or a hadronically decaying W or Z boson at $\sqrt{s}=13$ TeV, JHEP 07 (2017) 014, [arXiv:1703.01651].
  * [16] A. Collaboration, Search for dark matter at $\sqrt{s}=13$ TeV in final states containing an energetic photon and large missing transverse momentum with the ATLAS detector, Eur. Phys. J. C 77 (2017), no. 6 393, [arXiv:1704.03848].
  * [17] A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr, and G. Watt, LHAPDF6: parton density access in the LHC precision era, Eur. Phys. J. C 75 (2015) 132, [arXiv:1412.7420].
  * [18] A. Albert et al., Recommendations of the LHC Dark Matter Working Group: Comparing LHC searches for dark matter mediators in visible and invisible decay channels and calculations of the thermal relic density, Phys. Dark Univ. 26 (2019) 100377, [arXiv:1703.05703].