11institutetext: Computer Science and Artificial Intelligence Laboratory, MIT,
Cambridge, MA, USA
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# A Geometry-Sensitive Quorum Sensing Algorithm for the Best-of-N Site
Selection Problem

Grace Cai  0000-0001-5929-4347 Nancy Lynch  0000-0003-3045-265X

###### Abstract

The house hunting behavior of the Temnothorax albipennis ant allows the colony
to explore several nest choices and agree on the best one. Their behavior
serves as the basis for many bio-inspired swarm models to solve the same
problem. However, many of the existing site selection models in both insect
colony and swarm literature test the model’s accuracy and decision time only
on setups where all potential site choices are equidistant from the swarm’s
starting location. These models do not account for the geographic challenges
that result from site choices with different geometry. For example, although
actual ant colonies are capable of consistently choosing a higher quality,
further site instead of a lower quality, closer site, existing models are much
less accurate in this scenario. Existing models are also more prone to
committing to a low quality site if it is on the path between the agents’
starting site and a higher quality site. We present a new model for the site
selection problem and verify via simulation that is able to better handle
these geographic challenges. Our results provide insight into the types of
challenges site selection models face when distance is taken into account. Our
work will allow swarms to be robust to more realistic situations where sites
could be distributed in the environment in many different ways.

##  1 Introduction

Swarms of birds, bees, and ants are able to coordinate themselves to make
decisions using only local interactions [3, 14, 18]. Modelling these natural
swarms has inspired many successful swarm algorithms [7]. One such bio-
inspired algorithm comes from the house hunting behavior of ants. Models of
the ants’ behavior when selecting a new nest serve as the basis for swarm
algorithms which seek to select the best site out of a discrete number of
candidate sites in space [17].

Many variations of the best-of-N site selection problem have been studied for
swarms [22]. For example, when sites are of equal quality, choosing one is a
symmetry-breaking problem [9, 23]. Situations with asymmetric site qualities
and costs (where higher quality sites have a higher cost of being chosen) have
also been studied – for example, when one of two candidate sites is
significantly larger than the other (making it harder for agents to detect
other agents favoring the larger site, even when it is of higher quality) [4].

However, most site selection models are mainly tested on small numbers of
candidate nest sites that are equidistant from the agents’ starting location
(also known as the home nest) [6, 15]. In many applications of the site
selection problem such as shelter seeking, sites will not be distributed so
uniformly.

This equidistant setup fails to capture two important geographical details
that existing algorithms struggle with in making accurate decisions. Firstly,
nests that are closer to the home nest are advantaged because they are more
likely to be found. Even so, house hunting ants can still choose higher
quality sites that are much further than lower quality, closer sites. We have
found that existing site selection models often commit to the closer site even
when there is a better, further option. Secondly, using sites equidistant from
the home nest eliminates the possibility of some nests being in the way of
others. Site selection models often trigger consensus on a new site after a
certain quorum population of agents have been detected in it. If a low quality
nest is on the path from the home nest to a high quality nest, agents
travelling between the home nest and the high quality nest could saturate the
path and detect a quorum for the lower quality nest that is in the way instead
of the highest quality nest.

This paper aims to create a new algorithm that can successfully account for a
more varied range of nest distributions, allowing agents to successfully
choose higher quality nests even when they have the disadvantage of being
further from the agents’ starting location or there are other lower quality
nests in the way. The model should also perform with similar accuracy compared
to existing models on the default setup with equidistant candidate nest sites.
We show via simulation that incorporating a quorum threshold that decreases
with site quality allows for increased accuracy compared to previous models.
We also show that setups where candidate sites are in the way of each other or
are of similar quality can make it harder for site selection models to produce
accurate results.

Section 2 describes the house hunting process of ants and overviews existing
swarm models. Section 3 describes our model. We provide details on the
implementation of our model, test accuracy and decision time in different
geographic situations, and report the results in Section 4. We discuss these
results in Section 5. Lastly, we suggest future work in Section 6. The full
simulation code can be found at [1].

##  2 Background

###  2.1 Ant House Hunting

When the T. albipennis ants’ home nest is destroyed, the colony can find and
collectively move to a new, high quality nest. To do so, T. albipennis scouts
first scan the area, searching for candidate nests. When a nest is found, the
scouts wait a period of time inversely proportional to the nest quality before
returning to the home nest. There, they recruit others to examine the new site
in a process known as forward tandem running. Tandem runs allow more ants to
learn the path to a new site in case the ants decide to move there. When an
ant in a candidate nest encounters others in the site at a rate surpassing a
threshold rate (known as the quorum threshold), ants switch their behavior to
carrying other members of the colony to the new nest. Carrying is three times
faster than tandem runs and accelerates the move to the new nest [13, 14].

This decision-making process allows ants to not only agree on a new nest, but
also to choose the highest quality nest out of multiple nests in the
environment. This is true even if the high quality nest is much further from
the home nest than the low quality nest [8, 19]. Franks [8] found that with a
low quality nest 30 cm from home and a high quality nest 255cm from home, 88%
of ant colonies successfully chose the high quality nest even though it was 9
times further.

###  2.2 House Hunting and Site Selection Models

To better study the ants’ behavior, models have been designed to simulate how
ants change behavior throughout the house hunting process [15, 24]. These
models, initiated by Pratt [15], allow simulated ants to probabilistically
transition through four phases – the Exploration, Assessment, Canvassing, and
Transport phases. The Exploration represents when the ants are still exploring
their environment for new sites. When an ant discovers a site, it enters the
Assessment phase, in which it examines the quality of the site and determine
whether to accept or reject it. If the ant accepts the site, it enters the
Canvassing phase, which represents the process of recruiting other ants via
forward tandem runs. Finally, if a quorum is sensed, the ant enters the
Transport phase, which represents the carrying behavior used to move the
colony to the new site.

These models, however, assume that when an ant transitions from the
Exploration phase to the Assessment phase, it is equally likely to choose any
of the candidate nest sites to assess. This assumes that any nest is equally
likely to be found, which is unlikely in the real world because sites closer
to the home nest are more likely to be discovered. To our knowledge, house
hunting models have not tried to model situations where nests have different
likelihoods of being found, as is the case when nests have different distances
from the home nest [24].

The corresponding problem to house hunting in robot swarms is known as the
$N$-site selection problem [22]. Agents, starting at a central home site, must
find and choose among $N$ candidate nest sites in the environment and move to
the site with highest quality. Unlike house hunting models, which do not
physically simulate ants in space, swarm models set up ants in a simulated
arena and let them physically explore sites and travel between them.

Inspired by ant modelling, [6] and [17] have modeled swarm agents using four
main states – Uncommitted Latent, Uncommitted Interactive, Favoring Latent,
and Favoring Interactive (with [6] adding a fifth Committed state to emulate
having detected a quorum). Uncommitted Latent agents remain in the home nest
while Uncommitted Interactive agents explore the arena for candidate sites.
Favoring Interactive agents have discovered and are favoring a certain site
and recruit other agents to the site, while Favoring Latent agents remain in
the new site to try and build up quorum. Agents probabilistically transition
between these states based on environmental events (e.g. the discovery of a
new site) and eventually end up significantly favoring a new candidate nest or
committed to it. Other swarm models for N-site selection typically use a
similar progression through uncommitted, favoring, and committed type phases
[12].

One setup where a high quality site was twice as far as a low quality one was
successfully solved in [17], but for the most part these models and their
variations have mainly been tested in arenas with two candidate sites
equidistant from the home nest [2, 6, 11, 16]. Our model aims to analyze the
behavior of these models in more varied site setups and and improve upon them.

##  3 Model

We first describe our new discrete geographical model for modeling swarms.
Then we discuss the individual restrictions, parameters, and agent algorithms
needed for the house hunting problem specifically.

A more formal description of our general model (described in Section 3.1) can
be found in Appendix A.

###  3.1 General Model

We assume a finite set $R$ of agents, with a state set $SR$ of potential
states. Agents move on a discrete rectangular grid of size $n\times m$,
formally modelled as directed graph $G=(V,E)$ with $|V|=mn$. Edges are
bidirectional, and we also include a self-loop at each vertex. Vertices are
indexed as $(x,y)$, where $0\leq x\leq n-1$, $0\leq y\leq m-1$. Each vertex
also has a state set $SV$ of potential states.

We use a discrete model so the model can be simulated in a distributed fashion
on each vertex to reduce computation time.

####  3.1.1 Local Configurations:

A local configuration $C^{\prime}(v)$ captures the contents vertex $v$. It is
a triple $(sv,myagents,srmap)$, where $sv\in SV$ is the vertex state of $v$,
$myagents\subseteq R$ is the set of agents at $v$, and
$srmap:myagents\rightarrow SR$ assigns an agent state to each agent at $v$.

####  3.1.2 Local Transitions:

The transition of a vertex $v$ may be influenced by the local configurations
of nearby vertices. We define an influence radius $I$, which is the same for
all vertices, to mean that vertex indexed at $(x,y)$ is influenced by all
valid vertices $\\{(a,b)|a\in[x-I,x+I],b\in[y-I,y+I]\\}$, where $a$ and $b$
are integers. We can use this influence radius to create a local mapping
$M_{v}$ from local coordinates to the neighboring local configurations. For a
vertex $v$ at location $(x,y)$, we produce $M_{v}$ such that
$M_{v}(a,b)\rightarrow C^{\prime}(w)$ where $w$ is the vertex located at
$(x+a,y+b)$ and $-I<a,b<I$. This influence radius is representative of a
sensing and communication radius. Agents can use all information from vertices
within the influence radius to make decisions.

We have a local transition function $\delta$, which maps all the information
associated with one vertex and its influence radius at one time to new
information that can be associated with the vertex at the following time. It
also produces directions of motion for all the agents at the vertex.

Formally, for a vertex $v$, $\delta$ probabilistically maps $M_{v}$ to a
quadruple of the form $(sv_{1},myagents,srmap_{1},dirmap_{1})$, where
$sv_{1}\in SV$ is the new state of the vertex, $srmap_{1}:myagents\rightarrow
SR$ is the new agent state mapping for agents at the vertex, and
$dirmap_{1}:myagents\rightarrow\\{R,L,U,D,S\\}$ gives directions of motion for
agents currently at the vertex. Note that $R$, $L$, $U$, and $D$ mean right,
left, up, and down respectively, and $S$ means to stay at the vertex. The
local transition function $\delta$ is further broken down into two phases as
follows.

Phase One: Each agent in vertex $v$ uses the same probabilistic transition
function $\alpha$, which probabilistically maps the agent’s state $sr\in SR$,
location $(x,y)$, and the mapping $M_{v}$ to a new suggested vertex state
$sv^{\prime}$, agent state $sr^{\prime}$, and direction of motion
$d\in\\{R,L,U,D,S\\}$. We can think of $\alpha$ as an agent state machine
model.

Phase Two: Since agents may suggest conflicting new vertex states, a rule $L$
is used to select one final vertex state. The rule also determines for each
agent whether they may transition to state $sr^{\prime}$ and direction of
motion $d$ or whether they must stay at the same location with original state
$sr$.

####  3.1.3 Probabilistic Execution:

The system operates by probablistically transitioning all vertices $v$ for an
infinite number of rounds. During each round, for each vertex $v$, we obtain
the mapping $M_{v}$ which contains the local configurations of all vertices in
its influence radius. We then apply $\delta$ to $M_{v}$ to transition vertex
$v$ and all agents at vertex $v$. For each vertex $v$ we now have
$(sv_{v},myagents_{v},srmap_{v},dirmap_{v})$ returned from $\delta$.

For each $v$, we take $dirmap_{v}$, which specifies the direction of motion
for each agent and use it to map all agents to their new vertices. For each
vertex $v$, it’s new local configuration is just the new vertex state
$sv_{v}$, the new set of agents at the vertex, and the $srmap$ mapping from
agents to their new agent states.

###  3.2 House Hunting Model

The goal of the house hunting problem is for agents to explore the grid and
select the best site out of $N$ sites to migrate to collectively. We model
sites as follows.

A set $S$, $|S|=N$ of rectangular sites are located within this grid, where
site $s_{i}$ has lower left vertex $(x_{i}^{1},y_{i}^{1})$ and upper right
vertex $(x_{i}^{2},y_{i}^{2})$. Each site $s_{i}$ also has a quality
$s_{i}.q\in[0,1]$. To represent these sites, we let the vertex state set be
$SV=S\cup\\{\emptyset\\}$ for each vertex, indicating which site, if any, the
vertex belongs to. Furthermore, we denote the site $s_{0}$ to be the home
nest. In the initial configuration, all agents start out at a random vertex in
the home nest, chosen uniformly from among the vertices in that nest.

###  3.3 Agent States and Transition Function

The agent state set $SR$ is best described in conjunction with the agent
transition function $\alpha$. Agents can take on one of $6$ core states, each
a combination of one of three preference states (Uncommitted, Favoring,
Committed), and two activity states (Nest, Active). The state model can be
seen in Figure 1.

$U^{A}$ $U^{N}$ $F^{A}_{i}$ $F^{N}_{i}$ $F^{N}_{j\neq i}$
$P_{N}$$P_{A}$$P_{S_{i}}v_{i}(1-x)$$\beta$$P_{N}$$P_{A}$$Ds_{j}P_{ij}$$Ds_{i}P_{ji}$$P_{S_{i}}v_{i}x$

$U^{A}$ $U^{N}$ $F^{N}_{i}$ $F^{A}_{i}$ $C^{A}_{i}$ $C^{N}_{i}$
$0.5(P_{Q_{i}}+P_{Q})$$0.5(P_{Q_{i}}+P_{Q})$$0.5(P_{Q_{i}}+P_{Q})$$0.5(P_{Q_{i}}+P_{Q})$$0.5P_{Q_{i}}$$0.5P_{Q_{i}}$$0.5P_{Q_{i}}$$0.5P_{Q_{i}}$

Figure 1: State model. $\\{U,F,C\\}$ denote preference states. The superscript
$\\{N,A\\}$ denotes the activity state, and a subscript $i$ denotes that an
agent is favoring or committed to site $i$. The transitions for Uncommitted
and Favoring states are shown on the left, and transitions from Uncommitted
and Favoring to Committed states are on the right.

Uncommitted Nest ($U^{N}$) agents stay in the home nest to prevent too many
agents from flooding the environment. They have a chance of transitioning to
Uncommitted Active ($U^{A}$) agents, which try to explore the arena and
discover new sites. $U^{A}$ agents move according to the Levy flight random
walk111The Levy flight is expressed within our general model by breaking it
down into individual steps, and progress along the flight is maintained in the
agent’s state. See [1] for implementation details., which has been shown to be
used by foraging ants [21]. $U^{N}$ agents transition to $U^{A}$ with
probability $P_{A}$, and $U^{A}$ agents transition to $U^{N}$ agents with
probability $P_{N}$. This results in an expected $x=\frac{P_{A}}{P_{A}+P_{N}}$
percent of uncommitted agents are active, whereas $1-x$ agents remain in the
nest. Prior work [17] lets $P_{N}=9P_{A}=L$, where $L$ is the inverse of the
average site round trip time, chosen to promote sufficient mixing. This leads
to $10\%$ of the agent population being active.

Uncommitted Active agents have a chance $P_{S_{i}}$ of discovering a new nest,
which is 1 if a new nest is within influence radius and $0$ otherwise. If they
discover a nest $s_{i}$, they explore and accept it with probability $s_{i}.q$
(the quality of $s_{i}$). They then have an $x\%$ chance of transitioning to
Favoring Active, and a $(1-x)\%$ chance of transitioning to Favoring Nest.

Favoring agents $(F^{A}_{i},F^{N}_{i})$ prefer the site $s_{i}$ that they
discovered. Favoring Active $(F^{A}_{i})$ agents remain in site $s_{i}$ to
build quorum. Favoring Nest $(F^{N}_{i})$ agents return to the home nest to
recruit others to site $s_{i}$. Favoring Nest agents transition to Active with
the same probability $P_{A}$ and Favoring Active agents transition back to
Nest agents with probability $P_{N}$, creating the same effect where an
expected $90\%$ of the favoring agent population is $F^{A}_{i}$ while the rest
are $F^{N}_{i}$.

$F^{N}_{i}$ agents have a probability $\beta$ of abandoning their nest, which
is 1 if the time spent without seeing other agents surpasses $t_{\beta}$.
$F^{A}_{i}$ agents can be inhibited by other $F^{A}_{i}$ agents as follows.
The chance an agent favoring nest $i$ is converted to favoring nest $j$ is
$Dr_{j}P_{ij}$, where the factor of $D$ is the probability of agents messaging
each other (to prevent excessive messaging). $r_{j}$ is the number of agents
favoring $s_{j}$ that have the agent within their influence radius. After an
agent hears of the new site $s_{j}$, it visits the site to evaluate $s_{j}.q$
and changes its preference to $s_{j}$ if $s_{j}.q>s_{i}.q$. Thus, the
condition $P_{ij}$ is $1$ when $s_{j}.q>s_{i}.q$ and $0$ otherwise.

$U^{A}$ agents and $F^{N}_{i}$ agents can detect a quorum and commit to a site
when $q$ agents in the site are within their influence radius. The quorum size
scales with site value as $q=\lfloor(q_{MIN}-q_{MAX})*s_{i}.q+q_{MAX}\rfloor$,
where $q_{MAX}$ and $q_{MIN}$ are the maximum and minimum possible quorum
threshold respectively. The condition $P_{Q}$ is $1$ when quorum is satisfied
and $0$ otherwise. Agents in any Favoring or Uncommitted state will transition
to the committed state, if they encounter an agent already in quorum. The
condition $P_{Q_{i}}$ is $1$ when another quorum agent for $s_{i}$ is
encountered and $0$ otherwise. Furthermore, agents have an $\frac{1}{2}$
chance of transitioning to Committed Active $(C^{A}_{i})$ and a
$\frac{1}{2}\%$ chance of Committed Nest $C^{N}_{i}$ after having detected or
been notified of a quorum.

$C^{N}_{i}$ agents head to the home nest to inform others of the move, while
$C^{A}_{i}$ agents randomly wander the grid to find stragglers. Agents in
quorum states continue to wander until they have sensed quorum for $t_{Q}$
time steps, whereupon they return to the new selected site $s_{i}$.

The resulting agent state set $SR$ is a product of the $6$ core states needed
in the state model as well as a number of auxiliary variables such as an
agent’s destination, the names of the sites it favors or has sensed quorum
for, and parameters for an agent’s random walk when exploring the grid.

Since in the house hunting problem (unlike other problems like task
allocation), an agent never modifies the environment, an agent’s proposed new
vertex state is always the same as the old vertex state. Therefore, phase two
of $\delta$ is not needed to reconcile conflicting vertex state suggestions
from agents.

The transition function $\alpha$, which for each agent returns a proposed new
vertex state $sv^{\prime}$, agent state $sr^{\prime}$ and direction of motion
works as follows. The agent never modifies the grid, so $sv^{\prime}=sv$. The
agent state $sr^{\prime}$ and direction $d$ are calculated according to the
core transitions and the auxiliary variables needed to keep track of those
transitions. For example, when an agent is headed towards a site, the
direction $d$ is calculated to try to match a straight-line traversal to the
site from the agent’s starting location as close as possible. When an agent is
staying within a site, the direction $d$ is calculated to be a random walk
within the site boundaries.

The total set of variables parameters is
$\\{P_{A},P_{N},D,t_{Q},t_{\beta},q_{MIN},q_{MAX}\\}$, as well as the site
locations $(x_{i}^{1},y_{i}^{1}),(x_{i}^{2},y_{i}^{2})$ and quality $s_{i}.q$.
In Section 4, we explore how changes in $q_{MIN}$, $q_{MAX}$, and the site
locations and quality impact the accuracy, decision time, and split decisions
made by the model.

##  4 Results

The model was tested in simulation using Pygame, with each grid square
representing 1cm2. Agents moved at 1 cm/s, with one round representing one
second. We chose this speed because even the lowest cost robots are still able
to move at 1cm/s [20]. Agents had an influence radius of $2$. All simulations
were run using 100 agents, and a messaging rate of $1/15$. We let the
abandonment timeout $t_{\beta}=\frac{5}{L}$ and the quorum timeout
$t_{Q}=\frac{1}{L}$.

For each set of trials, we evaluated accuracy (measured as the fraction of
agents who chose the highest quality nest), decision time (measured as the
time it took for all agents to arrive at the nest they committed to), and
split decisions (measured as the number of trials where not all agents
committed to the same nest).

###  4.1 Further Nest of Higher Quality

House hunting ants are capable of choosing further, higher quality sites over
closer, lower quality ones [8]. When the far site and the near site are of
equal value, ants consistently choose the closer one. To test our model’s
ability to produce the same behavior, we replicated the experimental setups in
[8].

Three different distance comparisons were tested, with a further, higher
quality nest of quality $0.9$ being 2x, 3x, and 9x as far as a lower quality
nest of quality $0.3$ on the path from the high quality nest to the home nest.
We included a control setup for each of these distance comparisons where both
the far and close nest were quality $0.3$. The arena size was $N=16,M=80$ for
the 2x case, $N=18,M=180$ for the 3x case, and $N=18$, $M=300$ for the 9x
case.

We tested our model using two different quorum parameters. In one test, we had
$q_{MIN}=q_{MAX}=4$, intended to represent the behavior of previous models
with a fixed quorum threshold. In the other setup, $q_{MIN}=4$ and
$q_{MAX}=7$, allowing our model to use the new feature of scaling the quorum
threshold with site quality. We ran 100 trials for each set of parameters.

Figure 2: Decision Time and Accuracy for far nests 2, 3, and 9 times as far
from the home nest. Fixed quorum indicates the fixed threshold value of $4$,
and scaled quorum indicates $q_{MIN}=4,q_{MAX}=7$. The accuracy for the actual
ants is taken from [8].

As seen in Figure 2, using a scaled threshold significantly improved accuracy
from using a fixed one. In the control case, both the fixed and scaled quorum
threshold achieved high accuracy, with all accuracies being greater than
$99\%$. In cases where the far site was of higher quality, the decision time
for fixed and scaled quorum was comparable. However, the scaled quorum
threshold took significantly (Welch’s T-test, p=0.05) more time in the control
case to decide.

Furthermore, as seen in Figure 2, our model successfully chose the further
site with comparable (or significantly higher in the 9x case) accuracy than
ants themselves, indicating that our model is on par with the ants.

###  4.2 Effects of Lower Quality Nest Being In the Way

To isolate the effects of the low quality nest being in the way of the high
quality nest, we tested our model where the high quality nest (quality $0.9$)
was one of $\\{2,3,4,5,6,7,8,9\\}$ times further than the low quality nest
(quality $0.3$), but in opposite directions of the home nest. We compared
model performance when the low quality nest was in the way of the home nest.
We ran tests with $N=18$, $M=300$, with the low quality nest always $30$cm
from home. We again tested a fixed ( $q_{MIN}=q_{MAX}=4$) and
scaled($q_{MIN}=4,q_{MAX}=7$) quorum threshold on these setups. $100$ trials
were conducted for each set of parameters.

Figure 3: Decision Time and Accuracy for far nests $2-9$ times further than
the close nest for both fixed and scaled quorums. In the in-the-way setup, the
home nest, low quality nest, and high quality nest were lined up in that
order. In the out of way setup, the low quality nest, home nest, and high
quality nest were lined up in that order.

Figure 3 shows that for the out-of-way setup, the scaled quorum performs
significantly (Welch’s T-test, p=0.05) more accurately than the fixed quorum
on all far nest distances. For the in-the-way setup, the scaled quorum
performs significantly better (Welch’s T-test, p=0.05) when the far nest is
$3$x further or more. Note it is harder for the fixed quorum to solve the in-
the-way problem accurately compared to the out-of-way problem (Welch’s T-test,
p=0.05). It is likewise harder for the scaled quorum to solve the in-the-way
problem when the far nest is $\\{3,4,6,7,8,9\\}$ times further (Welch’s
T-test, p=0.05), showing that the in-the-way problem is harder to solve for
site selection algorithms.

For distances $3$x or further, there is no significant difference between the
decision times for the fixed out-of-way, scaled out-of-way, and scaled in-the-
way setups. For distances $5$x and further, the fixed quorum takes
significantly less time than the other setups but suffers in decision accuracy
(Welch’s T-test, p=0.05) compared to the other three setups.

###  4.3 Effects of Magnitude of Difference in Site Quality

Because site quality affects the quorum threshold, we expect it to be harder
for agents to correctly choose a high quality far site when it is only
slightly better than than nearby lower quality sites. This is because the
difference in quorum threshold is less pronounced for sites of similar
quality. For two equidistant nests, the algorithm should consistently choose
the best site as it has in past work, so the absolute difference in site
quality should not matter.

To test these effects, we used the setup in Section 4.1 where the further nest
was $2$x (60 cm) as far as the in-the-way close nest (30 cm), and compared it
to an equidistant setup where both candidate nests were $30$ cm away from the
home nest in opposite directions. We tested both a fixed quorum
$q_{MAX}=q_{MIN}=4$ and a scaled quorum on these setups $q_{MAX}=7,q_{MIN}=4$.
We varied the quality of the near nest in the set of potential values
$\\{0.3,0.6,0.9\\}$, corresponding to quorum thresholds of $\\{6,5,4\\}$
respectively, with the far nest having quality $1.0$. (In the equidistant
case, we varied the quality of one nest while the other had quality $1.0$.)
Figure 4 shows the resulting accuracy and decision time.

Figure 4: Decision Accuracy and Time given varying differences in site quality
between the near and the far nest.

As predicted, a smaller difference in site quality / quorum threshold led to
significantly (Welch’s T-test, p=0.05) lower decision accuracy for the non-
equidistant setup. In the equidistant setup, agents were able to achieve a
near-100% outcome regardless of magnitude of differences in site quality.
However, in the unbalanced setup, we confirmed that for larger differences in
site quality, the algorithm comes to a more accurate decision, showing that
non-equidistant candidate nest setups cause sensitivity to absolute site value
differences that can’t be seen in the equidistant setup.

##  5 Discussion

The results demonstrate our model’s ability to improve accuracy when choosing
from a higher quality, further site and a lower quality, closer site. This
improvement comes at the cost of a higher decision time when converging on a
lower quality site, because the quorum threshold for low quality sites is
higher in our model. This higher decision time is reasonable and represents
hesitance when committing to a poor quality option in the hopes of finding a
better one.

Our model also demonstrated the extra difficulty that comes with a lower
quality site being in the path from the home nest to a high quality site.
Qualitative observation showed that agents travelling back and forth between
the far site and the home nest often unintentionally contributed to a quorum
in the poor quality, in-the-way site as they travelled through it. We showed
that using a scaled quorum threshold as opposed to a fixed one is an effective
way of significantly increasing decision accuracy. However, even if the
closer, poor quality site is completely out of the way of the far, high
quality site, Figure 3 shows that using a scaled quorum can still help to
improve accuracy.

Figure 4 shows that our model is still successful when candidate sites are
equidistant from home, as is most commonly tested. We also show that an
equidistant setup is not influenced by the absolute difference between
candidate site qualities. Contrarily, in the setup with a further, high
quality nest, it is harder to make an accurate decision the smaller the
quality difference between the high and low quality nests. Note that it is
also less grievous of an error to choose the low quality nest when the quality
difference is small.

We observed a shorter decision time in conjunction with lower accuracy,
similar to the time-accuracy trade off in natural swarms [5, 10]. In each set
of 100 trials run on our setups had at most $2$ split decisions, indicating
our model is successful in keeping the swarm together even when migrating to
the further nest.

##  6 Future Work

While our model has made strides in being more accurate when choosing between
sites with different geographical distributions, many site setups have yet to
be tested. Future work could introduce obstacles to the environment, or try to
adapt the house hunting model to an arena with continuous site values.

Our model suggests that a quorum threshold that scales with site quality leads
to more accurate site selection. Future work could explore if actual ants do
the same and use this information to create more accurate models.

Lastly, while our model is hard to analyze without making simplifications
(because it involves agents physically moving in space), future work could try
to develop analytical bounds. One method we envision is simplifying the chance
of each site being discovered to a fixed probability and trying to model agent
population flow between the different model states, similar to [16], which
does this for candidate sites all with an equal chance of discovery.

## References

  * [1] Cai, G.: Geometric Swarm Modelling (4 2022). https://doi.org/10.5281/zenodo.6508646
  * [2] Cai, G., Sofge, D.: An urgency-dependent quorum sensing algorithm for n-site selection in autonomous swarms. In: Proceedings of the 18th International Conference on Autonomous Agents and MultiAgent Systems. pp. 1853–1855 (2019)
  * [3] Camazine, S., Visscher, P.K., Finley, J., Vetter, R.S.: House-hunting by honey bee swarms: collective decisions and individual behaviors. Insectes Sociaux 46(4), 348–360 (1999)
  * [4] Campo, A., Garnier, S., Dédriche, O., Zekkri, M., Dorigo, M.: Self-organized discrimination of resources. PLoS One 6(5), e19888 (2011)
  * [5] Chittka, L., Dyer, A.G., Bock, F., Dornhaus, A.: Bees trade off foraging speed for accuracy. Nature 424(6947), 388–388 (2003)
  * [6] Cody, J.R., Adams, J.A.: An evaluation of quorum sensing mechanisms in collective value-sensitive site selection. In: 2017 International Symposium on Multi-Robot and Multi-Agent Systems (MRS). pp. 40–47 (2017). https://doi.org/10.1109/MRS.2017.8250929
  * [7] Fan, X., Sayers, W., Zhang, S., Han, Z., Ren, L., Chizari, H.: Review and classification of bio-inspired algorithms and their applications. Journal of Bionic Engineering 17(3), 611–631 (2020)
  * [8] Franks, N.R., Hardcastle, K.A., Collins, S., Smith, F.D., Sullivan, K.M., Robinson, E.J., Sendova-Franks, A.B.: Can ant colonies choose a far-and-away better nest over an in-the-way poor one? Animal Behaviour 76(2), 323–334 (2008)
  * [9] Hamann, H., Schmickl, T., Wörn, H., Crailsheim, K.: Analysis of emergent symmetry breaking in collective decision making. Neural Computing and Applications 21(2), 207–218 (2012)
  * [10] Heitz, R.P.: The speed-accuracy tradeoff: history, physiology, methodology, and behavior. Frontiers in neuroscience 8, 150 (2014)
  * [11] Khurana, S., Sofge, D.: Quorum sensing re-evaluation algorithm for n-site selection in autonomous swarms. In: ICAART (1). pp. 193–198 (2020)
  * [12] Parker, C.A.C., Zhang, H.: Cooperative decision-making in decentralized multiple-robot systems: The best-of-n problem. IEEE/ASME Transactions on Mechatronics 14, 240–251 (2009)
  * [13] Pratt, S.C.: Behavioral mechanisms of collective nest-site choice by the ant temnothorax curvispinosus. Insectes Sociaux 52(4), 383–392 (2005)
  * [14] Pratt, S.C.: Quorum sensing by encounter rates in the ant temnothorax albipennis. Behavioral Ecology 16(2), 488–496 (2005)
  * [15] Pratt, S.C., Sumpter, D.J., Mallon, E.B., Franks, N.R.: An agent-based model of collective nest choice by the ant temnothorax albipennis. Animal Behaviour 70(5), 1023–1036 (2005)
  * [16] Reina, A., Marshall, J.A., Trianni, V., Bose, T.: Model of the best-of-n nest-site selection process in honeybees. Physical Review E 95(5), 052411 (2017)
  * [17] Reina, A., Valentini, G., Fernández-Oto, C., Dorigo, M., Trianni, V.: A design pattern for decentralised decision making. PloS one 10(10), e0140950 (2015)
  * [18] Reynolds, C.W.: Flocks, herds and schools: A distributed behavioral model. In: Proceedings of the 14th annual conference on Computer graphics and interactive techniques. pp. 25–34 (1987)
  * [19] Robinson, E.J., Smith, F.D., Sullivan, K.M., Franks, N.R.: Do ants make direct comparisons? Proceedings of the Royal Society B: Biological Sciences 276(1667), 2635–2641 (2009)
  * [20] Rubenstein, M., Ahler, C., Nagpal, R.: Kilobot: A low cost scalable robot system for collective behaviors. In: 2012 IEEE international conference on robotics and automation. pp. 3293–3298. IEEE (2012)
  * [21] Sims, D.W., Humphries, N.E., Bradford, R.W., Bruce, B.D.: Lévy flight and brownian search patterns of a free-ranging predator reflect different prey field characteristics. Journal of Animal Ecology 81(2), 432–442 (2012)
  * [22] Valentini, G., Ferrante, E., Dorigo, M.: The best-of-n problem in robot swarms: Formalization, state of the art, and novel perspectives. Frontiers in Robotics and AI 4, 9 (2017)
  * [23] Wessnitzer, J., Melhuish, C.: Collective decision-making and behaviour transitions in distributed ad hoc wireless networks of mobile robots: Target-hunting. In: European Conference on Artificial Life. pp. 893–902. Springer (2003)
  * [24] Zhao, J., Lynch, N., Pratt, S.C.: The power of social information in ant-colony house-hunting: A computational modeling approach. bioRxiv pp. 2020–10 (2021)

## Appendix A: Formal General Model

Below we present an alternative, more formal description of our general model
from section 3.1.

###  6.1 General Model

We assume a finite set $R$ of agents, all of which have a state set $SR$ of
potential states. Agents move on a discrete rectangular grid of size $n\times
m$, formally modelled as directed graph $G=(V,E)$ with $|V|=mn$. Edges are
bidirectional, and we also include a self-loop at each vertex. Vertices are
indexed as $(x,y)$, where $0\leq x,y\leq n-1$. Each vertex also has a state
set $SV$ of potential states.

###  6.2 Grid and Agent Setup

We define four kinds of configurations, global vs. local, and ordinary vs.
transitory. The transitory configurations are used as intermediate steps in
defining system executions.

##### (Global) configurations:

A _(global) configuration_ $C$ is a triple of mappings,
$(svmap,srmap,locmap)$, where:

  * •

$svmap:V\rightarrow SV$ is the vertex state mapping, which assigns a vertex-
state to each vertex,

  * •

$srmap:R\rightarrow SR$ is the agent state mapping, which assigns an agent-
state to each agent, and

  * •

$locmap:R\rightarrow V$ is the location mapping, which assigns a location to
each agent.

##### Local configurations:

A _local configuration_ $C^{\prime}$ is intended to capture the contents of
one vertex/square). It is a triple $(sv,myagents,srmap)$, where:

  * •

$sv\in SQ$ is the vertex-state of the given vertex,

  * •

$myagents\subseteq R$ is the set of agents at the vertex, and

  * •

$srmap:myagents\rightarrow SR$ assigns an agent-state to each agent at the
vertex.

Define an operator $project(C,v)$, which takes a global configuration $C$ and
projects it to give a local configuration of vertex $v$. Namely, if
$C=(svmap,srmap,locmap)$, and $v$ is any vertex, then $project(C,v)$ is the
local configuration $(svmap(v),R^{\prime},srmap\lceil R^{\prime})$, where
$R^{\prime}=\\{r\in R:locmap(r)=v\\}$. That is, we pick out the state of the
given vertex $v$, and the states for all the agents located at $v$.

We also define an inverse operation called $merge$, which combines local
configurations to yield a global configuration. We assume a collection of
"compatible" local configurations $\\{C^{\prime}(v):v\in V\\}$, where
"compatible" here means that the unique vertices below exist. The resulting
global configuration $C$ is defined by:

  * •

For every $v\in V$, $svmap(v)=C^{\prime}(v).sv$, that is, take the state of
$v$ from the local configuration belonging to vertex $v$.

  * •

For every $r\in R$, $locmap(r)$ is the unique vertex $v$ such that $r\in
C^{\prime}(v).myagents$, and $srmap(r)=C^{\prime}(locmap(r)).srmap(r)$.

###  6.3 Transitory Configurations

We also have a notion of transitory configuration, which is used as an
intermediate stage between two ordinary configurations, in constructing
executions. It represents agents in motion from one vertex to another.

A _transitory configuration_ $T$ is a triple of mappings
$(svmap,srmap,edgemap)$, where

  * •

$svmap$ and $srmap$ have the same types as for ordinary configurations, and

  * •

$edgemap:R\rightarrow E$ assigns a directed edge to each agent.

This represents newly-computed states for all vertex and agents, plus the
edges on which the agents are traveling to get to their next locations.

##### Local transitory configurations:

A _local transitory configuration_ $T^{\prime}$ is a quadruple
$(sv,myagents,srmap,dirmap)$, where

  * •

$sv$, $myagents$, and $srmap$ are as in the definition of a local
configuration, and

  * •

$dirmap:myagents\rightarrow\\{R,L,U,D,S\\}$ is the direction mapping, which
assigns a travel direction to each agent. Direction $S$ means to stay at the
vertex.

This represents newly-computed states for a single vertex and its agents, plus
directions of travel for the local agents.

As before, define an operator $project(T,v)$, which takes a global transitory
configuration $T$ and projects it to give a local transitory configuration
$T(v)$ for vertex $v$. Namely, if $T=(svmap,srmap,edgemap)$, and $v$ is any
vertex, then $project(T,v)$ is the local transitory configuration
$(svmap(v),myagents,srmap\lceil myagents,dirmap)$, where:

  * •

$myagents=\\{r\in R:edgemap(r)\mbox{ is an edge with }source=v\\}$, and

  * •

For each $r\in myagents,dirmap(r)$ is the direction of the edge $edgemap(r)$,
that is, $R$ if the edge goes from $v$ to its right neighbor, etc.

We also define an inverse, $merge$ operation, which combines local transitory
configurations to give a global transitory configuration. We assume a
collection of "compatible" local configurations $\\{T^{\prime}(v):v\in V\\}$,
where "compatible" means that the unique vertices below exist. The resulting
global transitory configuration $T$ is defined by:

  * •

For every $v\in V$, $svmap(v)=T^{\prime}(v).sv$,

  * •

For every $r\in R$, let $v$ be the unique vertex such that $r\in
T^{\prime}(v).myagents$; then define:

    * –

$srmap(r)=T^{\prime}(v).srmap(r)$.

    * –

$edgemap(r)$ is the directed edge from $v$ to $w$, where $w$ is the vertex
that is the target of the edge from $v$ in the direction given by
$T^{\prime}(v).dirmap(r)$.

Given a transitory configuration $T=(svmap,srmap,edgemap)$, we can map it to
an ordinary configuration using $C=move(T)$ , which is the same as $T$ except
that, instead of edgemap, we have locmap, where for each $r\in R$,
$locmap(r)=target(edgemap(r))$.

###  6.4 Local transitions

The transition of a vertex $v$ may be influenced by the local configurations
of nearby vertices in addition to itself. We define an influence radius $I$,
which is the same for all vertices, to mean that vertex indexed at $(x,y)$ is
influenced by all vertices $\\{(a,b)|a\in[x-I,x+I],b\in[y-I,y+I]\\}$, where
$a$ and $b$ are integers mod $n$. We can use this influence radius to create a
local mapping $M_{v}$ from local coordinates to the neighboring local
configurations. Thus, for a vertex $v$ at location $(x,y)$, we produce $M_{v}$
such that $M_{v}(a,b)\rightarrow C^{\prime}(w)$ where $w$ is the vertex
located at $(x+a,y+b)$ and $-I<a,b<I$. This influence radius is representative
of a sensing and communication radius for agents.

We have a local transition function $\delta$, which maps all the information
associated with one vertex and its influence radius at one time to new
information that can be associated with the vertex at the following time. It
also produces directions of motion for all the agents at the vertex.

Formally, for a vertex $v$, $\delta$ maps the mapping $M_{v}$ to a probability
distribution on local transitory configurations of the form
$(sv_{1},myagents,srmap_{1},dirmap_{1})$, where:

  * •

$sv_{1}\in SV$ is the new state of the vertex,

  * •

$srmap_{1}:myagents\rightarrow SR$ is the new agent state mapping, for agents
currently at the vertex, and

  * •

$dirmap_{1}:myagents\rightarrow\\{R,L,U,D,S\\}$ gives directions of motion for
all the agents currently at the vertex.

###  6.5 Local transition function $\delta$

The local transition function $\delta$ is further broken down into two phases
as follows.

##### Phase One:

Each agent in vertex $v$ uses the same probabilistic transition function
$\alpha$, which maps the agent’s state $sr\in SR$, location $(x,y)$, the
vertex state of the location $sv\in SV$, and the mapping $M_{v}$ to a
distribution over new suggested vertex state $sv^{\prime}$, agent state
$sr^{\prime}$, and direction of motion $d\in\\{R,L,U,D,S\\}$.

##### Phase Two:

A rule $L$ is used to reconcile the different vertex states suggested by each
agent at the vertex and select one final vertex state. The rule also
determines for each agent whether they may transition to state $sr^{\prime}$
and direction of motion $d$ or stay at the same location with original state
$sr$.

###  6.6 Probabilistic execution

The system executes probabilistically, moving through an infinite sequence of
configurations $C_{0},C_{1},C_{2},\ldots,$ each derived probabilistically from
the previous one.

We first describe how the system moves from any configuration $C_{t}$ to the
next configuration $C_{t+1}$. This is a two-phase process where the first
phase is probabilistic and the second is deterministic. The first phase yields
a transitory configuration, and the second phase converts that to an ordinary
configuration.

##### Phase 1:

For each vertex $v$ independently, use $project(C_{t},w)\ \forall w\in IR(v)$
to obtain the local mapping $M_{v}$. Then apply $\delta$ to $M_{v}$, and
sample the resulting distribution to select a local transitory configuration
for $v$. Apply $merge$ to combine these independent results for all $v$ as
described above, to obtain a global transitory configuration.

##### Phase 2:

For the transitory configuration $T$ resulting from Phase 1, compute $move(T)$
to get the new configuration.


††thanks<EMAIL_ADDRESS>0000-0002-7880-9454 ††thanks:
<EMAIL_ADDRESS>0000-0002-5872-519X ††thanks:
<EMAIL_ADDRESS>0000-0002-1287-8780

# Transmon Qubit Constraints on Dark Matter-Nucleon Scattering

Anirban Das  Center for Theoretical Physics, Department of Physics &
Astronomy, Seoul National University, Seoul 08826, South Korea  Noah Kurinsky
SLAC National Accelerator Laboratory, 2575 Sand Hill Rd, Menlo Park, CA 94025,
USA  Kavli Institute for Particle Astrophysics and Cosmology, Stanford
University, Stanford, CA 94035, USA  Rebecca K. Leane  SLAC National
Accelerator Laboratory, 2575 Sand Hill Rd, Menlo Park, CA 94025, USA  Kavli
Institute for Particle Astrophysics and Cosmology, Stanford University,
Stanford, CA 94035, USA

(April 30, 2024)

###### Abstract

We recently pointed out that power measurements of single quasiparticle
devices can be used to detect dark matter. These devices have the lowest known
energy thresholds, far surpassing standard direct detection experiments,
requiring energy deposition above only about an meV. We calculate dark matter
induced quasiparticle densities in transmon qubits, and use the latest
transmon qubit measurements that provide one of the strongest existing lab-
based bounds on dark matter-nucleon scattering below about 100 MeV. We
strongly constrain sub-component dark matter, using both a dark matter
population thermalized in the Earth as well as the dark matter wind from the
Galactic halo. We demonstrate future potential sensitivities using devices
with low quasiparticle densities.

††preprint: SLAC-PUB-17769

Direct detection experiments are currently the main detection effort for dark
matter scattering on Earth. As the dark matter wind impinges on Earth’s
surface, dark matter particles can deposit energy inside direct detection
experiments, which search for anomalous recoils of their Standard Model (SM)
targets. This provides a robust lab-based test of the dark matter-SM
scattering cross section. However, direct detection experiments are limited to
dark matter masses above about a GeV for standard analyses assuming nuclear
recoils [1], or MeV-scale masses when exploiting the Migdal effect [2, 3, 4,
5, 6] or electron recoils [7, 8, 9].

To date, no definitive dark matter induced recoil signal has been detected,
motivating a push for new approaches to search for dark matter scattering
outside the known sensitivity of direct detection experiments [10]. To probe
lighter dark matter, or dark matter with less energy that that coming from the
dark matter wind, lower threshold devices are required. New proposals have
been made in recent years on multiple fronts, with new detection possibilities
with superconductors [11, 12, 13, 14, 15, 16], superfluids [17, 18, 19], polar
crystals [20, 21, 22], topological materials [23], and Dirac materials [24,
25, 26, 27, 28].

In recent work [29], we showed that existing low-noise quantum devices provide
strong sensitivity to dark matter with low energy deposition. Construction of
such quantum devices with as low noise backgrounds as possible is currently
under rapid development in fields from quantum computing to astrophysics.
Quantum information is highly susceptible to decoherence, with the presence of
background noise leading to faults in quantum algorithms. While there is an
active field of research dedicated to solving this problem, called quantum
error correction, noise is difficult to remove and currently limits the
performance of quantum computers. Within astrophysics, the desire to image
increasingly cold objects in the far infrared requires new cryogenic space
telescopes that demand low noise equivalent power in their detectors. Broadly
speaking, the most powerful new opportunities for low-threshold dark matter
detection therefore are driven by advances in technology targeted at
applications such as these.

Previously, we specifically pointed out that power measurements of single
quasiparticle devices provide a new detection mechanism for dark matter with
low energy deposition, including both light dark matter, as well as heavier
dark matter that has become bound to the Earth and is therefore not as
energetic as the dark matter wind [29]. Such devices detect quasiparticle
charges as they tunnel out of a superconducting island. Quasiparticle
excitations are produced by broken Cooper pairs, which can be broken by SM
sources, but also in large rates by dark matter energy deposition. As the
Cooper-pair binding energy sets the energy deposition threshold, these devices
lead to sensitivities to dark matter with energy depositions as low as about
an meV, six orders of magnitude below the usual $\sim$keV threshold required
for detection with direct detection experiments.

In this work, we consider the latest quasiparticle measurements from transmon
qubits [30], and apply them to dark matter scattering predictions. A transmon
qubit is a type of superconducting qubit, which is a fundamental unit of
quantum information processing. It is based on a superconducting circuit
consisting of a Josephson junction shunted by a capacitor. Transmon qubits
have exceptionally reduced sensitivity to background noise, and therefore are
ideal low-background devices for low-threshold dark matter detection.

We will show that these latest quantum measurements provide the best lab-based
sensitivity to date to dark matter-nucleon scattering, for dark matter with
mass lighter than about 100 MeV, and offer excellent prospects for future
sensitivities. We will consider the implications for sub-component dark
matter, which applies to dark matter candidates that do not comprise the
entirety of the dark matter energy density observed in the Universe. Proposals
to detect sub-component dark matter have become of increasing interest
recently [31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43]. We will
calculate constraints on such a candidate stronger than any of the previous
proposals, and in some cases clearly surpass spin-independent direct detection
constraints on sub-component dark matter.

Our paper is organized as follows. In Sec. I, we first detail the dark matter
inputs relevant for our calculations. In Sec. II we discuss our transmon qubit
device and dark matter induced quasiparticle production and tunneling rates.
In Sec. III we discuss our new constraints on dark matter-nucleon scattering,
on both the usual dark matter which explains the full energy content of the
Universe, as well as sub-component dark matter candidates. We also calculate
projections for future sensitivities with transmon qubits. We conclude with a
summary of our findings in Sec. IV.

##  I Dark Matter Density and Velocity Inputs

There are two potential dark matter components that can be tested by direct
detection experiments or our transmon qubits. As these two components have
different dark matter density and velocity values, we detail them separately
below.

###  I.1 Dark Matter Wind from the Galactic Halo

Dark matter which scatters with the Earth as it enters from the Galactic halo
is the standard dark matter distribution tested in direct detection
experiments. As such, our calculations using this dark matter distribution
will be directly comparable to existing direct detection experiments, with no
additional assumptions required. We will call this component “halo dark
matter”.

For halo dark matter, we use a velocity distribution given by a truncated
Maxwell-Boltzmann distribution,

$\displaystyle f_{\chi}^{\rm
halo}(\bm{v})=\frac{1}{N_{0}}e^{-((\bm{v}+\bm{v}_{\oplus})/\bm{v}_{\rm
halo})^{2}}\Theta(v_{\mathrm{esc}}^{\rm halo}-(\bm{v}+\bm{v}_{\oplus}))\,,$
(1)

where $N_{0}$ normalizes the distribution, and $\bm{v}_{\rm halo}$ is the
average dark matter velocity in the halo $\bm{v}_{\rm halo}=230$ km/s. We
incorporate the relative velocity between the Earth and the dark matter, where
$|\bm{v}_{\oplus}|=240\,\mathrm{km/s}$ is the Earth’s velocity in the Galactic
rest frame. Above $v_{\rm esc}^{\rm halo}\approx 500$ km/s is the Galactic
escape velocity at the Sun’s position.

The local dark matter density is assumed to be the standard value of
$\rho_{\chi}=0.4\,\mathrm{GeV/cm^{3}}$ [44], unless otherwise stated in the
case of studying sub-components (as discussed at the end of this section).

###  I.2 Thermalized Dark Matter Population Bound to the Earth

Another dark matter component other than the usually tested halo dark matter
component, is the potential thermalized dark matter population bound to the
Earth. As the halo dark matter component enters the Earth and scatters, it may
continue to scatter and become trapped within the Earth if it loses sufficient
kinetic energy to be below the escape velocity of the Earth. For sufficiently
large dark matter-SM scattering cross sections, the dark matter can be trapped
very quickly after rapidly thermalizing and entering local thermal equilibrium
with the surrounding Earth matter. Under these conditions, the radial
distribution of dark matter within the Earth, $n_{\chi}$, is described by the
first-order differential equation [45]

$\displaystyle\frac{\nabla
n_{\chi}}{n_{\chi}}+\left(\kappa+1\right)\frac{\nabla
T}{T}+\frac{m_{\chi}g}{T}=\frac{\Phi}{n_{\chi}D_{\chi
N}}\frac{R_{\oplus}^{2}}{r^{2}}\,,$ (2)

where $T$ is the Earth’s radial temperature profile at position $r$,
$R_{\oplus}$ is Earth’s radius, $g$ is gravitational acceleration, $D_{\chi
N}\sim\lambda\,v_{\rm th}$ and $\kappa\sim-1/[2(1+m_{\chi}/m_{\rm SM})^{3/2}]$
are diffusion coefficients [45], with $\lambda$ the dark matter mean free
path, $v_{\rm th}$ the dark matter thermal velocity, and $m_{\chi}$ and
$m_{\rm SM}$ the dark matter mass and SM target mass respectively. $\Phi$ is
the incoming captured flux of dark matter particles from the Galactic halo,
and is given by [45]

$\Phi=v_{\chi}\sqrt{\frac{8}{3\pi}}\left[1+\frac{3}{2}\left(\frac{v_{\rm
esc}}{v_{\chi}}\right)^{2}\right]\frac{\rho_{\chi}f_{\rm cap}}{m_{\chi}},$ (3)

where $f_{\rm cap}$ is the fraction of dark matter particles captured that
pass through the Earth.

The dark matter density profile is normalized by enforcing that its volume
integral equals the total number of particles expected within the Earth, i.e.
[45]

$\displaystyle 4\pi\int_{0}^{R}r^{2}\,n_{\chi}\,dr=N_{\chi}\,,$ (4)

where the total number of dark matter particles is given by the product of the
dark matter capture rate $C$ and the lifetime of the Earth $\tau_{\rm
Earth}\sim 4.5$ Gyr,

$N_{\chi}=C\,\tau_{\rm Earth}\,.$ (5)

To obtain the capture rate $C$ for a range of dark matter masses and
interaction cross sections, we use the public package Asteria [46]. Asteria
takes into account multiple kinematic and interaction regimes of interest,
including reflection of strongly interacting light dark matter, and the
relation between the fraction of dark matter particles that pass through the
Earth that are captured, $f_{\rm cap}$, and the dark matter-nucleon scattering
cross section, $\sigma_{\chi N}$.

As was shown in Ref. [45], solving Eq. (2) for $n_{\chi}(r)$ reveals that the
thermalized dark matter population has a large enhancement at the surface of
the Earth, compared to the halo density value which locally for a $e.g.$ 1 GeV
dark matter particle is $\sim 0.4$ cm-3 (see also the earlier Refs. [31, 32,
33, 34, 36, 37, 38]). In contrast, depending on the scattering cross section,
the local dark matter density for a 1 GeV dark matter particle can be about
$\sim 10^{14}\leavevmode\nobreak\ $cm-3 for the population bound to the Earth.
This large enhancement in density is however somewhat offset by the fact that
as this dark matter population is thermalized within the Earth, its velocity
is much lower than the halo velocity. We model the thermalized dark matter
velocity distribution as a truncated Maxwell-Boltzmann distribution,

$\displaystyle f_{\chi}^{\rm
bound}(\bm{v})=\frac{1}{N_{0}}e^{-(\bm{v}/\bm{v}_{\rm
th})^{2}}\Theta(v_{\mathrm{esc}}^{\rm Earth}-v)\,,$ (6)

where $N_{0}$ normalizes the distribution, $v_{\rm th}^{2}=8T_{\chi}/\pi
m_{\chi}$ with $T_{\chi}\simeq 300\,\mathrm{K}$, and Earth’s escape velocity
is $v_{\mathrm{esc}}^{\rm Earth}\approx 11$ km/s. Dark matter being at room
temperature of about $T_{\chi}\simeq 300\,\mathrm{K}$ is expected as dark
matter is not expected to thermalize with the cold device itself, as its mean
free path is much larger than the size of our device of interest. The velocity
in Eq. (6) requires energy detection thresholds of less than about 0.05 eV,
which is not achievable with standard direct detection. As we pointed out in
Ref. [29], this population is already highly detectable with low-threshold
quantum devices. Our transmon qubit device will provide the best existing
sensitivity even for sub-component dark matter.

An important limitation for detection of the Earth-bound thermalized dark
matter distribution is a process called evaporation. In particular, when dark
matter is captured inside the Earth, it receives thermal kicks from the
Earth’s internal temperature. If these kicks cause dark matter to gain so much
energy that it can overcome the escape velocity of the Earth, it leaves the
Earth, or “evaporates”, such that there is no population remaining inside the
Earth to be detected. We take the evaporation process into account by
demanding that the dark matter number density is not depleted, such that total
number of dark matter particles is [47]

$\displaystyle N_{\chi}^{\rm tot}=\frac{C}{E}\,\left(1-e^{-E\,\tau_{\rm
Earth}}\right),$ (7)

where $C$ is the dark matter capture rate determined in Asteria [46] as
discussed above, and $E$ is the evaporation rate which is calculated following
Ref. [47]. Here we will focus on contact interactions only, and on the case
where there is no significant dark matter annihilation.

###  I.3 Testing Dark Matter Sub-Components

As well as studying dark matter which makes up all of the relic abundance of
dark matter, we will also consider dark matter sub-components. Dark matter
sub-components can be defined as some dark matter candidates which do not make
up full the relic abundance of dark matter observed in the Universe today.
That is, the candidate’s contribution to the dark matter density can be
expressed as some fraction $f_{\chi}$ of the usual
$\rho_{\chi}=0.4\,\mathrm{GeV\,cm^{-3}}$ observed at the local position, such
that the local mass density of the candidate will instead be the lower value
of

$\rho_{\chi}^{\rm sub}=f_{\chi}\times 0.4\,\mathrm{GeV\,cm^{-3}}.$ (8)

Sub-component dark matter can be motivated from the perspective that the dark
sector may be vast, with more than one particle candidate, not dissimilar to
the spectrum of particles observed within the SM. There are a range of
particle dark matter models considered in the literature that naturally can
provide dark matter candidates in our parameter space of interest, including
dark matter candidates with quantum chromodynamics (QCD) charge such as the
sexaquark [48, 49, 50, 51, 43]. Another example explored in $e.g.$ Ref. [37]
is the scenario where asymmetric dark matter is coupled to the dark photon,
with predictions in previously untested parameter space to which we are
sensitive. In general, sub-fractions of the full dark matter density may be
expected when the interaction rate is fairly large, leading to large
annihilation rates in the early Universe and therefore a depletion in the
relic dark matter abundance.

Detecting any potential sub-components of dark matter can be difficult. This
is because testing only a very small fraction of the dark matter density
requires higher-precision experiments. Proposals to detect sub-component dark
matter have become of increasing interest recently [31, 32, 33, 34, 35, 36,
37, 38, 39, 40, 41, 42, 43]. To probe sub-component dark matter, we will
consider values in Eq. (8) of $f_{\chi}=1,10^{-3},10^{-6},10^{-9},10^{-12}$,
to span an extreme range of sub-densities, and also to facilitate direct
comparison with previous sub-component dark matter studies, which often have
used these benchmarks in their projections.

Velocities of dark matter sub-components are the same for the halo or
thermalized dark matter cases as relevant as discussed above in the previous
subsections; only the density value is affected by being a sub-component.

##  II Dark Matter Detection with Transmon Qubits

###  II.1 Transmon Qubit Overview

A transmon qubit is a type of superconducting qubit used in quantum computing.
It is a modified version of the Cooper pair box qubit, which consists of a
superconducting island connected to a superconducting reservoir by a Josephson
junction. The transmon qubit is designed to have reduced sensitivity to charge
noise compared to the original Cooper pair box qubit. The key feature of a
transmon qubit is that it has a large shunting capacitor connected in parallel
with the Josephson junction. This capacitor effectively lowers the charging
energy of the qubit, reducing its sensitivity to charge noise. As a result,
transmon qubits can have longer coherence times, which is essential for
performing quantum computations. Transmon qubits are widely used in various
quantum computing architectures, including those based on circuit quantum
electrodynamics and superconducting quantum processors. They have been
instrumental in many experimental demonstrations of quantum algorithms and
quantum error correction codes. Their low-sensitivity to charge noise, and
therefore low backgrounds for the purposes of dark matter detection, makes
them ideal low-threshold dark matter detectors.

In order to investigate how quasiparticle tunneling impacts the decoherence of
a transmon qubit, a study by Ref. [30] devised a single-junction
superconducting qubit using aluminium and examined its decoherence by tracking
the rate of single-charge tunneling. By observing the relaxation rate of the
qubit, they determined a quasiparticle density of $0.04\pm
0.01\,\mu\mathrm{m}^{-3}$, characterized by a thermalized distribution [30].
We previously used this device to probe dark matter scattering in Ref. [29].
Here we extend our investigations by testing the sensitivity of this device to
subcomponents of dark matter which may interact with the device, as well as
mapping out future sensitivities. We now discuss dark matter scattering rates
within these devices, and how dark matter can induce quasiparticle production
and tunneling that is limited by this quasiparticle density measurement.

###  II.2 Calculation of Scattering Rates and Energy Deposition

We are interested in very low energy depositions as appropriate for light dark
matter, or dark matter with low velocities. Therefore, for our regime of
interest interatomic forces will be important, as the momentum transferred by
dark matter $\bm{q}$ will be comparable to the inverse size of the nuclear
wavefunction in the detector crystal. In this regime, lattice vibrations or
phonon excitations will provide the dark matter scattering rate, rather than
the higher energy nuclear recoil regime. For our device the total scattering
rate per unit target mass is [52, 53]

$\displaystyle\Gamma=\frac{\pi\sigma_{\chi N}n_{\chi}}{\rho_{T}\mu^{2}}\int
d^{3}vf_{\chi}(\bm{v})\int\frac{d^{3}q}{(2\pi)^{3}}F_{\mathrm{med}}^{2}(q)S(\bm{q},\omega_{\bm{q}}).$
(9)

Here, $f_{\chi}(\bm{v})$ is dark matter velocity distribution of interest as
detailed in the previous section, $\rho_{T}$ is the target density,
$\sigma_{\chi N}$ is the dark matter-nucleon scattering cross section, $\mu$
is the reduced mass of the dark matter-nucleon system, $F_{\mathrm{med}}(q)$
is a form-factor that depends on the mediator (we assume
$F_{\mathrm{med}}(q)=1$). The differential rate as a function of deposited
energy $\omega$ can be written by inserting a delta function,

$\displaystyle\frac{d\Gamma}{d\omega}=\frac{\pi\sigma_{\chi
N}n_{\chi}}{\rho_{T}\mu^{2}}$ $\displaystyle\int d^{3}v\,f_{\chi}(\bm{v})$
(10)
$\displaystyle\times\int\frac{d^{3}q}{(2\pi)^{3}}\,F_{\mathrm{med}}^{2}(q)\,S(\bm{q},\omega)\,\delta(\omega-\omega_{\bm{q}})\,,$

where the dark matter energy deposition is given by

$\omega_{\bm{q}}=\bm{q}\cdot\bm{v}-\frac{q^{2}}{2m_{\chi}}.$ (11)

In Eq. (9), $S(\bm{q},\omega_{\bm{q}})$ is the dynamic structure factor
containing the detector response to dark matter scattering and depends on the
crystal structure of the target material, which is given by

$\displaystyle S(\bm{q},\omega_{\bm{q}})$
$\displaystyle\approx\dfrac{2\pi}{V_{c}}\sum_{d}f_{d}^{2}\,e^{-2W_{d}(\bm{q})}\sum_{n}\left(\dfrac{q^{2}}{2m_{d}}\right)^{n}\dfrac{1}{n!}$
(12) $\displaystyle\times\left(\prod_{i=1}^{n}\int
d\omega_{i}\dfrac{D_{d}(\omega_{i})}{\omega_{i}}\right)\delta\left(\omega-\sum_{i}\omega_{i}\right)\,,$

Figure 1: Phonon density of state $D(\omega)$ for our target aluminium as a
function of energy $\omega$.

in the large momentum transfer regime [54]. Here, $V_{c}$ is the primitive
cell volume, $n$ is the number of phonons excited, $m_{d}$ is the mass of the
target atoms. The average phonon energy $\bar{\omega}_{d}$, determined by the
density of state of the material, determines the typical number of phonons
excited,

$n\approx\frac{q^{2}}{2m_{d}\bar{\omega}_{d}}.$ (13)

Here, as we have an aluminium target, $\bar{\omega}_{d}=25.5$ meV. As an
example benchmark, for thermalized dark matter with 1 GeV mass and a maximum
momentum transfer of $q=20$ keV, the typical phonon number is therefore
$n=0.29$.

Figure 1 shows $D(\omega)$, the phonon density of state of the single atom in
the primitive cell, for our target aluminium. The phonon density of state is
important for understanding at what energy resonant energy transfer occurs,
and therefore at what energies are there especially large scattering rates.
The peaks shown in the aluminium phonon density of state are from the
transverse and longitudinal acoustic phonon branches. Aluminium has a face-
centered-cubic crystal structure with only one atom in its primitive cell,
which implies that it only has three acoustic phonon branches.

In Eq. (12), we take $f_{d}=A_{d}$ as the coupling for spin-independent
interactions, where the scattering benefits from nuclear coherence. Above the
Debye-Waller function is given by

$\displaystyle W_{d}(\bm{q})=\dfrac{q^{2}}{4m_{d}}\int
d\omega\,\dfrac{D_{d}(\omega)}{\omega}\,.$ (14)

To implement scattering rates we use the publicly available code DarkELF [55,
56], with some additional modifications. To allow for tests of both halo dark
matter and Earth-bound dark matter, we update the local dark matter density
and velocity inputs as detailed in the previous section. We also modify the
code to apply to aluminium, which is our target material. As discussed above,
aluminium has only one atom in the primitive cell. Therefore, we modify
DarkELF to consider only one atom in the primitive cell, rather than for two
atoms per primitive cell as per the default DarkELF setup.

###  II.3 Dark Matter Induced Quasiparticle Density

In superconducting metals, the low temperature allows electrons to form Cooper
pairs. These Cooper pairs are bound via long-range interactions with phonons.
When low-energy dark matter scatters with the superconducting aluminium film
as described in the previous subsection, its kinetic energy can be deposited
in the form of phonons. These phonons can break the Cooper pairs if they have
energy exceeding the Cooper pair binding energy of $\sim$meV. This leads to an
excess of quasiparticles in the metal. These excess quasiparticles may tunnel
out of the superconducting island, and the total number of quasiparticles
present at a given time is sensitive to the rate in which they are produced,
are trapped, or recombine. This can be determined using the mean field results
[57]

$\displaystyle\dfrac{dn_{\mathrm{qp}}}{dt}$
$\displaystyle=-\Gamma_{R}-\Gamma_{T}+\Gamma_{G}$
$\displaystyle\approx-\bar{\Gamma}n_{\mathrm{qp}}^{2}-\bar{\Gamma}_{T}n_{\mathrm{qp}}+\Gamma_{G}\,,$
(15)

where $n_{\mathrm{qp}}$ is the steady-state quasiparticle density, and
$\Gamma_{R},\Gamma_{T},\Gamma_{G}$ are the recombination, trapping, and
generation rates, respectively. The quasiparticle generation rate is related
to the steady state power density injected by dark matter $P_{\rm DM}$ as

$\displaystyle\Gamma_{G}$ $\displaystyle=\frac{\epsilon_{\rm qp}}{2\Delta}\int
d\omega\leavevmode\nobreak\ \omega\,\frac{d\Gamma}{d\omega}$ (16)
$\displaystyle=\frac{P_{\rm DM}}{2\Delta},$ (17)

where $\Delta$ is the superconducting energy gap, $\epsilon_{\rm qp}$ is the
quasiparticle generation efficiency. We assume a quasiparticle generation
efficiency of 60% ($\epsilon_{\rm qp}=0.6$) [58, 14], and an energy gap
$\Delta\simeq 340\,\mathrm{\mu eV}$ for aluminium.

In equilibrium, $i.e.$ when the quasiparticle density is constant in time, Eq.
(15) then simplifies to

$\frac{P_{\rm
DM}}{2\Delta}=\bar{\Gamma}n_{\mathrm{qp}}^{2}+\bar{\Gamma}_{T}n_{\mathrm{qp}}\,.$
(18)

Under the assumption of no quasiparticle trapping, the expression is further
simplified, and the quasiparticle density is then related to the dark matter
power injection as

$n_{\mathrm{qp}}=\sqrt{\frac{P_{\rm DM}}{2\Delta\,\bar{\Gamma}}},$ (19)

where $\bar{\Gamma}=40\,\mathrm{s^{-1}\mu m^{3}}$ for aluminium. The
quasiparticle generation rate $\bar{\Gamma}$ has some associated uncertainty,
and so when making plots we will take $\overline{\Gamma}=4$ or
$400\,\mathrm{s^{-1}}\,\mu\mathrm{m}^{3}$ for conservative and optimistic
values respectively.

Figure 2: Dark matter induced quasiparticle density $n_{\rm qp}$ in aluminium
for three benchmark temperatures $T=50,160,180\,$mK and a benchmark dark
matter-nucleon scattering cross section of $\sigma_{\chi
N}=10^{-30}\,\mathrm{cm^{2}}$. Plot assumes halo dark matter.

The steady-state density induced by dark matter is therefore

$n_{\mathrm{qp}}\approx\left(\frac{P_{\rm DM}}{3.6\times
10^{-21}\mathrm{W}}\right)^{1/2}\,\mu\mathrm{m}^{-3}.$ (20)

For our transmon qubit, this implies that the upper limit of any residual
power injection is $3.92\times 10^{-24}\,\mathrm{W}\mu\mathrm{m}^{-3}$, for
dark matter or any other additional source. Given the scattering formalism
discussed in the previous subsection, this then translates directly into a
constraint on the dark matter scattering rate. Note that by considering this
as an upper limit, we are being conservative. This is because the source of
the quasiparticle density measured is not known, and as we pointed out in Ref.
[29], it may already include a possible positive dark matter signal.

The total quasiparticle density at a temperature $T$ is given by the sum of
the dark matter-induced component and the thermal population,

$\displaystyle
n_{\mathrm{qp}}^{\mathrm{tot}}=n_{\mathrm{qp}}+2\nu_{0}\Delta\sqrt{\dfrac{2\pi
k_{B}T}{\Delta}}\exp\left(-\dfrac{\Delta}{k_{B}T}\right)\,.$ (21)

Here, $\nu_{0}=1.2\times 10^{4}\,\mu\mathrm{m}^{-3}\,\mathrm{\mu eV}^{-1}$ is
the Cooper pair density of state at the Fermi level in aluminium.

Figure 2 shows this quasiparticle density in aluminium for three benchmark
temperatures $T=50,160,180\,$mK and a benchmark dark matter-nucleon scattering
cross section of $\sigma_{\chi N}=10^{-30}\,\mathrm{cm^{2}}$. These choices
are motivated as the 180 mK line is effectively the setup of Ref. [59], except
recombination limited, and 50 mK is chosen as another benchmark as it
corresponds to approaching the limit of what can be realistically achieved.
160 mK is an arbitrary intermediate value. This plot assumes the dark matter
inputs from the incoming Galactic halo dark matter; quasiparticle densities
can be higher for the Earth-bound thermalized dark matter.

###  II.4 Other Single Quasiparticle Devices

In our previous work [29], we considered two other low quasiparticle density
devices other than transmon qubits; a low-noise bolometer [60], and SuperCDMS-
CPD [61], whose power measurements had not previously been considered for a
dark matter search. We now briefly comment on these devices, which can produce
weaker but supporting bounds in addition to our transmon qubit.

Ref. [60] developed a quantum capacitance detector where photon-produced free
electrons in a superconductor tunnel into a small capacitive island. This
setup is embedded in a resonant circuit, and therefore can be referred to as a
“quantum resonator”. This quantum resonator measured excess power of $4\times
10^{-20}\,\mathrm{W}$ [60], making it the most sensitive existing far-infrared
detector. Given the power measurement of this device, it is a few orders of
magnitude less sensitive than our transmon qubit.

SuperCDMS detectors are transition-edge sensors (TES). The lowest bias power
measurement is from SuperCDMS-CPD [61], which is a TES is coupled to a large
silicon absorber [62]. SuperCMDS-CPD measured an excess power of 6 pW was
measured in the phonon sensor arrays, which is an excess substrate power of
$10^{-24}\,\mathrm{W}\mu\mathrm{m}^{-3}$. This therefore produces weaker or
comparable bounds to our transmon qubit, depending on the accuracy of the
quasiparticle generation rate.

Both of these measurements support our constraints, especially SuperCDMS-CPD’s
bias power measurements. Note that for both of these devices, we
conservatively only considered dark matter detection from the superconducting
films, rather than the substrate, see Ref. [29] for extended discussion on
this point.

##  III New Dark Matter Parameter Space Sensitivities

Figure 3: Limit on the dark matter-nucleon scattering cross section (gray-
shaded) for incoming Galactic dark matter, with the two lines being
conservative and optimistic values of the quasiparticle generation rate. We
show two projections (black dashed and dotted) with quasiparticle density
$n_{\mathrm{qp}}=10^{-3}\,\mu\mathrm{m}^{-3}$ for two temperatures as labeled,
see text for motivation.

Figure 3 shows the bounds and sensitivities we derive for spin-independent
scattering using a transmon qubit, using the incoming Galactic dark matter,
under the assumption the dark matter makes up the full relic abundance ($i.e.$
$f_{\chi}=1$). This is in direct one to one correspondence with the usual
direct detection bounds, that is, all lines here have about the same
assumptions of the local dark matter density and velocity. The limits in Fig.
3 can be simply rescaled linearly in $f_{\chi}$ for subcomponent dark matter.
We see that our bounds on halo dark matter are stronger than lab-based
experiments below about 100 MeV, as we already had shown in Ref. [29]. We also
show new projected sensitivities for a few benchmark scenarios, assuming
improved future single quasiparticle devices. We have taken a quasiparticle
density about an order of magnitude than our existing transmon qubit, as
labeled as the dashed line with $n_{\rm qp}=10^{-3}\,\mu\mathrm{m}^{-3}$, at
two quasiparticle temperatures of 180 mK and 50 mK as labeled. As noted in the
previous section, these choices are motivated as the 180 mK line is
effectively the setup of Ref. [59], except recombination limited, and 50 mK is
chosen as another benchmark as it corresponds to approaching the limit of what
can be realistically achieved in the near future. Note that all estimates
assume recombination-limited quasiparticle density modeling. Active work to
understand sources of quasiparticle poisoning for qubit error mitigation, both
from qubit drive [63] and radiation [64, 65], will contribute to even more
stringent sensitivities than the projections shown, as quasiparticle
backgrounds are reduced further.

In Fig. 3 we also compare with existing limits, including those from
astrophysical systems such as Milky Way satellites [66], Lyman-alpha [67], and
the Cosmic Microwave Background (CMB) [68]. There are also lab experiments
overlapping with part of our parameter space, namely CRESST [69], SuperCDMS
[70], Edelweiss [71], CDEX [72], XQC [73], and CRESST UG [74]. Note that for
cross sections above about $10^{-30}$ cm2, the Born approximation breaks down,
and the nuclear coherence across different detector materials is not well
defined without using a particle dark matter model, see $e.g.$ Refs. [75, 76].

Figure 4: New constraints on the dark matter-nucleon spin-independent
scattering cross section $\sigma_{\chi N}^{\rm SI}$ as a function of dark
matter mass $m_{\chi}$, using transmon qubit measurements for the Earth-bound
dark matter population, for varying fractions of the dark matter density
$f_{\chi}$ as labeled. Dashed line is the conservative bound, solid is the
optimistic bound (see text for details). Complementary limits are also shown
as labeled.

Figure 4 shows our new limits on the dark matter-nucleon spin-independent
scattering cross section as a function of dark matter mass, using transmon
qubit measurements for the Earth-bound dark matter population, for varying
fractions of the dark matter density $f_{\chi}$ as labeled. The solid line
corresponds to the optimistic bound based on quasiparticle generation rate
uncertainties, and the dashed line is conservative. We see that while the
$f_{\chi}=1$ case (dark matter candidate corresponds to the full dark matter
density) overlaps most with existing experiments, as soon as this density is
decreased other experiments weaken or disappear, while our bounds are not
largely changed in comparison.

In Fig. 4, note that we linearly rescale the complementary direct detection
limits only from below; the ceilings of the limits also decrease but non-
linearly with decreasing dark matter fraction, and require a dedicated
calculation [77]. The reason the ceilings scale non-linearly with decreasing
dark matter fraction, is that the dark matter velocity is also relevant when
the scattering limit is near the detection threshold [77]. As only CRESST has
had its limit ceilings recast for subfractions [77], we adopt the recast
CRESST bounds for both the lower and upper ends of the bounds, but only
rescale the other direct detection constraints from below. From the point of
view of how much new parameter space our transmon qubit probes, this is a
conservative choice, although the ceiling limit change with decreasing dark
matter subfractions in any case is not very large [77]. In the case where
$f_{\chi}=0.1$, already the astrophysical bounds disappear. This is because
limits relying on structure formation require the dark matter candidate to
make up most of the dark matter content, and as soon as dark matter is a
smaller component, it has a negligible effect on structure.

In Fig 4, we show an extreme range of sub-fractions of dark matter, to
facilitate direct comparison to the same values chosen in the literature, but
also to demonstrate the strong constraining power of transmon qubits. We
produce limits stronger than any other existing study on the Earth-bound
thermalized population, and even produce exclusions which outperform the
projections of most of the proposals to detect this population. We see that
even for subfractions that are as tiny as $10^{-12}$, we still have
sensitivity to dark matter parameters. On the other hand, existing direct
detection experiments completely lose sensitivity for sub-components below
about $10^{-9}$, as shown in Fig. 4.

Notice that all the thermalized dark matter bounds in Fig. 4 are truncated by
the same S-shaped curve. This line is not due to inherent sensitivity in the
transmon qubit, but rather due to evaporation. As discussed earlier, when dark
matter is captured inside the Earth, it receives thermal kicks from the
Earth’s internal temperature. If these kicks cause dark matter to gain so much
energy that it can overcome the escape velocity of the Earth, it leaves the
Earth, such that there is no population remaining inside the Earth to be
detected. This is actually our key limiting factor across most of the
parameter space for the Earth-bound dark matter. This evaporation line is
obtained for non-annihilating dark matter following the original treatment of
Ref. [78] (see Eq. (7)); note that other works studying the Earth-bound
population have often instead used the surface of last scattering, which is
less accurate at lower scattering cross sections [78]. Any difference in the
location of the evaporation-limited truncation of our bounds compared to other
works studying this Earth-bound population should be applied equally; $i.e.$
our transmon qubit sensitivities themselves are not actually weaker, instead
generally our treatment of evaporation is more conservative, and if we adopt
the surface of last scattering treatment as per other Earth-bound works, we
would obtain the same limit truncation as other works.

Figure 5 shows the strength of the transmon qubit constraints on the
thermalized dark matter population (with $f_{\chi}=1$), assuming no
evaporation. Here this underscores the point discussed in the previous
paragraph, that our results are largely evaporation limited, not detector
sensitivity limited. The dashed line for comparison shows where evaporation
rapidly removes the dark matter population from the Earth (under the
assumption of contact interactions); any other treatment of dark matter
evaporation may simply cut our plot with where evaporation becomes relevant
for the given scenario. One scenario in which this is interesting, is where
evaporation is inhibited due to the presence of long-range forces [79].
However it is important to note in that case that long-range forces which
serve to avoid evaporation generally also alter the dark matter distribution
inside the Earth. In some cases this can be significant, as avoiding
evaporation requires significantly reinforcing gravity, and gravity is a key
input for the dark matter distribution. We have not altered the distribution
inside the Earth in this plot, as the precise way this happens will depend on
the particle model. Instead here we aim to make the point that (i) our
detector has strong sensitivity outside the evaporation range, and evaporation
is largely what limits our constraints, as well as (ii) that our transmon
qubit may provide strong sensitivities to particle dark matter models which do
not evaporate at these masses.

Note that in Fig. 5 we have restricted the cross section range to be the same
as the other figures, as this is the regime where the dark matter mean free
path is much shorter than the size of the Earth. The transmon qubit certainly
has sensitivity at even lower cross sections, but for simplicity and the point
we want to make here, we do not study the long-mean free path regime.

Fig. 5 also helps make sense of the shape of the halo dark matter limits in
Fig. 3 compared to the thermalized dark matter limits: the halo limits do not
care about evaporation, because even dark matter which will later evaporate
does first enter the Earth at the detector position, which is the stage in
which the limit is set for regular halo dark matter. The halo limits therefore
extend to lower dark matter masses than the thermalized dark matter limits,
and do not feature the S-like shape at their boundary. It also explains why
the Earth-bound thermalized dark matter limits do not weaken linearly with the
dark matter density sub-fraction value; it is not until the sub-fraction
causes the sensitivity to weaken above the evaporation threshold that the
linear scaling with dark matter density is revealed, as is seen in the bottom
two panels of Fig. 4.

Figure 5: Constraints on the dark matter-nucleon scattering cross section as a
function of dark matter mass, using transmon qubit data assuming that the
thermalized dark matter does not evaporate (cyan shaded, full plot range). The
evaporation line for dark matter in the Earth is shown for comparison as the
dashed line, and complementary constraints are also shown.

##  IV Summary and Outlook

Low-energy dark matter interactions, such as those arising from light dark
matter or dark matter thermalized and bound to the Earth, require new
techniques for detection. In this work we have investigated the sensitivity of
a transmon qubit with world-record low quasiparticle density to low-energy
dark matter interactions. This transmon qubit is advantageous over traditional
direct detection experiments, as its energy deposition detection threshold is
as low as about an meV, six orders of magnitude lower than standard direct
detection setups. The meV threshold is set by the Cooper-pair binding energy
in the superconducting island, as when the Cooper pairs are broken by dark
matter interactions, they release detectable quasiparticles. The transmon
qubit can be therefore substantially more sensitive than direct detection to
low-energy dark matter interactions.

We investigated current bounds and future sensitivities for two low-energy
dark matter components: light dark matter incoming from the Galactic halo, as
well as the dark matter distribution thermalized and bound to the Earth. For
the thermalized population, there have recently been a growing number of
proposals to detect it, especially if dark matter which scatters with the
detector makes up some sub-component of the full dark matter relic density. We
set new constraints on dark matter sub-fractions as low as $10^{-12}$, to
which no existing experiment is currently sensitive. For many other
subfractions, our transmon qubit provides the strongest constraints to date on
thermalized dark matter. For the thermalized population, we also showed how
the sensitivity depends on the evaporation threshold, with the main limitation
of our sensitivity arising from evaporation rather than the transmon qubit’s
sensitivity itself. We showed that these devices may be fruitful moving
forward for detecting dark matter particle models that do not evaporate at
very light dark matter masses.

For the Galactic halo dark matter population, transmon qubits provide
constraints on dark matter-nucleon scattering stronger than direct detection
for dark matter masses below about 100 MeV. We showed the future promise for
transmon qubits or any low-quasiparticle density device to detect light dark
matter, considering both lower quasiparticle densities, as well as devices
cooled to lower temperatures. We expect that low quasiparticle density devices
can reach cross sections as low as about $10^{-31}$ cm2 for dark matter masses
below about 100 MeV for Galactic halo dark matter, and many orders of
magnitude below this for thermalized dark matter. Active work to understand
sources of quasiparticle poisoning for qubit error mitigation, both from qubit
drive [63] and radiation [64, 65] will led to reduced quasiparticle
backgrounds, and even better sensitivities in the future as quasiparticle
backgrounds are reduced further. Overall, new low-quasiparticle density
devices as they are actively created for applications such as quantum
computing, infrared telescopes, or radiation detectors, offer substantial
potential to probe dark matter with low energy deposition, and are a promising
path forward to detect light dark matter.

## Acknowledgments

We thank Chris Cappiello and Tom Connolly for helpful discussions. AD was
supported by Grant Korea NRF-2019R1C1C1010050. NK and RKL are supported by the
U.S. Department of Energy under Contract DE-AC02-76SF00515. NK is also
supported by the US. Department of Energy Early Career Research Program (ECRP)
under FWP 100872.

## References

  * Aalbers  _et al._ [2022] J. Aalbers _et al._ (LUX-ZEPLIN), (2022), arXiv:2207.03764 [hep-ex] .
  * Vergados and Ejiri [2005] J. D. Vergados and H. Ejiri, Phys. Lett. B 606, 313 (2005), arXiv:hep-ph/0401151 .
  * Moustakidis  _et al._ [2005] C. C. Moustakidis, J. D. Vergados, and H. Ejiri, Nucl. Phys. B 727, 406 (2005), arXiv:hep-ph/0507123 .
  * Bernabei  _et al._ [2007] R. Bernabei _et al._ , Int. J. Mod. Phys. A 22, 3155 (2007), arXiv:0706.1421 [astro-ph] .
  * Ibe  _et al._ [2018] M. Ibe, W. Nakano, Y. Shoji, and K. Suzuki, JHEP 03, 194, arXiv:1707.07258 [hep-ph] .
  * Aprile  _et al._ [2019a] E. Aprile _et al._ (XENON), Phys. Rev. Lett. 123, 241803 (2019a), arXiv:1907.12771 [hep-ex] .
  * Tiffenberg  _et al._ [2017] J. Tiffenberg, M. Sofo-Haro, A. Drlica-Wagner, R. Essig, Y. Guardincerri, S. Holland, T. Volansky, and T.-T. Yu (SENSEI), Phys. Rev. Lett. 119, 131802 (2017), arXiv:1706.00028 [physics.ins-det] .
  * Barak  _et al._ [2020] L. Barak _et al._ (SENSEI), Phys. Rev. Lett. 125, 171802 (2020), arXiv:2004.11378 [astro-ph.CO] .
  * Aprile  _et al._ [2019b] E. Aprile _et al._ (XENON), Phys. Rev. Lett. 123, 251801 (2019b), arXiv:1907.11485 [hep-ex] .
  * Hochberg  _et al._ [2022] Y. Hochberg, Y. F. Kahn, R. K. Leane, S. Rajendran, K. Van Tilburg, T.-T. Yu, and K. M. Zurek, Nature Rev. Phys. 4, 637 (2022).
  * Hochberg  _et al._ [2016a] Y. Hochberg, Y. Zhao, and K. M. Zurek, Phys. Rev. Lett. 116, 011301 (2016a), arXiv:1504.07237 [hep-ph] .
  * Hochberg  _et al._ [2016b] Y. Hochberg, M. Pyle, Y. Zhao, and K. M. Zurek, JHEP 08, 057, arXiv:1512.04533 [hep-ph] .
  * Hochberg  _et al._ [2021a] Y. Hochberg, Y. Kahn, N. Kurinsky, B. V. Lehmann, T. C. Yu, and K. K. Berggren, Phys. Rev. Lett. 127, 151802 (2021a), arXiv:2101.08263 [hep-ph] .
  * Hochberg  _et al._ [2021b] Y. Hochberg, E. D. Kramer, N. Kurinsky, and B. V. Lehmann, (2021b), arXiv:2109.04473 [hep-ph] .
  * Hochberg  _et al._ [2019] Y. Hochberg, I. Charaev, S.-W. Nam, V. Verma, M. Colangelo, and K. K. Berggren, Phys. Rev. Lett. 123, 151802 (2019), arXiv:1903.05101 [hep-ph] .
  * Chiles  _et al._ [2022] J. Chiles _et al._ , Phys. Rev. Lett. 128, 231802 (2022), arXiv:2110.01582 [hep-ex] .
  * Schutz and Zurek [2016] K. Schutz and K. M. Zurek, Phys. Rev. Lett. 117, 121302 (2016), arXiv:1604.08206 [hep-ph] .
  * Knapen  _et al._ [2017] S. Knapen, T. Lin, and K. M. Zurek, Phys. Rev. D 95, 056019 (2017), arXiv:1611.06228 [hep-ph] .
  * Caputo  _et al._ [2019] A. Caputo, A. Esposito, and A. D. Polosa, Phys. Rev. D 100, 116007 (2019), arXiv:1907.10635 [hep-ph] .
  * Griffin  _et al._ [2018] S. Griffin, S. Knapen, T. Lin, and K. M. Zurek, Phys. Rev. D 98, 115034 (2018), arXiv:1807.10291 [hep-ph] .
  * Knapen  _et al._ [2018] S. Knapen, T. Lin, M. Pyle, and K. M. Zurek, Phys. Lett. B 785, 386 (2018), arXiv:1712.06598 [hep-ph] .
  * Cox  _et al._ [2019] P. Cox, T. Melia, and S. Rajendran, Phys. Rev. D 100, 055011 (2019), arXiv:1905.05575 [hep-ph] .
  * Sánchez-Martínez  _et al._ [2019] M.-A. Sánchez-Martínez, I. n. Robredo, A. Bidauzarraga, A. Bergara, F. de Juan, A. G. Grushin, and M. G. Vergniory, Materials 3, 014001 (2019), arXiv:1905.04805 [cond-mat.mtrl-sci] .
  * Hochberg  _et al._ [2018] Y. Hochberg, Y. Kahn, M. Lisanti, K. M. Zurek, A. G. Grushin, R. Ilan, S. M. Griffin, Z.-F. Liu, S. F. Weber, and J. B. Neaton, Phys. Rev. D 97, 015004 (2018), arXiv:1708.08929 [hep-ph] .
  * Geilhufe  _et al._ [2018] R. M. Geilhufe, B. Olsthoorn, A. Ferella, T. Koski, F. Kahlhoefer, J. Conrad, and A. V. Balatsky, Phys. Status Solidi RRL 12, 1800293 (2018), arXiv:1806.06040 [cond-mat.mtrl-sci] .
  * Geilhufe  _et al._ [2020] R. M. Geilhufe, F. Kahlhoefer, and M. W. Winkler, Phys. Rev. D 101, 055005 (2020), arXiv:1910.02091 [hep-ph] .
  * Coskuner  _et al._ [2021] A. Coskuner, A. Mitridate, A. Olivares, and K. M. Zurek, Phys. Rev. D 103, 016006 (2021), arXiv:1909.09170 [hep-ph] .
  * Das  _et al._ [2023] A. Das, J. Jang, and H. Min, (2023), arXiv:2312.00866 [hep-ph] .
  * Das  _et al._ [2024] A. Das, N. Kurinsky, and R. K. Leane, Phys. Rev. Lett. 132, 121801 (2024), arXiv:2210.09313 [hep-ph] .
  * Ristè  _et al._ [2013] D. Ristè, C. C. Bultink, M. J. Tiggelman, R. N. Schouten, K. W. Lehnert, and L. Dicarlo, Nature Communications 4, 1913 (2013), arXiv:1212.5459 [cond-mat.mes-hall] .
  * Neufeld  _et al._ [2018] D. A. Neufeld, G. R. Farrar, and C. F. McKee, Astrophys. J. 866, 111 (2018), arXiv:1805.08794 [astro-ph.CO] .
  * Pospelov and Ramani [2021] M. Pospelov and H. Ramani, Phys. Rev. D 103, 115031 (2021), arXiv:2012.03957 [hep-ph] .
  * Pospelov  _et al._ [2020] M. Pospelov, S. Rajendran, and H. Ramani, Phys. Rev. D 101, 055001 (2020), arXiv:1907.00011 [hep-ph] .
  * Rajendran and Ramani [2021] S. Rajendran and H. Ramani, Phys. Rev. D 103, 035014 (2021), arXiv:2008.06061 [hep-ph] .
  * Xu and Farrar [2021] X. Xu and G. R. Farrar, (2021), arXiv:2112.00707 [hep-ph] .
  * Budker  _et al._ [2022] D. Budker, P. W. Graham, H. Ramani, F. Schmidt-Kaler, C. Smorra, and S. Ulmer, PRX Quantum 3, 010330 (2022), arXiv:2108.05283 [hep-ph] .
  * McKeen  _et al._ [2022] D. McKeen, M. Moore, D. E. Morrissey, M. Pospelov, and H. Ramani, (2022), arXiv:2202.08840 [hep-ph] .
  * Billard  _et al._ [2022] J. Billard, M. Pyle, S. Rajendran, and H. Ramani, (2022), arXiv:2208.05485 [hep-ph] .
  * Li  _et al._ [2023] Y. Li, Z. Liu, and Y. Xue, JCAP 05, 060, arXiv:2209.04387 [hep-ph] .
  * Pospelov and Ray [2024] M. Pospelov and A. Ray, JCAP 2024 (01), 029, arXiv:2309.10032 [hep-ph] .
  * McKeen  _et al._ [2023] D. McKeen, D. E. Morrissey, M. Pospelov, H. Ramani, and A. Ray, Phys. Rev. Lett. 131, 011005 (2023), arXiv:2303.03416 [hep-ph] .
  * Ema  _et al._ [2024] Y. Ema, M. Pospelov, and A. Ray, (2024), arXiv:2402.03431 [hep-ph] .
  * Moore and Slatyer [2024] M. Moore and T. R. Slatyer, (2024), arXiv:2403.03972 [hep-ph] .
  * Pato  _et al._ [2015] M. Pato, F. Iocco, and G. Bertone, Journal of Cosmology and Astroparticle Physics 2015 (12), 001–001.
  * Leane and Smirnov [2023a] R. K. Leane and J. Smirnov, JCAP 10, 057, arXiv:2209.09834 [hep-ph] .
  * Leane and Smirnov [2023b] R. K. Leane and J. Smirnov, JCAP 12, 040, arXiv:2309.00669 [hep-ph] .
  * Gould [1987] A. Gould, Astrophys. J. 321, 560 (1987).
  * Jaffe [1977] R. L. Jaffe, Phys. Rev. Lett. 38, 195 (1977).
  * Farrar and Zaharijas [2003] G. R. Farrar and G. Zaharijas, Phys. Lett. B 559, 223 (2003), [Erratum: Phys.Lett.B 575, 358–358 (2003)], arXiv:hep-ph/0302190 .
  * Farrar and Zaharijas [2006] G. R. Farrar and G. Zaharijas, Phys. Rev. Lett. 96, 041302 (2006), arXiv:hep-ph/0510079 .
  * Farrar [2017] G. R. Farrar, (2017), arXiv:1708.08951 [hep-ph] .
  * Kahn and Lin [2022] Y. Kahn and T. Lin, Rept. Prog. Phys. 85, 066901 (2022), arXiv:2108.03239 [hep-ph] .
  * Trickle  _et al._ [2020] T. Trickle, Z. Zhang, K. M. Zurek, K. Inzani, and S. M. Griffin, JHEP 03, 036, arXiv:1910.08092 [hep-ph] .
  * Campbell-Deem  _et al._ [2020] B. Campbell-Deem, P. Cox, S. Knapen, T. Lin, and T. Melia, Phys. Rev. D 101, 036006 (2020), [Erratum: Phys.Rev.D 102, 019904 (2020)], arXiv:1911.03482 [hep-ph] .
  * Knapen  _et al._ [2022] S. Knapen, J. Kozaczuk, and T. Lin, Phys. Rev. D 105, 015014 (2022), arXiv:2104.12786 [hep-ph] .
  * Campbell-Deem  _et al._ [2022] B. Campbell-Deem, S. Knapen, T. Lin, and E. Villarama, (2022), arXiv:2205.02250 [hep-ph] .
  * Bespalov  _et al._ [2016] A. Bespalov, M. Houzet, J. S. Meyer, and Y. V. Nazarov, Physical Review Letters 117, 10.1103/physrevlett.117.117002 (2016).
  * Kaplan  _et al._ [1976] S. B. Kaplan, C. C. Chi, D. N. Langenberg, J. J. Chang, S. Jafarey, and D. J. Scalapino, Phys. Rev. B 14, 4854 (1976).
  * Connolly  _et al._ [2023] T. Connolly, P. D. Kurilovich, S. Diamond, H. Nho, C. G. L. Bøttcher, L. I. Glazman, V. Fatemi, and M. H. Devoret, (2023), arXiv:2302.12330 [quant-ph] .
  * Echternach  _et al._ [2018] P. M. Echternach, B. J. Pepper, T. Reck, and C. M. Bradford, Nature Astronomy 2, 90 (2018).
  * Fink  _et al._ [2021] C. W. Fink, S. L. Watkins, T. Aramaki, P. L. Brink, J. Camilleri, X. Defay, S. Ganjam, Y. G. Kolomensky, R. Mahapatra, N. Mirabolfathi, W. A. Page, R. Partridge, M. Platt, M. Pyle, B. Sadoulet, B. Serfass, and S. Z. and, Applied Physics Letters 118, 022601 (2021).
  * Ren  _et al._ [2021] R. Ren, C. Bathurst, Y. Chang, R. Chen, C. Fink, Z. Hong, N. Kurinsky, N. Mast, N. Mishra, V. Novati, G. Spahn, H. M. zu Theenhausen, S. Watkins, Z. Williams, M. Wilson, A. Zaytsev, D. Bauer, R. Bunker, E. Figueroa-Feliciano, M. Hollister, L. Hsu, P. Lukens, R. Mahapatra, N. Mirabolfathi, B. Nebolsky, M. Platt, F. Ponce, M. Pyle, T. Reynolds, and T. Saab, Physical Review D 104, 10.1103/physrevd.104.032010 (2021).
  * Amin  _et al._ [2024] K. R. Amin, A. M. Eriksson, M. Kervinen, L. Andersson, R. Rehammar, and S. Gasparinetti, (2024), arXiv:2404.01277 [quant-ph] .
  * Fink  _et al._ [2023] C. W. Fink, C. P. Salemi, B. A. Young, D. I. Schuster, and N. A. Kurinsky, (2023), arXiv:2310.01345 [physics.ins-det] .
  * Yelton  _et al._ [2024] E. Yelton _et al._ , (2024), arXiv:2402.15471 [quant-ph] .
  * Nadler  _et al._ [2021] E. O. Nadler _et al._ (DES), Phys. Rev. Lett. 126, 091101 (2021), arXiv:2008.00022 [astro-ph.CO] .
  * Rogers  _et al._ [2022] K. K. Rogers, C. Dvorkin, and H. V. Peiris, Phys. Rev. Lett. 128, 171301 (2022), arXiv:2111.10386 [astro-ph.CO] .
  * Xu  _et al._ [2018] W. L. Xu, C. Dvorkin, and A. Chael, Phys. Rev. D 97, 103530 (2018), arXiv:1802.06788 [astro-ph.CO] .
  * Angloher  _et al._ [2017] G. Angloher _et al._ (CRESST), Eur. Phys. J. C 77, 637 (2017), arXiv:1707.06749 [astro-ph.CO] .
  * Alkhatib  _et al._ [2021] I. Alkhatib _et al._ (SuperCDMS), Phys. Rev. Lett. 127, 061801 (2021), arXiv:2007.14289 [hep-ex] .
  * Armengaud  _et al._ [2019] E. Armengaud _et al._ (EDELWEISS), Phys. Rev. D 99, 082003 (2019), arXiv:1901.03588 [astro-ph.GA] .
  * Liu  _et al._ [2022] Z. Z. Liu _et al._ (CDEX), Phys. Rev. D 105, 052005 (2022), arXiv:2111.11243 [hep-ex] .
  * Mahdawi and Farrar [2018] M. S. Mahdawi and G. R. Farrar, JCAP 10, 007, arXiv:1804.03073 [hep-ph] .
  * Petricca  _et al._ [2020] F. Petricca _et al._ (CRESST), J. Phys. Conf. Ser. 1342, 012076 (2020), arXiv:1711.07692 [astro-ph.CO] .
  * Digman  _et al._ [2019] M. C. Digman, C. V. Cappiello, J. F. Beacom, C. M. Hirata, and A. H. G. Peter, Phys. Rev. D 100, 063013 (2019), arXiv:1907.10618 [hep-ph] .
  * Xu and Farrar [2020] X. Xu and G. R. Farrar, (2020), arXiv:2101.00142 [hep-ph] .
  * Cappiello [2023] C. V. Cappiello, Phys. Rev. Lett. 130, 221001 (2023), arXiv:2301.07728 [hep-ph] .
  * Gould [1990] A. Gould, Astrophys. J. 356, 302 (1990).
  * Acevedo  _et al._ [2024] J. F. Acevedo, R. K. Leane, and J. Smirnov, JCAP 04, 038, arXiv:2303.01516 [hep-ph] .