$(\mathbf{x},\mathbf{p})\leftarrow\mathsf{Algorithm~{}\ref{alg:R1}(\mathbf{x},\mathbf{p})}$

7: else if $|\mathbf{x}_{a}|<|\mathbf{x}_{c}|$ then

8: $k_{c}\leftarrow j\in\mathbf{x}_{c}\cap K$

9: $\mathsf{Transfer}(k_{c},a)$

10: else if $|\mathbf{x}_{b}|<|\mathbf{x}_{a}|$ or $|\mathbf{x}_{b_{1}}|\geq
3$ then

11: if $\exists j_{a}\in\mathbf{x}_{a}\text{ s.t. }j_{a}\in\mathsf{MPB}_{b}$
then

12: $\mathsf{Transfer}(j_{a},b)$

13: else

14: raise payments of $\mathbf{x}_{a}$ by a factor of $k$

15: $j_{a}\leftarrow j\in\mathbf{x}_{a}$

16: $\mathsf{Transfer}(j_{a},b)$

17: else

18: if $\exists j_{c}\in\mathbf{x}_{c}\text{ s.t. }d_{b}(j_{c})=1$ then

19: $\mathsf{Transfer}(j_{c},b)$

20: else

21: $k_{c}\leftarrow j\in\mathbf{x}_{c}\cap K$

22: $j_{b}\leftarrow\\{j\in\mathbf{x}_{b}\mid d_{c}(j)=1$ }

23: $\mathsf{Swap}(k_{c},j_{b})$

24:return $(\mathbf{x},\mathbf{p})$

We show that:

######  Lemma 26.

Algorithm 8 returns an EFX competitive equilibrium.

###### Proof.

We first consider some properties of any input of Algorithm 8. When $R=2$ and
$|N_{1}|=1$, $a$ cannot EFX-envy $b$ or $c$, as they will always have at least
as many $k$-chores as $a$, and $b$ and $c$ cannot EFX-envy $a$ as they see all
of her chores as $k$. Furthermore, without loss of generality we may assume
that $c$ receives her first $k$-chore before $b$ does, so $c$ will always have
at least as many $k$-chores as $b$. Thus, the EFX-envy must exist from $c$
towards $b$. We now cover the conditions in Lines 5, 7, 10, and 17.

In Line 5, when $c$ has exactly as many 1-chores as $a$, we have that all
agents have equal number 1-chores and balanced $k$-chores. Then, we may use
Algorithm 5 to find an EFX allocation.

In Line 7, we now suppose that $c$ has fewer 1-chores than $a$, and
$|\mathbf{x}_{a}|<|\mathbf{x}_{c}|$. Note that the latter condition only holds
when no transfers occur in Algorithm 4. In this case, we transfer a universal
$k$-chore from $c$ to $a$. The number of chores remains balanced between all
agents. Now $b$ and $c$ possess the same number of $k$-chores, so they cannot
EFX-envy each other, and it is clear that they remain unenvious of $a$. We
need only show that $a$ does not begin to EFX-envy either $b$ or $c$. Since
$c$ has fewer 1-chores than $a$, yet $|\mathbf{x}_{a}|<|\mathbf{x}_{c}|$, it
must be that $c$ has at least two more $k$-chores than $a$ and $b$ has at
least one more $k$-chore than $a$. Then, even after $c$ loses a $k$-chore to
$a$, $a$ has at most as many $k$-chores as $b$ and $c$, so she does not EFX-
envy them.

Now in Line 10, we first consider the case when
$|\mathbf{x}_{b}|<|\mathbf{x}_{a}|$. We transfer a chore from $a$ to $b$,
which both $b$ and $c$ value as $k$. The number of chores remains balanced
between all agents, and $b$ and $c$ now have the same number of $k$-chores.
Since they still see all of $a$’s chores as $k$, they remain unenvious of $a$,
and $a$ herself clearly remains unenvious as she has only lost a chore. Thus,
the allocation is EFX. If it is not the case that
$|\mathbf{x}_{b}|<|\mathbf{x}_{a}|$, we have
$|\mathbf{x}_{c}|\leq|\mathbf{x}_{a}|\leq|\mathbf{x}_{b}|$. Since
$|\mathbf{x}_{c}|<|\mathbf{x}_{b}|$ implies that $b$ has as many $k$-chores as
$c$, and thus no EFX-envy between them, we must in fact have
$|\mathbf{x}_{c}|=|\mathbf{x}_{a}|=|\mathbf{x}_{b}|=p$. From
$|\mathbf{x}_{c}|=|\mathbf{x}_{b}|$ and the fact that $c$ EFX-envies $b$ we
also conclude that $b$ has one more 1-chore than $c$, so it must be that $k>2$
for $c$ to EFX-envy $b$. Now suppose $b$ has at least three 1-chores and she
receives a chore from $a$ with value $k$. Clearly $a$ remains unenvious of $b$
and $c$. By Lemma 19 $b$ and $c$ do not EFX-envy each other, as they have
equal number $k$-chores and balanced total chores. In addition, $c$ cannot
EFX-envy $a$ as the two still have balanced total chores and $c$ sees all of
$a$’s chores as $k$. We need only check that $b$ does not begin to EFX-envy
$a$. Since $b$ has at least three 1-chores and at most $p-2$ $k$-chores, and
we know $k>2$, we have
$d_{b}(\mathbf{x}_{b})-1\leq(p-2)k+2\leq(p-1)k=d_{b}(\mathbf{x}_{a})$. So $b$
does not EFX-envy $a$ and the allocation is EFX.

In Line 17, as none of the previous conditions hold, we must have
$|\mathbf{x}_{c}|=|\mathbf{x}_{a}|=|\mathbf{x}_{b}|=p$,
$|\mathbf{x}_{b}|_{1}<3$, and $k>2$. We argue that in fact
$|\mathbf{x}_{b}|_{1}=2$ and $|\mathbf{x}_{c}|_{1}=1$. That is, $b$ must have
exactly two 1-chores and $c$ must have exactly one 1-chore. Note that for $c$
to have EFX-envy for $b$, $c$ must have at least one 1-chore. Then we must
have $|\mathbf{x}_{b}|_{1}\neq|\mathbf{x}_{c}|_{1}$, as otherwise
$|\mathbf{x}_{b}|=|\mathbf{x}_{c}|$ would imply that $b$ and $c$ have the same
number of $k$-chores and no EFX-envy between them. So, we must have
$|\mathbf{x}_{b}|_{1}>|\mathbf{x}_{c}|_{1}\geq 1$ and in conjunction with
$|\mathbf{x}_{b}|_{1}<3$ the desired result follows. Now suppose there exists
a 1-chore of $c$’s which is also value 1 for $b$. Upon transferring this chore
to $b$, we can see that $c$ no longer EFX-envies $b$, as $c$ now only has
$p-1$ $k$-chores, so $d_{c}(\mathbf{x}_{c})-k=(p-2)k<(p-2)k+3\leq
d_{c}(\mathbf{x}_{b})$. We also check that $b$ does not start EFX-envying $c$.
Since $k>2$ we have
$d_{b}(\mathbf{x}_{b})-1=(p-2)k+2\leq(p-1)k=d_{b}(\mathbf{x}_{c})$. It is easy
to verify that $a$ still does not EFX-envy $b$ or $c$, as $a$ still has
balanced total chores with respect to either and fewer $k$-chores than either.
Thus, the allocation is EFX.

Note that all transfers consist of only MPB chores, so it remains that we have
a CE and the allocation is fPO. ∎

###  A.8 Algorithm 9: $\mathbf{R=2,|N_{1}|=2}$

Algorithm 9 When $R=2$ and $|N_{1}|=2$

Input: EF1 competitive equilibrium $(\mathbf{x},\mathbf{p})$
Output: EFX competitive equilibrium $(\mathbf{x},\mathbf{p})$

1:if $\mathbf{x}$ is not EFX then

2: $a\leftarrow\mathsf{argmax}_{i\in N_{1}}|\mathbf{x}_{i_{1}}|$

3: $b\leftarrow i\in N_{1}\setminus\\{a\\}$

4: $c\leftarrow i\in N_{2}$

5: if $|\mathbf{x}_{c}|>|\mathbf{x}_{a}|$ then

6: if $|\mathbf{x}_{c_{k}}|>|\mathbf{x}_{b_{k}}|$ then

7: $k_{c}\leftarrow j\in\mathbf{x}_{c}\cap K$

8: $\mathsf{Transfer}(k_{c},a)$

9: else

10:
$(\mathbf{x},\mathbf{p})\leftarrow\mathsf{Algorithm~{}\ref{alg:R1}(\mathbf{x},\mathbf{p})}$

11: else

12: if $|\mathbf{x}_{c_{1}}|\geq|\mathbf{x}_{b_{1}}|$ then

13:
$(\mathbf{x},\mathbf{p})\leftarrow\mathsf{Algorithm~{}\ref{alg:R1}(\mathbf{x},\mathbf{p})}$

14: else

15: if $\exists j_{c}\in\mathbf{x}_{c}\text{ s.t. }d_{b}(j_{c})=1$ then

16: $k_{b}\leftarrow j\in\mathbf{x}_{b}\cap K$

17: $\mathsf{Swap}(j_{c},k_{b})$

18: else

19: $k_{c}\leftarrow j\in\mathbf{x}_{c}\cap K$

20: $j_{a}\leftarrow\\{j\in\mathbf{x}_{a}\mid d_{b}(j)=1\\}$

21: $\mathsf{Transfer}(k_{c},a)$

22: $\mathsf{Transfer}(j_{a},b)$

23:return $(\mathbf{x},\mathbf{p})$

We show that:

######  Lemma 27.

Algorithm 9 returns an EFX competitive equilibrium.

###### Proof.

We have shown in the proof of Lemma 21 that when $R=2$ and $|N_{1}|=2$, if any
chore transfers occur in the execution of Algorithm 4 the output of Algorithm
4 is already EFX. Thus, before Line 5 of Algorithm 9, we must have a cost-
minimizing CE where each agent has only 1-chores or $K$-chores. We define
agent $a$ to be the agent with weakly greater 1-chores in $N_{1}$, $b$ to be
the remaining agent in $N_{1}$, and $c$ to be the agent in $N_{2}$. Without
loss of generality, we may assume that $b$ receives a $k$-chore before $a$ in
Algorithm 4 so that $|\mathbf{x}_{b_{k}}|\geq|\mathbf{x}_{a_{k}}|$.
Furthermore, if EFX-envy exists in $\mathbf{x}$, it must be that
$|\mathbf{x}_{b_{k}}|=|\mathbf{x}_{a_{k}}|+1$ and $b$ is the EFX-envious
agent.

In Line 5, we check that $c$ has strictly more total chores than $a$. If so,
we check in Line 6 that $c$ has strictly more $K$-chores than $b$. If both
conditions are satisfied, we directly transfer a $K$-chore from $c$ to $a$.
Since $|\mathbf{x}_{c_{k}}|>|\mathbf{x}_{b_{k}}|\geq|\mathbf{x}_{a_{k}}|+1$,
it remains after the transfer that $c$ has at least as many $K$-chores as $a$
and $b$, so they cannot EFX-envy $c$. It is also the case that $a$ and $b$ now
have equal $K$-chores, so they cannot EFX-envy each other. It remains that $c$
cannot EFX-envy $a$ or $b$, so the allocation is EFX. Else, if
$|\mathbf{x}_{c}|>|\mathbf{x}_{a}|$ but
$|\mathbf{x}_{c_{k}}|\leq|\mathbf{x}_{b_{k}}|$, we call Algorithm 5. Algorithm
5 either finds an EFX allocation immediately or calls Algorithm 6 or Algorithm
7. We have seen from Lemma 24 and Lemma 25 that these algorithms return an EFX
competitive equilibrium. The same argument applies to the execution of Line
13.

We then check the conditions in Line 14. It must be that
$|\mathbf{x}_{c}|\leq|\mathbf{x}_{a}|$ and
$|\mathbf{x}_{c_{1}}|<|\mathbf{x}_{b_{1}}|$. Note that $b$ is given a
$K$-chore. By the tiebreaker rules, it be then that
$|\mathbf{x}_{c}|\geq|\mathbf{x}_{b}|$. In addition, Invariant 4 guarantees
that $a$ and $b$ have fully balanced chores. Since
$|\mathbf{x}_{b_{k}}|=|\mathbf{x}_{a_{k}}|+1$, we have that
$|\mathbf{x}_{b}|\geq|\mathbf{x}_{a}|$. Combining these results we see that
$|\mathbf{x}_{c}|\geq|\mathbf{x}_{b}|\geq|\mathbf{x}_{a}|\geq|\mathbf{x}_{c}|$,
so in fact $|\mathbf{x}_{c}|=|\mathbf{x}_{b}|=|\mathbf{x}_{a}|$. It follows
then that $|\mathbf{x}_{a_{1}}|=|\mathbf{x}_{b_{1}}|+1>|\mathbf{x}_{c_{1}}+1$,
i.e., $a$ has exactly one more 1-chore than $b$, and $b$ has at least one
1-chore than $c$. Similarly, we have
$|\mathbf{x}_{c_{k}}|>|\mathbf{x}_{b_{k}}|=|\mathbf{x}_{a_{k}}|+1$, i.e., $c$
has at least one more $K$-chore than $b$, and $b$ has exactly one more
$K$-chore than $a$. In this case, we first check if $c$ has a chore $j_{c}$
which is disutility 1 for $b$. If so, we swap a $K$-chore of $b$ with $j_{c}$.
Then $a$ and $b$ cannot EFX-envy each other as they have equal $K$-chores now,
and they clearly continue to not EFX-envy $c$. At the same time $c$ can still
not EFX-envy $a$ as $c$ has disutility $k$ for all of her chores. Regarding
$b$, $c$ now views all but one of $b$’s chores as $k$. Then, since every agent
has exactly $q$ chores, we have
$\max_{j\in\mathbf{x}_{c}}d_{c}(\mathbf{x}_{c}\setminus
j)\leq(q-1)k\leq(q-1)k+1=d_{c}(\mathbf{x}_{b})$. So $c$ does not EFX-envy $b$
and the allocation is EFX.

Else, in Line 18 we still have that
$|\mathbf{x}_{c}|=|\mathbf{x}_{b}|=|\mathbf{x}_{a}|=q$ but now $b$ must have
disutility $k$ for all of $c$’s 1-chores. In this case, we transfer a
$K$-chore from $c$ to $a$, and a 1-chore from $a$ to $b$. Notably, since $c$
cannot be EFX-envied by $b$ here, it must be that $b$ EFX-envies $a$ all of
$a$’s chores are MPB for $b$. After this transfer $c$ continues to not envy
anybody. Now $a$ and $b$ have balanced chores and equal $K$-chores, so they do
not EFX-envy each other. It remains that $a$ and $c$ have balanced chores, and
$c$ still has at least as many $K$-chores as $a$, so $a$ does not EFX-envy
$c$. We now examine $b$ and $c$. It remains that $c$ has at least as many
$K$-chores as $b$, but $b$ now has two more total chores than $c$. Thus, after
removing a chore it can EFX-envy $c$ by at most 1. However, since we had to
have had EFX-envy between $a$ and $b$ when
$|\mathbf{x}_{b}|=|\mathbf{x}_{a}|$, by Lemma 20 it must be that $k>2$ in this
instance. Then, since $b$ has disutility $k>2$ for all of $c$’s 1-chores, $b$
does not EFX-envy $c$. Note, it must be that $c$ has a 1-chore since it could
not form a group by itself with no initial chores. Thus, the allocation is
EFX.

As only MPB chores are transferred, it must be that we maintain a CE and the
allocation returned is fPO. ∎

###  A.9 Algorithm 5: Obtaining EFX+fPO

######  Lemma 28.

Given a bivalued chore allocation instance $(N,M,D)$ with three agents,
Algorithm 5 computes an EFX + fPO allocation in strongly polynomial-time.

###### Proof.

We show that Algorithm 5 computes such an allocation in strongly-polynomial
time. Algorithm 5 first calls Algorithm 4. Then,

  * •

If $R=1$, Algorithm 5 is invoked. Then, it must be that either Algorithm 5
finds an EFX+fPO allocation, or one of Algorithms 6 and 7 is called. From
Lemmas 24 and 25 we know that an EFX+fPO allocation is found in either case.

  * •

If $R=2$, then either Algorithm 8 or Algorithm 9 is invoked. From Lemmas 26
and 27 we know that both algorithms give an EFX+fPO allocation.

  * •

Finally, if $R=3$, we show that the output of Algorithm 4 is already EFX+fPO.
Suppose the groups are such that $N_{1}=\\{a\\}$, $N_{2}=\\{a\\}$, and
$N_{3}=\\{c\\}$. From Lemma 8 we have that the allocation is balanced. Then
from Lemma 18 it cannot be that $c$ EFX-envies $a$ or $b$, and it cannot be
that $b$ EFX-envies $a$. If an agent belongs to a raised group, then by Lemma
14 they have only 1-chores, and so by Lemma 17 agents in such groups cannot
EFX-envy anyone. We then need only consider an agent in an unraised group who
may EFX-envy an agent in a group below her. By construction of the groups, we
have after $\mathsf{MakeInitGroups}$ is called that
$|\mathbf{x}_{a_{1}}|\geq|\mathbf{x}_{b_{1}}|\geq|\mathbf{x}_{c_{1}}|$. By
Condition (3) of Lemma 10, agents in unraised groups cannot have lost a chore.
By Lemma 11 they may only gain $k$-chores. Then, by the tiebreaker rules, we
have that
$|\mathbf{x}_{c_{k}}|\geq|\mathbf{x}_{b_{k}}|\geq|\mathbf{x}_{a_{k}}|$. Since
any unraised agent sees every $k$-chore of every other unraised agent as $k$
by Lemma 15 and the allocation is balanced, it cannot be that $a$ EFX-envies
$b$ or $c$, or that $b$ EFX-envies $c$. Thus, the allocation is EFX, and we
have seen from Lemma 8 that it is also fPO.

Thus, Algorithm 5 finds an EFX+fPO allocation. We now show that Algorithm 5
terminates in strongly-polynomial time. By Lemma 8 Algorithm 4 terminates in
$\mathsf{poly}(n,m)$ time. If Algorithm 5 is not called, it is clear that
Algorithm 8 and Algorithm 9 both terminate in $\mathsf{poly}(n,m)$ time as
they perform only constantly many chore transfers. If Algorithm 5 is called,
it performs only constantly many chore transfers before calling either
Algorithm 6 or Algorithm 7. Algorithm 6 performs only constantly many
transfers outside of its two while loops in Line 12 and Line 36. In both of
these cases, we show in Lemma 24 that the loop terminates after at most $m$
transfers. The same is done for Algorithm 7, which performs constantly many
tranfers outside of its two while loops in Lines 19 and 17. Lemma 25 shows
that they too terminate in at most $m$ transfers. Thus, it must be that
Algorithm 5 terminates in strongly-polynomial time. ∎

##  Appendix B Illustrative Examples

######  Example 6.

The algorithm of [26] does not necessarily return a PO allocation.

Zhou and Wu [26] use the $\mathsf{RoundRobin}$ procedure as a subroutine in
their algorithm. We show that $\mathsf{RoundRobin}$ need not give a PO
allocation. Consider the following instance with three agents $\\{a,b,c\\}$
and three chores $\\{j_{1},j_{2},j_{3}\\}$.

| $j_{1}$ | $j_{2}$ | $j_{3}$
---|---|---|---
$a$ | 1 | 1 | 1
$b$ | $k$ | 1 | $k$
$c$ | 1 | $k$ | $k$

Consider an execution of the $\mathsf{RoundRobin}$ algorithm using the order
$a,b,c$. First $a$ picks $j_{1}$, then $b$ picks $j_{2}$, then $c$ picks
$j_{3}$. This gives a disutility vector of $v_{1}=(1,1,k)$ for agents $a,b,c$.
However allocating $j_{3}$ to $a$, $j_{2}$ to $b$ and $j_{1}$ to $c$ gives a
utility vector of $v_{2}=(1,1,1)$, which Pareto-dominates $v_{1}$. Hence the
outcome of $\mathsf{RoundRobin}$ is not necessarily PO, thus neither is
outcome of the algorithm of [26].

######  Example 7.

Execution of Algorithm $\mathsf{MakeInitGroups}$

Consider the chore allocation instance below (with only $L$-chores).

| $j_{1}$ | $j_{2}$ | $j_{3}$ | $j_{4}$ | $j_{5}$
---|---|---|---|---|---
$a$ | 1 | 1 | 1 | 1 | $k$
$b$ | 1 | 1 | 1 | $k$ | 1
$c$ | $k$ | $k$ | $k$ | $k$ | 1

$\mathsf{MakeInitGroups}$ begins with some cost-minimizing allocation
$\mathbf{x}$, say:

  * •

$\mathbf{x}_{a}=\\{j_{1},j_{2},j_{3},j_{4}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{5}\\}$,

  * •

$\mathbf{x}_{c}=\\{\\}$.

Note that all chore payments are 1. Then, $a$ is the big earner, and she
pEF1-envies an agent who is in her component, such as $b$. Indeed,
$(b,j_{1},a)$ is a shortest path from $b$ to $a$ in the MPB graph. We thus
transfer $j_{1}$ from $a$ to $b$, giving us:

  * •

$\mathbf{x}_{a}=\\{j_{2},j_{3},j_{4}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{5},j_{1}\\}$,

  * •

$\mathbf{x}_{c}=\\{\\}$.

It remains that $a$ is the big earner but she no longer pEF1-envies $b$, only
$c$. As $c$ is in her component via shortest path $(c,j_{5},b,j_{2},a)$, we
transfer $j_{5}$ to $c$, so we have:

  * •

$\mathbf{x}_{a}=\\{j_{2},j_{3},j_{4}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{1}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{5}\\}$.

Still $a$ is the big earner. While she pEF1-envies both $b$ and $c$, we note
that $c$ is no longer in her component. We do however have shortest path
($b,j_{2},a$). We transfer $j_{2}$ to $b$ to get:

  * •

$\mathbf{x}_{a}=\\{j_{3},j_{4}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{1},j_{2}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{5}\\}$.

Now we may take either $a$ or $b$ to be the big earner. Suppose we choose $b$.
Only $a$ is in $b$’s component, and $b$ does not pEF1-envy her. We are thus
ready to form agent group $N_{1}=\\{b,a\\}$. Then, $c$ is left as the only
agent not in a group, so she cannot pEF1-envy any other agent without a group.
So, $c$ forms her own group $N_{2}=\\{c\\}$. As all agents now belong to a
group, $\mathsf{MakeInitGroups}$ outputs the following:

  * •

$\mathbf{x}_{a}=\\{j_{3},j_{4}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{1},j_{2}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{5}\\}$.

  * •

$N_{1}=\\{b,a\\}$

  * •

$N_{2}=\\{c\\}$

It is easy to verify that this output satisfies properties (i), (ii), and
(iii) of Section 5.2.

######  Example 8.

Execution of Algorithm 4, Algorithm 5, Algorithm 5, and Algorithm 7

Consider the chore allocation instance below where $k=5$.

| $j_{1}$ | $j_{2}$ | $j_{3}$ | $j_{4}$ | $j_{5}$ | $j_{6}$ | $j_{7}$ | $j_{8}\ldots j_{11}$
---|---|---|---|---|---|---|---|---
$a$ | 1 | 1 | 1 | 1 | 1 | 1 | $k$ | $k$
$b$ | 1 | 1 | 1 | 1 | 1 | 1 | 1 | $k$
$c$ | 1 | 1 | 1 | 1 | $k$ | 1 | 1 | $k$

Suppose Algorithm 4 outputs one agent group (so no transfers are required)
with the following allocation:

  * •

$\mathbf{x}_{a}=\\{j_{1},j_{2},j_{3},j_{8}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{4},j_{5},j_{9},j_{10}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{6},j_{7},j_{11}\\}$.

Then, Algorithm 5 calls Algorithm 5. Since $b$ EFX-envies both $a$ and $c$, we
swap $j_{9}$ and $j_{3}$, giving us the following allocation where the only
EFX-envy is from $a$ towards $b$:

  * •

$\mathbf{x}_{a}=\\{j_{1},j_{2},j_{9},j_{8}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{4},j_{5},j_{3},j_{10}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{6},j_{7},j_{11}\\}$.

Finally Algorithm 7 is called. We first execute on Line 14, where we swap
$j_{4}$ and $j_{7}$, giving us:

  * •

$\mathbf{x}_{a}=\\{j_{1},j_{2},j_{9},j_{8}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{7},j_{5},j_{3},j_{10}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{6},j_{4},j_{11}\\}$.

We see that $a$ still EFX-envies $c$, so we execute Line 17, swapping $j_{9}$
with $j_{4}$. We now have:

  * •

$\mathbf{x}_{a}=\\{j_{1},j_{2},j_{4},j_{8}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{7},j_{5},j_{3},j_{10}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{6},j_{9},j_{11}\\}$.

Still we do not have EFX, as $c$ now EFX-envies $a$. Then, we execute on Line
23, transferring $j_{6}$ to $b$. This gives us our final EFX allocation:

  * •

$\mathbf{x}_{a}=\\{j_{1},j_{2},j_{4},j_{8}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{6},j_{7},j_{5},j_{3},j_{10}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{9},j_{11}\\}$.

######  Example 9.

Execution of Algorithm 4, Algorithm 5, Algorithm 5, and Algorithm 6

Consider the chore allocation instance below where $k=5$.

| $j_{1}$ | $j_{2}$ | $j_{3}$ | $j_{4}$ | $j_{5}$ | $j_{6}\ldots j_{10}$
---|---|---|---|---|---|---
$a$ | 1 | 1 | 1 | $k$ | $k$ | $k$
$b$ | 1 | 1 | 1 | 1 | 1 | $k$
$c$ | 1 | 1 | 1 | 1 | 1 | $k$

Suppose Algorithm 4 outputs one agent group (so no transfers are required)
with the following allocation:

  * •

$\mathbf{x}_{a}=\\{j_{1},j_{2},j_{6}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{3},j_{7},j_{8}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{4},j_{5},j_{9},j_{10}\\}$.

Algorithm 5 again calls Algorithm 5. Since both $b$ and $c$ EFX-envy $a$, and
$|\mathbf{x}_{c}|>|\mathbf{x}_{a}|$, we transfer a $K$-chore from $c$ to $a$,
giving us:

  * •

$\mathbf{x}_{a}=\\{j_{1},j_{2},j_{6},j_{10}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{3},j_{7},j_{8}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{4},j_{5},j_{9}\\}$.

Now only $b$ EFX-envies $c$ and we move to Algorithm 6. Here we execute at
Line 37 and first transfer a 1-chore from $a$ to $c$, then transfer a 1-chore
from $b$ to $c$. This gives us the EFX allocation:

  * •

$\mathbf{x}_{a}=\\{j_{2},j_{6},j_{10}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{7},j_{8}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{4},j_{5},j_{9},j_{1},j_{3}\\}$.

######  Example 10.

Execution of Algorithm 4, Algorithm 5, and Algorithm 8

Consider the chore allocation instance below where $k=5$.

| $j_{1}$ | $j_{2}$ | $j_{3}$ | $j_{4}$ | $j_{5}$ | $j_{6}$ | $j_{7}$ | $j_{8}$ | $j_{9}\ldots j_{11}$
---|---|---|---|---|---|---|---|---|---
$a$ | 1 | 1 | 1 | 1 | 1 | 1 | $k$ | 1 | $k$
$b$ | $k$ | $k$ | $k$ | $k$ | $k$ | $k$ | 1 | $k$ | $k$
$c$ | $k$ | $k$ | $k$ | $k$ | $k$ | $k$ | 1 | 1 | $k$

Here Algorithm 4 creates two agent groups $\\{a\\}$ and $\\{b,c\\}$ with the
following initial allocation:

  * •

$\mathbf{x}_{a}=\\{j_{1},j_{2},j_{3},j_{4},j_{5},j_{6}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{7},j_{9}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{8},j_{10},j_{11}\\}$.

Then, to balance the chores, we raise the payments of chores belonging to $a$
and transfer a chore to both $b$ and $c$, giving us:

  * •

$\mathbf{x}_{a}=\\{j_{1},j_{2},j_{3},j_{4}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{7},j_{9},j_{5}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{8},j_{10},j_{11},j_{6}\\}$.

Now the chores are balanced and Algorithm 8 is called. Since $c$ still EFX-
envies $b$, we next execute Line 14, raising payments of $\mathbf{x}_{a}$ and
then transferring a chore, say $j_{4}$, to $b$. This gives us the EFX
allocation:

  * •

$\mathbf{x}_{a}=\\{j_{1},j_{2},j_{3}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{7},j_{9},j_{5},j_{4}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{8},j_{10},j_{11},j_{6}\\}$.

######  Example 11.

Execution of Algorithm 4, Algorithm 5, and Algorithm 9

Consider the chore allocation instance below where $k=5$.

| $j_{1}$ | $j_{2}$ | $j_{3}$ | $j_{4}$ | $j_{5}$ | $j_{6}$ | $j_{7}\ldots j_{12}$
---|---|---|---|---|---|---|---
$a$ | 1 | 1 | 1 | 1 | 1 | $k$ | $k$
$b$ | 1 | 1 | 1 | 1 | 1 | $k$ | $k$
$c$ | $k$ | $k$ | $k$ | $k$ | $k$ | 1 | $k$

Here Algorithm 4 creates two agent groups $\\{a,b\\}$ and $\\{c\\}$ with the
following initial allocation:

  * •

$\mathbf{x}_{a}=\\{j_{1},j_{2},j_{3},j_{7}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{4},j_{5},j_{8},j_{9}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{6},j_{10},j_{11},j_{12}\\}$.

Since the chores are balanced, no transfers are needed. Algorithm 5 then calls
Algorithm 9. We satisfy the conditions in Lines 11, 14 and 18 and execute
accordingly. We transfer a $K$-chore from $c$ to $a$ and a 1-chore from $a$ to
$b$. This gives us the final EFX allocation:

  * •

$\mathbf{x}_{a}=\\{j_{2},j_{3},j_{7},j_{12}\\}$,

  * •

$\mathbf{x}_{b}=\\{j_{4},j_{5},j_{8},j_{9},j_{1}\\}$,

  * •

$\mathbf{x}_{c}=\\{j_{6},j_{10},j_{11}\\}$.


# Prospects of detection of subsolar mass primordial black hole and white
dwarf binary mergers

Takahiro S. Yamamoto<EMAIL_ADDRESS>Ryoto Inui
<EMAIL_ADDRESS>Department of Physics, Nagoya University,
Furo-cho Chikusa-ku, Nagoya 464-8602, Japan  Yuichiro Tada
<EMAIL_ADDRESS>Institute for Advanced Research, Nagoya
University, Furo-cho Chikusa-ku, Nagoya 464-8601, Japan  Department of
Physics, Nagoya University, Furo-cho Chikusa-ku, Nagoya 464-8602, Japan
Shuichiro Yokoyama<EMAIL_ADDRESS>Kobayashi Maskawa Institute, Nagoya
University, Chikusa, Aichi 464-8602, Japan  Kavli IPMU (WPI), UTIAS, The
University of Tokyo, Kashiwa, Chiba 277-8583, Japan

###### Abstract

The subsolar mass primordial black hole (PBH) attracts attention as robust
evidence of its primordial origin against the astrophysical black hole. Not
only with themselves, PBHs can also form binaries with ordinary astrophysical
objects, catching them by gravitational wave (GW) bremsstrahlung. We discuss
the detectability of the inspiral GWs from binaries consisting of a PBH and a
white dwarf (WD) by using space-borne gravitational wave interferometers like
DECIGO. The conservative assessment shows the expected event number in three
years by DECIGO is $\mathcal{O}(10^{-6})$ for $M_{\textrm{PBH}}\sim
0.1M_{\odot}$. Possible enhancement mechanisms of WD-PBH binary formation may
amplify this event rate. We discuss how large enhancement associated with WDs
is required to detect WD-PBH merger events without violating the existing
constraints on the PBH-PBH merger by the ground-based detector.

##  I introduction

PBH are hypothetical black holes that could have formed in the early universe
before star formation [1, 2, 3, 4]. They have attracted attention more and
more as a candidate for dark matter [5], seeds of supermassive black holes [6,
7] (see, e.g., Refs. [8, 9] for implications of supermassive PBH by high-
redshift luminous galaxies found by James Webb Space Telescope), planetary-
mass microlensing objects towards the Galactic bulge detected by the Optical
Gravitational Lensing Experiment [10], hypothetical Planet 9 in our solar
system [11, 12], triggers of faint supernovae called Calcium-rich gap
transients [13], etc. (see, e.g., Refs. [14, 15] for recent reviews on
astrophysical implications of PBHs). A significant property of PBHs is that
they can have a wide range of mass even smaller than the solar mass, unlike
the astrophysical black holes that are formed through the gravitational
collapse of massive stars. Thus, the detection of subsolar mass black holes is
a smoking gun for the existence of PBH.

While subsolar mass PBHs have been searched mainly by the gravitational
microlensing111So far any microlensing search for PBHs has not detected a PBH,
but several candidate events are reported, represented by the single event by
the Subaru Hyper Suprime-Camera [16]. so far (see Ref. [17] for a summary of
the existing observational constraints on the PBH abundance), GWs have been
attracting much attention as a new tool to probe the subsolar mass PBH. In the
early universe, PBH can form binaries by the tidal torque exerted by a
neighboring PBH [18, 19, 20]. After forming the binaries, their orbits
gradually shrink as they emit the GW. Depending on the separation and the
component masses of binaries, their GWs can have the frequency to which the GW
detectors are sensitive. LIGO–Virgo–KAGRA (LVK) collaboration has already
obtained the constraints on the PBH abundance from the search of the subsolar
mass black hole binaries [21, 22, 23, 24]. Several groups independently report
candidates of GW from subsolar mass black hole binaries [25, 26, 27].
Furthermore, Pujolas et al. [28] estimated the detectability of the PBH
binaries in various future observations, such as Einstein Telescope [29], LISA
[30], and DECIGO [31, 32].

While the dominant formation channel of the PBH binaries is the three-body
interaction in the early universe, the GW bremsstrahlung can lead to PBH
binaries in the late universe [33, 34, 20]. In this channel, an unbounded PBH
is captured by another PBH through the close encounter due to the energy loss
by the GW emission. This channel can form not only PBH-PBH binaries but also
the binaries of a PBH and another astrophysical compact object. For example,
Tsai et al. [35] and Sasaki et al. [36] estimate the merger rate of binaries
consisting of a PBH and a neutron star (NS). Tsai et al. discussed the
scenario in which GW170817 is an NS-PBH merger event. Sasaki et al. shows that
the NS-PBH binaries are subdominant, and their merger rate is not enough to
explain LVK’s NS-BH merger candidates.

In this work, as another system, we consider the binaries consisting of a
subsolar mass PBH and a WD. The merger frequencies of the WD-PBH binaries are
$\mathcal{O}(10^{-2}\text{--}10^{-1})\,$\mathrm{H}\mathrm{z}$$, which are
interesting targets of DECIGO. However, numerical simulations of the WD-PBH
merger are yet to be done. We then focus on the inspiral GW and investigate
the detectability of the WD-PBH merger event via the maximum detectable
distance of inspiral events with a sufficiently large signal-to-noise ratio
(SNR) in the DECIGO observation. As WD-PBH inspiral GWs can be mimicked by WD-
WD ones, we also discuss whether they can be distinguished only by the
differences in the merger frequencies. Finally, we estimate the expected event
number of the WD-PBH binaries in the three-year DECIGO mission.

The paper is organized as follows. In Sec. II, we discuss the detectability
and the distinguishability of WD-PBH binaries by assuming three years of
observation of DECIGO. In Sec. III, we present the merger rate estimation, and
Sec. IV is devoted to the conclusions and future perspective. The details of
the merger rate calculation are summarized in Appendix. A. Throughout the
paper, we suppose that the PBH mass function is monochromatic.

##  II Detectability by DECIGO

###  II.1 Signal-to-noise ratio for the inspiraling WD-PBH binary

Figure 1: Noise PSD of DECIGO (10).

In this subsection, we evaluate the detectability of WD-PBH binaries only
through the inspiral GWs. The most optimistic scenario, which we hereafter
assume, is that DECIGO keeps observing the binary inspiral for three years
from the beginning of the observation run, and the binary merges just after
the end of the observation run. The detectability is quantified by the SNR of
the inspiral waveform from the initial frequency $f_{\mathrm{min}}$ to the
final frequency $f_{\mathrm{max}}$ in the detector frame. Note that we
basically suppose quasi-circular orbits for WD-PBH binaries for brevity though
Cholis et al. [37] suggested that the binaries formed in the GW bremsstrahlung
have high eccentricities and merge within a few years, typically keeping their
eccentricities even at the merger.

The measured signal in the frequency domain is obtained by

$\tilde{h}(f)=F_{+}\tilde{h}_{+}(f)+F_{\times}\tilde{h}_{\times}(f)\,,$ (1)

where $F_{+}$ and $F_{\times}$ are the antenna pattern functions, and
$\tilde{h}_{+/\times}$ are plus/cross modes expressed by

$\tilde{h}_{+}(f)=(1+\cos^{2}\iota)\tilde{h}_{0}(f)\,,\quad\tilde{h}_{\times}(F)=-2i\cos\iota\tilde{h}_{0}(f)\,,$
(2)

with the inclination angle $\iota$. $\tilde{h}_{0}$ is given by

$\tilde{h}_{0}(f)=\sqrt{\frac{5}{96}}\pi^{-2/3}\quantity(\frac{(1+z)G\mathcal{M}_{\textrm{c}}}{c^{3}})^{5/6}\frac{c}{d_{\mathrm{L}}}f^{-7/6}\mathrm{e}^{i\Psi(f)}\,,$
(3)

where $\Psi(f)$ is the GW phase, $G$ is the Newtonian constant of gravity, and
$c$ is the speed of light. $\mathcal{M}_{\textrm{c}}$ is the chirp mass
defined by

$\mathcal{M}_{\textrm{c}}\coloneqq\frac{(m_{1}m_{2})^{3/5}}{(m_{1}+m_{2})^{1/5}}\,,$
(4)

with the component masses $m_{1}$ and $m_{2}$ in the source frame. $z$ and
$d_{\mathrm{L}}$ are the source redshift and the luminosity distance,
respectively. Their relation is given by

$d_{\textrm{L}}=(1+z)\int^{z}_{0}\frac{c\differential{z}}{H(z)}\,.$ (5)

We employ astropy [38] to calculate the luminosity distance and use the
cosmological parameters obtained by Planck18 [39]. Given the detector noise
spectrum $S_{\mathrm{n}}(f)$, we define the SNR by

$\varrho_{0}=\quantity[4\int^{\infty}_{0}\frac{|\tilde{h}(f)|^{2}}{S_{\mathrm{n}}(f)}\differential{f}]^{1/2}\,.$
(6)

Following Ref. [40], we estimate the SNR averaged over the extrinsic
parameters (e.g., the source position in the sky) by the formula

$\varrho=\frac{8}{5}\frac{A}{d_{\textrm{L}}}\sqrt{I_{7}}\,,$ (7)

where $I_{7}$ and $A$ are given by

$\displaystyle
I_{7}=\int^{f_{\mathrm{max}}}_{f_{\mathrm{min}}}\differential{f}\
\frac{f^{-7/3}}{S_{\textrm{n}}(f)}\,,$ (8) $\displaystyle
A=\sqrt{\frac{5}{96}}\pi^{-2/3}c\left(\frac{(1+z)G\mathcal{M}_{\textrm{c}}}{c^{3}}\right)^{5/6}\,.$
(9)

The integration range $[f_{\mathrm{min}},f_{\mathrm{max}}]$ in Eq. (8) is
obtained below as a frequency evolution during the three-year observation of
DECIGO.

In the present work, we use the DECIGO’s design power spectral density (PSD)
[41]

$S^{\mathrm{DECIGO}}_{\textrm{n}}(f)=\bqty{7.05\times
10^{-48}\quantity(1+\left(\frac{f}{f_{\textrm{p}}}\right)^{2})+4.8\times
10^{-51}\quantity(\frac{f}{$1\text{\,}\mathrm{H}\mathrm{z}$})^{-4}\frac{1}{1+\quantity(\frac{f}{f_{\textrm{p}}})^{2}}+5.33\times
10^{-52}\quantity(\frac{f}{$1\text{\,}\mathrm{H}\mathrm{z}$})^{-4}}\,$\mathrm{H}\mathrm{z}^{-1}$,$
(10)

with $f_{\textrm{p}}=$7.36\text{\,}\mathrm{H}\mathrm{z}$$. Figure 1 shows the
noise PSD of DECIGO.

Here, we focus only on GW from the inspiral phase. We assume that
$f_{\mathrm{max}}$ is given by

$\displaystyle f_{\textrm{max}}=\frac{f_{\textrm{merge}}^{\textrm{WD-
PBH}}}{1+z}\,,$ (11)

where $f_{\mathrm{merge}}^{\text{WD-PBH}}$ is a typical merger frequency in
the source frame. In the case of WD-PBH binaries, WD would be tidally
disrupted in the final stage of the inspiral phase, and then
$f_{\mathrm{merge}}^{\text{WD-PBH}}$ is determined by the binary separation
equal to the tidal radius, $R_{\textrm{tidal}}$, which is estimated by
balancing the tidal force and the WD’s self-gravity as

$\displaystyle F_{\textrm{tidal}}=F_{\textrm{grav}}\quad\text{at
$r=R_{\mathrm{tidal}}$}$ $\displaystyle\Rightarrow\quad$
$\displaystyle\frac{2GM_{\textrm{PBH}}}{r^{2}}\left(\frac{R_{\mathrm{WD}}}{r}\right)=\frac{GM_{\textrm{WD}}}{R_{\mathrm{WD}}^{2}}\quad\text{at
$r=R_{\mathrm{tidal}}$}$ $\displaystyle\Rightarrow\quad$ $\displaystyle
R_{\mathrm{tidal}}=\left(\frac{2M_{\textrm{PBH}}}{M_{\textrm{WD}}}\right)^{1/3}R_{\mathrm{WD}}\,,$
(12)

where $M_{\textrm{WD}}$ and $M_{\textrm{PBH}}$ are respectively masses of the
WD and the PBH, and $R_{\mathrm{WD}}$ is the radius of the WD.We use the
following mass-radius relation:

$R_{\textrm{WD}}=0.0126R_{\odot}\left(\frac{M_{\textrm{WD}}}{M_{\odot}}\right)^{-1/3}\,,$
(13)

which is proposed in Ref. [42]. The corresponding frequency is given by

$f_{\textrm{tidal}}=\frac{1}{\pi}\sqrt{\frac{G(M_{\textrm{WD}}+M_{\textrm{PBH}})}{R_{\textrm{tidal}}^{3}}}\,.$
(14)

Equation (II.1) implies that $R_{\textrm{tidal}}$ is smaller than the radius
of the WD, $R_{\mathrm{WD}}$, if $2M_{\textrm{PBH}}<M_{\textrm{WD}}$. In such
cases, as the merger frequency, $f_{\mathrm{merge}}^{\text{WD-PBH}}$, we
should use the Kepler frequency at the time when the orbital radius becomes
the sum of the radii of the WD and the PBH ($R_{\mathrm{PBH}}$):

$f_{\textrm{Kepler}}=\frac{1}{\pi}\sqrt{\frac{G(M_{\textrm{WD}}+M_{\textrm{PBH}})}{(R_{\mathrm{WD}}+R_{\mathrm{PBH}})^{3}}}.$
(15)

The merger frequency is determined by the one realized earlier:

$f_{\mathrm{merge}}^{\text{WD-
PBH}}=\min[f_{\textrm{Kepler}},f_{\textrm{tidal}}]\,.$ (16)

As for $f_{\mathrm{min}}$, we introduce $f_{\textrm{start}}$ as the
corresponding frequency in the source frame:

$f_{\textrm{min}}=\frac{f_{\textrm{start}}}{1+z}\,.$ (17)

The duration that the binary requires to evolve from the frequency $f_{1}$ to
$f_{2}$ ($>f_{1}$) is given by (see Refs. [43, 44, 45])

$t(f_{2})-t(f_{1})=\frac{5}{256}\left(\frac{G\mathcal{M}_{\textrm{c}}}{c^{3}}\right)^{-5/3}\pi^{-8/3}(f_{1}^{-8/3}-f_{2}^{-8/3})\,.$
(18)

From this expression, by assuming that the binary merges at the end of the
observation period in our optimistic scenario, we can evaluate the frequency
at the onset of the observation run (i.e., three years ($=T_{\mathrm{obs}}$)
before the merger), which is corresponding to $f_{\textrm{start}}$, as

$T_{\textrm{obs}}=t(f_{\mathrm{merge}}^{\text{WD-
PBH}})-t(f_{\textrm{start}})\,,$ (19)

and it gives

$\displaystyle f_{\textrm{start}}$
$\displaystyle=\left\\{\quantity(f_{\textrm{merge}}^{\textrm{WD-
PBH}})^{-8/3}+T_{\textrm{obs}}\frac{256\pi^{8/3}}{5}\left(\frac{G\mathcal{M}_{\textrm{c}}}{c^{3}}\right)^{5/3}\right\\}^{-3/8}.$
(20)

Figure 2: Horizon distance as a function of $M_{\textrm{PBH}}$ and
$M_{\textrm{WD}}$. The black solid line indicates
$1\text{\,}\mathrm{G}\mathrm{p}\mathrm{c}$. There are breaks around
$M_{\textrm{PBH}}\sim(2\text{--}5)\times 10^{-1}M_{\odot}$. Around there, the
merger frequency, defined by Eq. (16), changes drastically because the tidal
frequency $f_{\mathrm{tidal}}$ and the Kepler frequency $f_{\mathrm{Kepler}}$
have different mass dependencies.

We set the detection threshold on the SNR at $8$ [44] and refer to the maximum
distance of the detectable binaries as the horizon distance. Figure 2 shows
the horizon distance depending on the masses of a PBH and a WD. From this
figure, we can find that we would detect the GW from a WD-PBH binary at
$\sim$1\text{\,}\mathrm{G}\mathrm{p}\mathrm{c}$$ for $M_{\textrm{WD}}\gtrsim
0.8M_{\odot}$ and $M_{\textrm{PBH}}\gtrsim 0.5M_{\odot}$.

|  |
---|---|---
Figure 3: Comparison of the relevant frequencies for WD-PBH and WD-WD
binaries. Black solid and dotted lines are maximum and minimum frequencies
realized by WD-PBH inspirals, respectively. Orange dashed lines show the
merger frequency of the WD-WD binary whose component masses are
$M_{\textrm{WD,1}}$ and $M_{\textrm{WD,2}}$. If the orange dashed line is
below the black dotted line, WD-PBH and WD-WD binaries can be distinguished
only by the inspiral waveform. We assume three years of observation.

###  II.2 Distinguishability between WD-PBH and WD-WD

In the previous subsection, we evaluated the detectability of WD-PBH binaries.
The GW waveform from a compact binary is characterized by the chirp mass
$\mathcal{M}_{\textrm{c}}$. Therefore, if a WD-WD binary and a WD-PBH binary
have similar component masses, they can not be distinguished only from the
chirp mass. However, the merger frequency can be different because of the
difference in the radii of WD and PBH. We discuss, here, how we can
distinguish WD-WD binaries and WD-PBH binaries by the difference in the merger
frequencies.

The merger frequency of WD-WD binaries, which is denoted by
$f_{\textrm{merge}}^{\textrm{BWD}}$, is estimated as well as the merger
frequency of WD-PBH binaries. We consider the WD-WD binary with the masses of
$M_{\textrm{WD,1}}$ and $M_{\textrm{WD,2}}$. We use $R_{\textrm{tidal,1}}$ to
denote the tidal radius in which the first WD ($M_{\textrm{WD,1}}$ and
$R_{\textrm{WD,1}}$) is disrupted by the companion WD ($M_{\textrm{WD,2}}$).
It is estimated by

$R_{\textrm{tidal,1}}=\left(\frac{2M_{\textrm{WD,2}}}{M_{\textrm{WD,1}}}\right)^{1/3}R_{\textrm{WD,1}}\,.$
(21)

By changing the labels $1\leftrightarrow 2$, the tidal radius
$R_{\textrm{tidal,2}}$ can be defined. When the two WDs approach as close as
the radius

$R_{\textrm{tidal}}^{\textrm{BWD}}=\max[R_{\textrm{tidal,1}},R_{\textrm{tidal,2}}]\,,$
(22)

the lighter one is tidally disrupted. We define the frequency corresponding to
$R_{\textrm{tidal}}^{\textrm{BWD}}$ by

$f_{\textrm{tidal}}^{\textrm{BWD}}=\frac{1}{\pi}\sqrt{\frac{G(M_{\mathrm{WD,1}}+M_{\mathrm{WD,2}})}{\quantity(R_{\textrm{tidal}}^{\textrm{BWD}})^{3}}}\,.$
(23)

When $M_{\textrm{WD,1}}\simeq M_{\textrm{WD,2}}$,
$R_{\textrm{tidal}}^{\textrm{BWD}}$ could be smaller than the sum of the WD’s
radii. In such cases, as the merger frequency, we should employ the Kepler
frequency at the time when the binary separation is equal to the sum of the
WD’s radii:

$f_{\textrm{Kepler}}^{\textrm{BWD}}=\frac{1}{\pi}\sqrt{\frac{G(M_{\textrm{WD,1}}+M_{\textrm{WD,2}})}{(R_{\textrm{WD,1}}+R_{\textrm{WD,2}})^{3}}}\,.$
(24)

Thus, as in the case of the WD-PBH binary, the merger frequency of the WD-WD
binary is given by the lower one:

$f_{\textrm{merge}}^{\textrm{BWD}}=\min[f_{\textrm{Kepler}}^{\textrm{BWD}},f_{\textrm{tidal}}^{\textrm{BWD}}]\,.$
(25)

Figure 3 shows the comparisons between the merger frequencies, $f^{\textrm{WD-
PBH}}_{\textrm{merge}}$ and $f^{\textrm{BWD}}_{\textrm{merge}}$, and the
initial frequency $f^{\textrm{WD-PBH}}_{\textrm{start}}$ of the WD-PBH binary
as functions of $M_{\textrm{PBH}}$. If the frequency satisfies
$f_{\textrm{merge}}^{\textrm{BWD}}\leq f_{\textrm{start}}^{\textrm{WD-PBH}}$,
the WD-PBH binary merges at the frequency that WD-WD binaries with the same
masses never reach. So, we can distinguish the WD-PBH binary and the WD-WD
binary by the difference in the merger frequencies.

Figure 4 zooms in on the right panel of Fig. 3, focusing on the regions where
$0.5~{}M_{\odot}\lesssim M_{\mathrm{WD},1}\lesssim 0.6~{}M_{\odot}$. For some
mass ranges,the merger frequency of the WD-WD binary is located between
$f^{\textrm{WD-PBH}}_{\textrm{start}}$ and $f^{\textrm{WD-
PBH}}_{\textrm{merge}}$. This means that the waveforms of the WD-PBH binary
and the WD-WD binary overlap in the frequency domain. If this overlap is
significant, it might be difficult to distinguish between the WD-PBH and WD-WD
binaries. We define the quantity named “distinguishability” by

$x_{\mathrm{dis}}\coloneqq\frac{\varrho^{2}-\varrho^{\prime
2}}{\varrho^{2}}\,,$ (26)

where $\varrho^{\prime}$ is defined by

$\varrho^{\prime}=\frac{8}{5}\frac{A}{d_{\textrm{L}}}\sqrt{I^{\prime}_{7}}\,,$
(27)

with “wrongly” replacing the WD-PBH merger frequency with that of the WD-WD
one:

$I^{\prime}_{7}=\begin{cases}\displaystyle\int^{f_{\textrm{merge}}^{\textrm{BWD}}/(1+z)}_{f_{\textrm{start}}^{\textrm{WD-
PBH}}/(1+z)}\differential{f}\frac{f^{-7/3}}{S_{\textrm{n}}(f)}&\text{for
$f_{\textrm{start}}^{\textrm{WD-PBH}}<f_{\textrm{merge}}^{\textrm{BWD}}$,}\\\
0&\text{for $f_{\textrm{start}}^{\textrm{WD-PBH}}\geq
f_{\textrm{merge}}^{\textrm{BWD}}$.}\end{cases}$ (28)

If $f^{\textrm{BWD}}_{\textrm{merge}}$ is equal or very close to
$f^{\textrm{WD-PBH}}_{\textrm{merge}}$, we could not judge whether the binary
is a WD-PBH binary or a WD-WD binary only from their inspiral parts because
they have the same frequency evolution during the observational period. In
such case, $\varrho^{\prime}$ equals $\varrho$. In this sense, the
“distinguishability” $x_{\mathrm{dis}}$ is a quantity representing how likely
we can distinguish a WD-PBH binary from a WD-WD binary. That is, for
$x_{\mathrm{dis}}\simeq 1$, we can distinguish a WD-PBH binary and a WD-WD
binary, but for $x_{\mathrm{dis}}\simeq 0$, we can not distinguish a WD-PBH
binary and a WD-WD binary. Figure 5 shows the distinguishability
$x_{\mathrm{dis}}$ as a function of the WD’s mass with
$M_{\textrm{PBH}}=1.2M_{\odot}$ and $1.4M_{\odot}$. Let us focus on the case
of $M_{\textrm{PBH}}=1.4M_{\odot}$ (thick black line in Fig. 5). For
$M_{\textrm{WD}}\gtrsim 0.55M_{\odot}$, the distinguishability
$x_{\mathrm{dis}}$ of a WD-PBH binary and a WD-WD binary is almost unity. This
is because the merger frequency of a WD-WD binary is lower than the initial
frequency of a WD-PBH binary, as Fig. 4 shows. On the other hand, for
$M_{\textrm{WD}}\lesssim 0.55M_{\odot}$, the merger frequencies of a WD-WD
binary and a WD-PBH binary are the same. This leads to no distinguishability
in SNR. Even in this case, WD-PBH binaries and WD-WD binaries could be
distinguished once we go beyond the conservative approach. We can consider
various effects related to the coalescence (e.g., tidal deformability). We
will discuss it in more detail in the concluding section.

Figure 4: Enlarged figure of the right panel of Fig. 3. While binaries are
distinguishable for $M_{\mathrm{WD},1}\gtrsim 0.56M_{\odot}$, it is hard for
$M_{\mathrm{WD},1}\lesssim 0.56M_{\odot}$ and impossible for
$M_{\mathrm{WD},1}\lesssim 0.54M_{\odot}$ to distinguish them only by the
inspiral waveform and hence the merger waveform is required which is beyond
the scope of the paper. Figure 5: Distinguishability (26) in black lines. The
thick (thin) line shows the PBH mass of 1.4 $M_{\odot}$ (1.2 $M_{\odot}$). The
distance to the source is chosen such that the reference SNR $\varrho=30$.
Although the integrals $I_{7}$ and $I^{\prime}_{7}$ depend on the redshift,
the result is not much affected by the value of the reference SNR. The blue
dashed lines represent the tidal radius of the WD-PBH binary. The orange dash-
dotted lines represent the contact radius of the WD-WD binary whose component
masses are $M_{\textrm{WD}}$ and $1.4M_{\odot}$ ($1.2M_{\odot}$).

##  III Expected number of events

In this section, we estimate the merger rate of the WD-PBH binaries. The
expected number of events is given by

$N_{\textrm{event}}^{\text{WD-
PBH}}=T_{\textrm{obs}}\int^{M_{\textrm{WD,max}}}_{M_{\textrm{WD,min}}}\frac{\differential{M_{\textrm{WD}}}}{M_{\textrm{WD}}}\
\int^{z_{\textrm{up}}}_{0}\differential{z}\derivative{V}{z}\int^{M_{\textrm{h,max}}}_{M_{\textrm{h,min}}}\differential{M_{\textrm{h}}}\derivative{n_{\textrm{halo}}}{M_{\textrm{h}}}\derivative{R_{\textrm{WD-
PBH}}}{\ln M_{\textrm{WD}}}\,,$ (29)

with the comoving volume $V(z)$ within the redshift $z$, the halo mass
function $\displaystyle\derivative{n_{\textrm{halo}}}{M_{\textrm{h}}}$, and
the merger rate density per halo per logarithmic WD mass
$\displaystyle\derivative{R_{\textrm{WD-PBH}}}{\ln M_{\textrm{WD}}}$. We below
neglect the redshift dependence of the halo mass function and the merger rate
and use the current values for simplicity. We set the minimum and maximum
masses of halos at $M_{\mathrm{h},\mathrm{min}}=400M_{\odot}$ and
$M_{\mathrm{h},\mathrm{max}}=10^{16}M_{\odot}$, respectively. We have
confirmed that the event rate of WD-PBH mergers is not very sensitive to the
choice of $M_{\mathrm{h,min}}$ and $M_{\mathrm{h,max}}$.

Figure 6: Expected number of WD-PBH mergers detected by three-year observation
of DECIGO.

The detailed description of each term in the integrand is given in Appendix A.
Figure 6 shows the resultant event number in three years as a function of the
PBH mass. Note that we assume the fraction of PBH to the dark matter,
$f_{\textrm{PBH}}$, is unity though astrophysical observations have already
put the constraint $f_{\textrm{PBH}}\lesssim\mathcal{O}(10^{-1})$ for
$M_{\textrm{PBH}}=0.1M_{\odot}$. Figure 6 implies that the merger rate of WD-
PBH binaries is too low to detect them by three years of observation of
DECIGO.

The small event number is because of our conservative estimation. The expected
event number would be increased by several orders of magnitude, depending on
the assumptions on the number densities, the velocity dispersions, etc., of
the target objects. In the context of the NS-PBH merger, several possibilities
(e.g., formation in ultra-faint dwarf galaxies or globular clusters and PBH
clustering at their formation) provoking such enhancement are presented [35].
Before concluding this section, we briefly see how much enhancement is
required (and allowed) to detect the WD-PBH binaries. Let us parametrize the
enhancement of the event number by [35]

$\tilde{N}^{\textrm{WD-
PBH}}_{\textrm{event}}=\alpha_{\textrm{WD}}\times\alpha_{\textrm{PBH}}\times
f_{\textrm{PBH}}\times N^{\textrm{WD-PBH}}_{\textrm{event}}\,,$ (30)

where $\alpha_{\textrm{WD}}$ and $\alpha_{\textrm{PBH}}$ are the enhancement
factors induced by some mechanisms related to the WD and the PBH,
respectively.

Even forgetting about the realistic astrophysical processes,
$\alpha_{\mathrm{PBH}}$ cannot be enlarged arbitrarily because it would
conflict with the current constraint on the PBH-PBH merger rate achieved by
the ground-based GW detectors. As we explained in Sec. I, PBH-PBH binaries are
formed in two channels, and the early-time formation basically dominates the
late-time formation. However, if the enhancement factor
$\alpha_{\mathrm{PBH}}$ for WD-PBH is too large, it would also enhance the
late-time PBH-PBH formation, making it dominant and leading to the violation
of the current constraint on the event rate of subsolar mass PBH binary
coalescences. To discuss it quantitatively, we calculate the event number of
PBH-PBH mergers in the late-time formation scenario. Because
$\alpha_{\mathrm{PBH}}$ phenomenologically parametrizes the enhancement due to
the PBH properties, the PBH-PBH merger rate would also be enhanced in a
similar way as

$\tilde{N}^{\textrm{PBH-
PBH}}_{\textrm{event}}\sim(\alpha_{\textrm{PBH}}f_{\textrm{PBH}})^{2}N^{\textrm{PBH-
PBH}}_{\textrm{event}}$ (31)

where $N^{\textrm{PBH-PBH}}_{\textrm{event}}$ is the expected number of PBH-
PBH mergers estimated with $f_{\textrm{PBH}}=1$ and $\alpha_{\textrm{PBH}}=1$.
So far, binary mergers with masses of less than $1M_{\odot}$ have not been
detected. Therefore, $\tilde{N}^{\textrm{PBH-PBH}}_{\textrm{event}}$ should be
less than unity. For $M_{\textrm{PBH}}=10^{-1}M_{\odot}$ as our interest, we
find $N^{\textrm{PBH-PBH}}_{\textrm{event}}\sim 1\times 10^{-5}$ (see Fig. 10
in the appendix). Therefore, we get the upper bound

$\tilde{N}^{\textrm{PBH-PBH}}_{\textrm{event}}\lesssim
1\quad\Rightarrow\quad\alpha_{\textrm{PBH}}f_{\textrm{PBH}}\lesssim 3\times
10^{2}\,.$ (32)

Substituting this into $\tilde{N}^{\textrm{WD-PBH}}_{\textrm{event}}$ and
taking the value $N^{\textrm{WD-PBH}}_{\textrm{event}}=2\times 10^{-6}$, the
required enhancement factor $\alpha_{\textrm{WD}}$ for the detectable event
rate of WD-PBH binaries would be

$\tilde{N}^{\textrm{WD-PBH}}_{\textrm{event}}\gtrsim
1\quad\Rightarrow\quad\alpha_{\textrm{WD}}\gtrsim 2\times 10^{3}\,.$ (33)

If such an enhancement would be realized, one or more WD-PBH merger events
could be detected by DECIGO without contradicting the current constraint on
the PBH-PBH merger rate.

##  IV Conclusions and future perspective

We have investigated the detectability of WD and subsolar mass PBH binaries
formed by the GW bremsstrahlung. To this end, in Sec. II, we first
investigated the maximum distance of the detectable binaries as the horizon
distance assuming the simple waveform model given by Eqs. (7)–(9) and the
three-year observation of DECIGO with its design sensitivity. We have shown
that the GW from WD-PBH binaries at
$\sim$1\text{\,}\mathrm{G}\mathrm{p}\mathrm{c}$$ are detectable when their
mass satisfies $M_{\textrm{WD}}\gtrsim 0.8M_{\odot}$ and
$M_{\textrm{PBH}}\gtrsim 0.5M_{\odot}$. Second, we discussed the
distinguishability between WD-PBH and WD-WD mergers by estimating the merger
frequency from the comparison of their tidal and contact radii. We found that
WD-PBH mergers can be distinguished from other binaries just by looking at the
merger frequency unless the WD mass is less than $\sim 0.55M_{\odot}$.
Finally, We evaluated the merger rate and the expected number of WD-PBH merger
events in Sec. III. We assumed the standard halo properties and simplified WD
distributions. The conservative expectation of the number of WD-PBH merger
events is too small ($\sim\mathcal{O}(10^{-6})$) to be detected for
$M_{\textrm{PBH}}\sim 10^{-1}M_{\odot}$.

In a specific situation such as clustering distribution or low velocity-
dispersion of target objects, the expected number of WD-PBH merger events
would be much enhanced. We should note that an arbitrarily large enhancement
due to some mechanisms related to PBHs would contradict the current constraint
on the PBH-PBH merger rate in the ground-based detectors because the same
mechanisms would also promote the PBH-PBH merger. We conclude that
$\mathcal{O}(10^{3})$ enhancement due to WD-related mechanisms could enable
the detection of $\mathcal{O}(1)$ WD-PBH merger events by three-year DECIGO
observation without violating the current PBH-PBH constraint.

One may include additional effects to improve the detectability and
distinguishability of PBH-WD binaries. For example, the waveform modeling can
be improved by, e.g., the tidal deformability of WD. As a binary becomes
closer, the WD is gradually deformed, resulting in the modification in the
gravitational waveform. This modification is significant around the merger.
Another example is the effect of mass transfer. The WD-WD binary under the
stable mass transfer is discussed in the literature [46, 47]. It is known that
there is the possibility that the binary separation becomes wider as the mass
transfers from the smaller star to the heavier star (called _outspiral_). The
GW from such a binary are much different from that of the case without the
mass transfer. The synergy between GW and electromagnetic emissions can also
be beneficial in identifying the source. WD-WD binary mergers are candidates
of type Ia supernovae [48, 49]. WD-PBH binary mergers could emit
electromagnetic emissions that are different from those of WD-WD binary
mergers. On the analogy of NS-BH binaries, tidal disruption events can also be
expected from the WD-PBH binary mergers. Numerical simulations for these
merger events would be interesting (see [50] for a numerical simulation of a
binary coalescence of a NS and a subsolar mass PBH). As well as concrete
mechanisms for the merger rate enhancement discussed in Sec. III, these points
should be worth studying in future work.

###### Acknowledgements.

We thank Hiroki Takeda and Tomoya Kinugawa for the beneficial discussions.
This work is supported by JSPS KAKENHI Grants No. JP21K13918 (YT), JP20H01899,
JP23H04502 and JP23K13099 (TSY), JP20K03968 (SY), and JST SPRING Grant No.
JPMJSP2125 (RI).

##  Appendix A Details of merger rate calculation

###  A.1 Halo density profile

We employ the definition of the halo mass by the mass contained in a sphere
where the averaged density is $200$ times larger than the critical density
$\rho_{\textrm{crit}}$. The virial radius $r_{\textrm{vir}}$ is defined by the
radius of the sphere which is used in the definition of the halo mass. In
other words, $r_{\textrm{vir}}$ satisfies

$M_{\textrm{h}}=\frac{4\pi}{3}r^{3}_{\textrm{vir}}\times
200\rho_{\textrm{crit}}=4\pi\int_{0}^{r_{\textrm{vir}}}\differential r\
r^{2}\rho_{\textrm{DM}}(r)\,,$ (34)

where $\rho_{\textrm{DM}}(r)$ is the dark matter density profile. We use the
Navarro–Frenk–White (NFW) profile [51, 52] for dark matter mass density
profile (see Fig. 7 for examples of this profile),

$\rho_{\textrm{DM}}(r)=\frac{\rho_{s}}{(r/r_{\mathrm{s}})(1+r/r_{\mathrm{s}})^{2}}\,,$
(35)

where $\rho_{\mathrm{s}}$ and $r_{\mathrm{s}}$ characterize the density
profile and are determined from the halo mass $M_{\textrm{h}}$ and the
concentration factor $C\coloneqq r_{\textrm{vir}}/r_{\textrm{s}}$. The virial
radius is given by the halo mass as shown in Eq. (34). So, we get the
$r_{\textrm{s}}$ from $M_{\textrm{h}}$ and $C$. The density parameter
$\rho_{\textrm{s}}$ is given by

$\frac{\rho_{\mathrm{s}}}{\rho_{\textrm{crit}}}=\frac{200}{3}\frac{C^{3}}{\ln(1+C)-C/(1+C)}\,.$
(36)

Figure 7: Examples of NFW profiles. We fix the redshift to $z=0$.

Ludlow et al. [53] provides us with the fitting formula of the concentration
parameter as a function of the halo mass and the redshift,

$C(\nu)=c_{0}\left(\frac{\nu}{\nu_{0}}\right)^{-\gamma_{1}}\left[1+\left(\frac{\nu}{\nu_{0}}\right)^{1/\beta}\right]^{-\beta(\gamma_{2}-\gamma_{1})}\,,$
(37)

where $\nu(z)\coloneqq\delta_{\textrm{sc}}/\sigma_{\mathrm{m}}(M,z)$ is the
dimensionless peak height, $\delta_{\textrm{sc}}=1.686$ is the spherical top-
hat collapse threshold, and the parameters
$\\{c_{0},\beta,\gamma_{1},\gamma_{2},\nu_{0}\\}$ are

$\displaystyle c_{0}$ $\displaystyle=3.395\times x^{-0.215}\,,$ (38)
$\displaystyle\beta$ $\displaystyle=0.307\times x^{0.540}\,,$ (39)
$\displaystyle\gamma_{1}$ $\displaystyle=0.628\times x^{-0.047}\,,$ (40)
$\displaystyle\gamma_{2}$ $\displaystyle=0.317\times x^{-0.893}\,,$ (41)
$\displaystyle\nu_{0}$ $\displaystyle=(4.135-0.564x-0.210x^{2}$
$\displaystyle\quad+0.0557x^{3}-0.00348x^{4})\times D(z)^{-1}\,,$ (42)

with $x=1+z$ and the linear growth factor $D(z)^{-1}$. These relation is valid
for $0\leq\log(1+z)\leq 1$ and
$-8\leq\log_{10}\bqty{M/\pqty{h^{-1}M_{\odot}}}\leq 16.5$ . Here, $h$ is the
Hubble parameter in the unit of
$100\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}\,\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$.
The function $\sigma_{\mathrm{m}}(M,z)$ is the root mean square of the linear
density field smoothed on the scale
$R(M)=\bqty{3M/\pqty{4\pi\bar{\rho}_{\mathrm{m}}}}^{1/3}$ and defined by

$\sigma_{\mathrm{m}}^{2}(M_{\mathrm{h}},z)=\frac{1}{2\pi^{2}}\int^{\infty}_{0}\differential{k}k^{2}P(k,z)\absolutevalue{\hat{W}(kR(M_{\mathrm{h}}))}^{2}\,,$
(43)

where $P(k)$ stands for the linear matter power spectrum and $\hat{W}(kR)$ is
the Fourier transformation of the top-hat window function with radius $R$
defined in the real space. $\sigma_{\mathrm{m}}$ is approximated by

$\sigma_{\mathrm{m}}(M,z)=D(z)\frac{22.26\xi^{0.292}}{1+1.53\xi^{0.275}+3.36\xi^{0.198}}\,,$
(44)

where

$\xi\coloneqq\left(\frac{M_{\textrm{h}}}{10^{10}h^{-1}M_{\odot}}\right)^{-1}\,.$
(45)

The linear growth factor is approximated by

$D(z)=\frac{\Omega_{\textrm{m}}(z)}{\Omega_{\textrm{m,0}}}\frac{\Psi(0)}{\Psi(z)}(1+z)^{-1}\,,$
(46)

with

$\displaystyle\Psi(z)=\begin{multlined}\Omega_{\textrm{m}}(z)^{4/7}-\Omega_{\Lambda}(z)\\\
+\left(1+\frac{\Omega_{\textrm{m}}(z)}{2}\right)\left(1+\frac{\Omega_{\Lambda}(z)}{70}\right),\end{multlined}\Omega_{\textrm{m}}(z)^{4/7}-\Omega_{\Lambda}(z)\\\
+\left(1+\frac{\Omega_{\textrm{m}}(z)}{2}\right)\left(1+\frac{\Omega_{\Lambda}(z)}{70}\right),$
(47)
$\displaystyle\Omega_{\Lambda}(z)=\frac{\Omega_{\Lambda,0}}{\Omega_{\Lambda,0}+\Omega_{\textrm{m},0}(1+z)^{3}},$

and $\Omega_{\textrm{m}}(z)=1-\Omega_{\Lambda}(z)$.

###  A.2 Average of capture rate $\sigma v_{\mathrm{rel}}$

The cross-section $\sigma$ (not to confuse it with the perturbation amplitude
$\sigma_{\mathrm{m}}$ (43)) of a PBH and a WD with the relative velocity
$v_{\mathrm{rel}}$ is

$\sigma=2\pi\left(\frac{85\pi}{6\sqrt{2}}\right)^{2/7}\frac{G^{2}(M_{\textrm{PBH}}+M_{\textrm{WD}})^{10/7}M_{\textrm{WD}}^{2/7}M_{\textrm{PBH}}^{2/7}}{c^{10/7}v_{\textrm{rel}}^{18/7}}\,.$
(48)

The probability distribution of $v_{\mathrm{rel}}$ would be approximated by
the truncated Maxwell–Boltzmann distribution [54]:222Strictly speaking, Bird
et al. [54] employed this distribution for the PBH-PBH relative velocity.

$p(v_{\mathrm{rel}};v_{\mathrm{vir}},v_{\mathrm{DM}})\\\
=F_{0}\quantity[\mathrm{e}^{-v_{\mathrm{rel}}^{2}/v_{\mathrm{DM}}^{2}}-\mathrm{e}^{-v_{\mathrm{vir}}^{2}/v_{\mathrm{DM}}^{2}}]\Theta(v_{\mathrm{vir}}-v_{\mathrm{rel}})\,,$
(49)

where $v_{\mathrm{vir}}$ is the virial velocity, $v_{\mathrm{DM}}$ is the
velocity dispersion of the dark matter (PBHs), and $F_{0}$ is the constant
determined by the normalization

$1=4\pi\int^{v_{\mathrm{vir}}}_{0}\differential{v}v^{2}p(v;v_{\mathrm{vir}},v_{\mathrm{DM}})\,.$
(50)

$v_{\mathrm{DM}}$ and $v_{\mathrm{vir}}$ can be related with halo parameters
via Eq. (6) of Ref. [54]:

$v_{\mathrm{DM}}=\sqrt{\frac{GM(r<r_{\mathrm{max}})}{r_{\mathrm{max}}}}=\frac{v_{\mathrm{vir}}}{\sqrt{2}}\sqrt{\frac{C}{x_{\mathrm{max}}}\frac{g(x_{\mathrm{max}})}{g(C)}}\,,$
(51)

where $g(C)=\ln(1+C)-C/(1+C)$, and
$x_{\mathrm{max}}=r_{\mathrm{max}}/r_{s}=2.1626$ (equal to $C_{m}$ of Ref.
[54]). $M(r<r_{\textrm{max}})$ is the mass contained in the sphere of the
radius $r_{\textrm{max}}$. Using the NFW profile (35), we get

$\displaystyle M(r<r_{\textrm{max}})$
$\displaystyle=4\pi\int^{r_{\textrm{max}}}_{0}\differential r\
r^{2}\rho_{\textrm{DM}}(r)$ $\displaystyle=4\pi
r_{\textrm{s}}^{3}\rho_{\textrm{s}}\int^{x_{\textrm{max}}}_{0}\differential x\
\frac{x}{(1+x)^{2}}$ $\displaystyle=4\pi
r_{\textrm{s}}^{3}\rho_{\textrm{s}}\left[\ln(1+x_{\textrm{max}})-\frac{x_{\textrm{max}}}{1+x_{\textrm{max}}}\right]\,.$
(52)

Once we get the parameters $\\{r_{\mathrm{s}},\rho_{\mathrm{s}},C\\}$, we
obtain the velocity dispersion $v_{\mathrm{DM}}$ from the first equality of
Eq. (51) and the virial velocity $v_{\mathrm{vir}}$ from the second equality.
The average of $\sigma v_{\mathrm{rel}}$ is given by

$\langle\sigma
v_{\mathrm{rel}}\rangle=4\pi\int^{v_{\mathrm{vir}}}_{0}\differential{v}v^{2}p(v;v_{\mathrm{vir}},v_{\mathrm{DM}})\sigma
v\,.$ (53)

Using Eqs. (49) and (51), we get

$\displaystyle\langle\sigma v_{\textrm{rel}}\rangle$
$\displaystyle\simeq\mathsf{N}\times
2\pi\left(\frac{85\pi}{6\sqrt{2}}\right)^{2/7}$
$\displaystyle\quad\times\frac{G^{2}(M_{\textrm{PBH}}+M_{\textrm{WD}})^{10/7}M_{\textrm{WD}}^{2/7}M_{\textrm{PBH}}^{2/7}}{c^{10/7}v_{\textrm{DM}}^{11/7}}\,,$
(54)

where

$\mathsf{N}=\hat{F}_{0}\int^{\hat{v}_{\mathrm{vir}}}_{0}\differential{\hat{v}}\hat{v}^{3/7}\quantity(e^{-\hat{v}^{2}}-e^{-\hat{v}_{\mathrm{vir}}^{2}})\,,$
(55)

with $\hat{v}=v/v_{\mathrm{DM}}$,
$\hat{v}_{\mathrm{vir}}=v_{\mathrm{vir}}/v_{\mathrm{DM}}$, and
$\hat{F}_{0}=4\pi v_{\mathrm{DM}}^{3}F_{0}$.

###  A.3 Merger rate in a particular halo

The merger rate of a WD with the mass of $M_{\mathrm{WD}}$ and a PBH in a
particular halo with the mass of $M_{\textrm{h}}$ is

$\derivative{R_{\textrm{WD-PBH}}}{\ln
M_{\textrm{WD}}}=4\pi\int^{r_{\textrm{vir}}}_{0}\differential{r}r^{2}\derivative{n_{\mathrm{WD}}}{\ln
M_{\textrm{WD}}}n_{\mathrm{PBH}}\langle\sigma v_{\textrm{rel}}\rangle\,,$ (56)

where $n_{\mathrm{WD}}$ and $n_{\mathrm{PBH}}$ are the number density profiles
of WD and PBH, respectively.

For PBH, we simply assume that the density profile is proportional to the halo
mass density,

$\rho_{\textrm{PBH}}(r)=f_{\textrm{PBH}}\rho_{\textrm{DM}}(r)\,,$ (57)

and the mass spectrum of PBH is monochromatic,

$n_{\mathrm{PBH}}(r)=\frac{\rho_{\mathrm{PBH}}(r)}{M_{\mathrm{PBH}}}=f_{\textrm{PBH}}\frac{\rho_{\textrm{DM}}(r)}{M_{\mathrm{PBH}}}\,.$
(58)

Sasaki et al. [36] uses the exponential density profile for NS. We assume that
the number density profile of the WD also is an exponential form, i.e.,

${\derivative{n_{\mathrm{WD}}}{\ln
M_{\textrm{WD}}}}(r)=\mathcal{N}_{\mathrm{WD}}\mathrm{e}^{-r/r_{\mathrm{WD}}}\,,$
(59)

where $\mathcal{N}_{\mathrm{WD}}$ is the normalization constant depending on
the mass $M_{\textrm{WD}}$, and $r_{\mathrm{WD}}$ is the characteristic radius
that is assumed to be independent of the mass. We assume [36]

$r_{\mathrm{WD}}=0.1r_{\textrm{s}}\,.$ (60)

The normalization $\mathcal{N}_{\mathrm{WD}}$ is calculated as follows. First,
we start with the derivation of the galaxy’s stellar mass. We use the relation
between the halo mass $M_{\mathrm{h}}$ and the stellar mass $M_{\ast}$ of the
galaxy. Behroozi et al. [55] gives the fitting function

$\log_{10}[M_{\mathrm{h}}/M_{\odot}]=\log_{10}M_{1}+\beta\log_{10}\left(\frac{M_{\ast}}{M_{\ast,0}}\right)\\\
+\frac{(M_{\ast}/M_{\ast,0})^{\delta}}{1+(M_{\ast}/M_{\ast,0})^{-\gamma}}-\frac{1}{2}\,,$
(61)

with

$\displaystyle\log_{10}M_{1}$ $\displaystyle=M_{1,0}+M_{1,a}(a-1)\,,$ (62)
$\displaystyle\log_{10}M_{\ast,0}$
$\displaystyle=M_{\ast,0,0}+M_{\ast,0,a}(a-1)\,,$ (63) $\displaystyle\beta$
$\displaystyle=\beta_{0}+\beta_{a}(a-1)\,,$ (64) $\displaystyle\gamma$
$\displaystyle=\gamma_{0}+\gamma_{a}(a-1)\,,$ (65) $\displaystyle\delta$
$\displaystyle=\delta_{0}+\delta_{a}(a-1)\,,$ (66)

where $a$ is the scale factor. The parameters are taken from Table 2 of Ref.
[55]: $M_{1,0}=12.35$, $M_{1,a}=0.28$, $M_{\ast,0,0}=10.72$,
$M_{\ast,0,a}=0.55$, $\beta_{0}=0.44$, $\beta_{a}=0.18$, $\gamma_{0}=1.56$,
$\gamma_{a}=2.51$, $\delta_{0}=0.57$, and $\delta_{a}=0.17$. The galaxy’s
stellar mass determines the normalization constant of the mass spectrum of the
main sequence stars through the condition,

$\displaystyle M_{\ast}$
$\displaystyle=V\int^{M_{\textrm{MS,max}}}_{M_{\textrm{MS,min}}}\differential{M_{\textrm{MS}}}M_{\textrm{MS}}\derivative{n}{M_{\textrm{MS}}}$
$\displaystyle=V\int^{M_{\textrm{MS,max}}}_{M_{\textrm{MS,min}}}\differential{M_{\textrm{MS}}}\derivative{n}{\ln
M_{\textrm{MS}}}\,,$ (67)

where $V$ is defined by

$V=\frac{4\pi}{3}r_{\textrm{max}}^{3}\,,$ (68)

and we set $M_{\mathrm{MS,min}}=0.01M_{\odot}$ and
$M_{\mathrm{MS,max}}=100M_{\odot}$. We use the following mass spectrum [56],

$\displaystyle\frac{\mathrm{d}n}{\mathrm{d}\ln M_{\mathrm{MS}}}$
$\displaystyle=\begin{cases}A_{\mathrm{BD}}\left(\frac{M_{\mathrm{MS}}}{0.08M_{\odot}}\right)^{1-\alpha_{\mathrm{BD}}}&\left(0.01\leq
M_{\mathrm{MS}}/M_{\odot}\leq 0.08\right)\\\
A_{\mathrm{MS}}\left(\frac{M_{\mathrm{MS}}}{0.5M_{\odot}}\right)^{1-\alpha_{\mathrm{MS}1}}&\left(0.08\leq
M_{\mathrm{MS}}/M_{\odot}\leq 0.5\right)\\\
A_{\mathrm{MS}}\left(\frac{M_{\mathrm{MS}}}{0.5M_{\odot}}\right)^{1-\alpha_{\mathrm{MS}2}}&\left(M_{\mathrm{MS}}/M_{\odot}\geq
0.5\right)\end{cases}$ (69)

where $\alpha_{\textrm{BD}}=0.8$, $\alpha_{\textrm{MS1}}=1.3$, and
$\alpha_{\textrm{MS2}}=2.0$ [57]. The constant $A_{\textrm{MS}}$ is determined
by the continuity at $M_{\mathrm{MS}}=0.08M_{\odot}$, and $A_{\textrm{BD}}$ is
chosen so that Eq. (67) holds. We assume that the main sequence stars with
masses of $1M_{\odot}\leq M_{\mathrm{MS}}\leq 8M_{\odot}$ collapse to WD.
Taking into account the mass loss of the stars, the mass of a WD is estimated
by the empirical relation [58],

$\frac{M_{\textrm{WD}}}{M_{\odot}}=0.339+0.129\frac{M_{\textrm{MS}}}{M_{\odot}}\,,$
(70)

where $M_{\textrm{MS}}$ is a mass of the progenitor of the WD. Using this
relation, we can convert the mass spectrum of the main sequence stars to that
of the WD which is denoted by

$\derivative{\bar{n}_{\mathrm{WD}}}{\ln M_{\textrm{WD}}}\,,$

that has no radial dependence. This number density is understood as the
averaged one defined by

$\derivative{\bar{n}_{\mathrm{WD}}}{\ln
M_{\textrm{WD}}}=\frac{1}{V}\int_{0}^{r_{\textrm{max}}}\mathcal{N}_{\mathrm{WD}}\mathrm{e}^{-r/r_{\mathrm{WD}}}\times
4\pi r^{2}\differential{r}.$ (71)

The normalization constant is hence determined by

$\displaystyle\mathcal{N}_{\mathrm{WD}}$
$\displaystyle=\frac{V}{4\pi\int^{r_{\textrm{max}}}_{0}\differential{r}r^{2}\mathrm{e}^{-r/r_{\mathrm{WD}}}}\derivative{\bar{n}_{\mathrm{WD}}}{\ln
M_{\textrm{WD}}}$ $\displaystyle\simeq\frac{0.5V}{4\pi
r_{\mathrm{WD}}^{3}}\derivative{\bar{n}_{\mathrm{WD}}}{\ln
M_{\textrm{WD}}}\,.$ (72)

###  A.4 Merger rate per comoving volume and halo mass function

Figure 8: Halo mass function (74) in the form of
$\log_{10}\pqty{\frac{M_{\mathrm{h}}^{2}}{\bar{\rho}_{\mathrm{m}}}\derivative{n_{\mathrm{halo}}}{M_{\mathrm{h}}}}$.

Integrating the merger rate in a particular halo with a weight due to the halo
mass function (i.e., the number of halos with the mass of $M_{\textrm{h}}$ per
comoving volume), we get the merger rate of a particular mass of WD per unit
time per comoving volume by

$\int^{M_{\textrm{h,max}}}_{M_{\textrm{h,min}}}\differential{M_{\textrm{h}}}\derivative{n_{\textrm{halo}}}{M_{\textrm{h}}}\derivative{R_{\textrm{WD-
PBH}}(M_{\textrm{h}})}{\ln M_{\textrm{WD}}}\,,$ (73)

where $\derivative*{n_{\textrm{halo}}}{M_{\textrm{h}}}$ is the halo mass
function. The fitting function of the halo mass function is given by Tinker et
al. [59]:

$\derivative{n_{\textrm{halo}}}{M_{\mathrm{h}}}=f(\sigma_{\mathrm{m}})\frac{\bar{\rho}_{\mathrm{m}}}{M_{\mathrm{h}}}\derivative{\ln[\sigma_{\mathrm{m}}^{-1}]}{M_{\mathrm{h}}}\,,$
(74)

with the matter density
$\bar{\rho}_{\mathrm{m}}(z)=\Omega_{\mathrm{m}}(z)\rho_{\textrm{crit}}(z)$ and

$f(\sigma_{\mathrm{m}})=A\left[\left(\frac{\sigma_{\mathrm{m}}}{b}\right)^{-a}+1\right]\mathrm{e}^{-c/\sigma_{\mathrm{m}}^{2}}\,.$
(75)

For $\\{A,a,b,c\\}$, we use the values for $\Delta=200$ (corresponding to the
definition of the virial radius (34)) in Table 2 in Ref. [59], i.e.,

$A=0.186\,,\quad a=1.47\,,\quad b=2.57\,,\quad c=1.19\,.$ (76)

Figure 8 shows the halo mass function with different redshifts.

###  A.5 Rate of detectable mergers of PBH and WD with a particular mass
$M_{\mathrm{WD}}$

Figure 9: Unit comoving volume $\derivative*{V}{z}$ (77).

The unit comoving volume is defined by

$\derivative{V}{z}=\frac{c}{H(z)}4\pi\left[\int^{z}_{0}\frac{c\differential{z^{\prime}}}{H(z^{\prime})}\right]^{2}\,.$
(77)

Figure 9 shows $\derivative*{V}{z}$.

Using $\derivative*{V}{z}$, we get the rate of the detectable mergers of PBH
and WD with a particular mass $M_{\mathrm{WD}}$ by

$\int^{z_{\textrm{up}}}_{0}\differential{z}\derivative{V}{z}\int^{M_{\textrm{h,max}}}_{M_{\textrm{h,min}}}\differential{M_{\textrm{h}}}\derivative{n_{\textrm{halo}}}{M_{\textrm{h}}}\derivative{R_{\textrm{WD-
PBH}}(M_{\textrm{h}})}{\ln M_{\textrm{WD}}}\,.$ (78)

Here, $z_{\textrm{up}}$ is determined by the horizon distance,

$D_{\textrm{horizon}}=(1+z_{\textrm{up}})\int^{z_{\textrm{up}}}_{0}\frac{c\differential{z^{\prime}}}{H(z^{\prime})}\,.$
(79)

The horizon distance depends on the WD and PBH masses and the detector’s
sensitivity, and so does $z_{\textrm{up}}$.

###  A.6 Merger rate of PBH-PBH binaries formed in the late universe

To calculate the merger rate of PBH-PBH binaries, we follow the same process
as we did for PBH-WD binaries. The expected event number of PBH-PBH mergers is
given by

$N^{\textrm{PBH-PBH}}_{\textrm{event}}\\\
=T_{\textrm{obs}}\int^{z_{\textrm{up}}}_{0}\differential{z}\derivative{V}{z}\int^{M_{\textrm{h,max}}}_{M_{\textrm{h,min}}}\differential{M_{\textrm{h}}}\derivative{n_{\textrm{halo}}}{M_{\textrm{h}}}R_{\textrm{PBH-
PBH}}\,.$ (80)

The difference is the merger rate density per one halo. Because we assume the
monochromatic mass distribution for PBH, we do not carry out the integration
over $M_{\textrm{PBH}}$. The merger rate density per one halo $R_{\mathrm{PBH-
PBH}}$ is given by

$R_{\textrm{PBH-
PBH}}=4\pi\int^{r_{\textrm{vir}}}_{0}\differential{r}r^{2}f_{\textrm{PBH}}^{2}n_{\textrm{DM}}^{2}(r)\langle\sigma
v_{\textrm{rel}}\rangle\,.$ (81)

When we estimate $z_{\mathrm{up}}$, we use the PSD achieved in the first
observation run of the LIGO-Virgo Collaboration,

$S^{\mathrm{LIGO}}_{\mathrm{n}}(f)=\Bigg{[}\left(\frac{18\mathrm{Hz}}{0.1\mathrm{Hz}+f}\right)^{4}\times
10^{-44}+4.9\times 10^{-47}\\\ +1.6\times
10^{-45}\left(\frac{f}{2000\mathrm{Hz}}\right)^{2}\Bigg{]}\mathrm{Hz}^{-1}\,,$
(82)

which is taken from Ref. [60]. We assume the one-year observation. Figure 10
shows the event rate of PBH-PBH binary mergers detected by one-year
observation of the LIGO detector. The result is consistent with preceding
studies such as Bird et al. [54].

Figure 10: Expected number of PBH-PBH mergers detected by one-year observation
by LIGO.

## References

  * Zel’dovich and Novikov [1967] Y. B. Zel’dovich and I. D. Novikov, The Hypothesis of Cores Retarded during Expansion and the Hot Cosmological Model, Soviet Astron. AJ (Engl. Transl. ), 10, 602 (1967).
  * Hawking [1971] S. Hawking, Gravitationally Collapsed Objects of Very Low Mass, Monthly Notices of the Royal Astronomical Society 152, 75 (1971), https://academic.oup.com/mnras/article-pdf/152/1/75/9360899/mnras152-0075.pdf .
  * Carr and Hawking [1974] B. J. Carr and S. W. Hawking, Black holes in the early Universe, Mon. Not. Roy. Astron. Soc. 168, 399 (1974).
  * Carr [1975] B. J. Carr, The Primordial black hole mass spectrum, Astrophys. J. 201, 1 (1975).
  * Chapline [1975] G. F. Chapline, Cosmological effects of primordial black holes, Nature 253, 251 (1975).
  * Bean and Magueijo [2002] R. Bean and J. Magueijo, Could supermassive black holes be quintessential primordial black holes?, Phys. Rev. D 66, 063505 (2002), arXiv:astro-ph/0204486 .
  * Duechting [2004] N. Duechting, Supermassive black holes from primordial black hole seeds, Phys. Rev. D 70, 064015 (2004), arXiv:astro-ph/0406260 .
  * Liu and Bromm [2022] B. Liu and V. Bromm, Accelerating Early Massive Galaxy Formation with Primordial Black Holes, Astrophys. J. Lett. 937, L30 (2022), arXiv:2208.13178 [astro-ph.CO] .
  * Hütsi  _et al._ [2023] G. Hütsi, M. Raidal, J. Urrutia, V. Vaskonen, and H. Veermäe, Did JWST observe imprints of axion miniclusters or primordial black holes?, Phys. Rev. D 107, 043502 (2023), arXiv:2211.02651 [astro-ph.CO] .
  * Niikura  _et al._ [2019a] H. Niikura, M. Takada, S. Yokoyama, T. Sumi, and S. Masaki, Constraints on Earth-mass primordial black holes from OGLE 5-year microlensing events, Phys. Rev. D 99, 083503 (2019a), arXiv:1901.07120 [astro-ph.CO] .
  * Scholtz and Unwin [2020] J. Scholtz and J. Unwin, What if Planet 9 is a Primordial Black Hole?, Phys. Rev. Lett. 125, 051103 (2020), arXiv:1909.11090 [hep-ph] .
  * Witten [2020] E. Witten, Searching for a Black Hole in the Outer Solar System, arXiv:2004.14192 [astro-ph.EP] (2020).
  * Smirnov  _et al._ [2022] J. Smirnov, A. Goobar, T. Linden, and E. Mörtsell, White Dwarfs in Dwarf Spheroidal Galaxies: A New Class of Compact-Dark-Matter Detectors, arXiv:2211.00013 [astro-ph.CO] (2022).
  * Escrivà  _et al._ [2022] A. Escrivà, F. Kuhnel, and Y. Tada, Primordial Black Holes, arXiv:2211.05767 [astro-ph.CO] (2022).
  * Carr  _et al._ [2023] B. Carr, S. Clesse, J. Garcia-Bellido, M. Hawkins, and F. Kuhnel, Observational Evidence for Primordial Black Holes: A Positivist Perspective, arXiv:2306.03903 [astro-ph.CO] (2023).
  * Niikura  _et al._ [2019b] H. Niikura _et al._ , Microlensing constraints on primordial black holes with Subaru/HSC Andromeda observations, Nature Astron. 3, 524 (2019b), arXiv:1701.02151 [astro-ph.CO] .
  * Carr  _et al._ [2021] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Constraints on primordial black holes, Rept. Prog. Phys. 84, 116902 (2021), arXiv:2002.12778 [astro-ph.CO] .
  * Nakamura  _et al._ [1997] T. Nakamura, M. Sasaki, T. Tanaka, and K. S. Thorne, Gravitational waves from coalescing black hole MACHO binaries, Astrophys. J. Lett. 487, L139 (1997), arXiv:astro-ph/9708060 .
  * Ioka  _et al._ [1998] K. Ioka, T. Chiba, T. Tanaka, and T. Nakamura, Black hole binary formation in the expanding universe: Three body problem approximation, Phys. Rev. D 58, 063003 (1998), arXiv:astro-ph/9807018 .
  * Ali-Haïmoud  _et al._ [2017] Y. Ali-Haïmoud, E. D. Kovetz, and M. Kamionkowski, Merger rate of primordial black-hole binaries, Phys. Rev. D 96, 123523 (2017), arXiv:1709.06576 [astro-ph.CO] .
  * Abbott  _et al._ [2018] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Search for Subsolar-Mass Ultracompact Binaries in Advanced LIGO’s First Observing Run, Phys. Rev. Lett. 121, 231103 (2018), arXiv:1808.04771 [astro-ph.CO] .
  * Abbott  _et al._ [2019] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Search for Subsolar Mass Ultracompact Binaries in Advanced LIGO’s Second Observing Run, Phys. Rev. Lett. 123, 161102 (2019), arXiv:1904.08976 [astro-ph.CO] .
  * Abbott  _et al._ [2022a] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), Search for Subsolar-Mass Binaries in the First Half of Advanced LIGO’s and Advanced Virgo’s Third Observing Run, Phys. Rev. Lett. 129, 061104 (2022a), arXiv:2109.12197 [astro-ph.CO] .
  * Abbott  _et al._ [2022b] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), Search for subsolar-mass black hole binaries in the second part of Advanced LIGO’s and Advanced Virgo’s third observing run, arXiv:2212.01477 [astro-ph.HE] (2022b).
  * Phukon  _et al._ [2021] K. S. Phukon, G. Baltus, S. Caudill, S. Clesse, A. Depasse, M. Fays, H. Fong, S. J. Kapadia, R. Magee, and A. J. Tanasijczuk, The hunt for sub-solar primordial black holes in low mass ratio binaries is open, arXiv:2105.11449 [astro-ph.CO] (2021).
  * Morrás  _et al._ [2023] G. Morrás, J. F. Nuño Siles, J. García-Bellido, E. Ruiz Morales, A. Menéndez-Vázquez, C. Karathanasis, K. Martinovic, K. S. Phukon, S. Clesse, M. Martínez, and M. Sakellariadou, Analysis of a subsolar-mass compact binary candidate from the second observing run of Advanced LIGO, Physics of the Dark Universe 42, 101285 (2023), arXiv:2301.11619 [gr-qc] .
  * Prunier  _et al._ [2023] M. Prunier, G. Morrás, J. F. N. n. Siles, S. Clesse, J. García-Bellido, and E. Ruiz Morales, Analysis of the subsolar-mass black hole candidate SSM200308 from the second part of the third observing run of Advanced LIGO-Virgo, arXiv:2311.16085 [gr-qc] (2023).
  * Pujolas  _et al._ [2021] O. Pujolas, V. Vaskonen, and H. Veermäe, Prospects for probing gravitational waves from primordial black hole binaries, Phys. Rev. D 104, 083521 (2021), arXiv:2107.03379 [astro-ph.CO] .
  * Maggiore  _et al._ [2020] M. Maggiore _et al._ , Science Case for the Einstein Telescope, JCAP 03, 050, arXiv:1912.02622 [astro-ph.CO] .
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane _et al._ (LISA), Laser Interferometer Space Antenna, arXiv:1702.00786 [astro-ph.IM] (2017).
  * Seto  _et al._ [2001] N. Seto, S. Kawamura, and T. Nakamura, Possibility of direct measurement of the acceleration of the universe using 0.1-Hz band laser interferometer gravitational wave antenna in space, Phys. Rev. Lett. 87, 221103 (2001), arXiv:astro-ph/0108011 .
  * Kawamura  _et al._ [2021] S. Kawamura _et al._ , Current status of space gravitational wave antenna DECIGO and B-DECIGO, PTEP 2021, 05A105 (2021), arXiv:2006.13545 [gr-qc] .
  * Quinlan and Shapiro [1989] G. D. Quinlan and S. L. Shapiro, Dynamical Evolution of Dense Clusters of Compact Stars, Astrophys. J. 343, 725 (1989).
  * Mouri and Taniguchi [2002] H. Mouri and Y. Taniguchi, Runaway merging of black holes: analytical constraint on the timescale, Astrophys. J. Lett. 566, L17 (2002), arXiv:astro-ph/0201102 .
  * Tsai  _et al._ [2021] Y.-D. Tsai, A. Palmese, S. Profumo, and T. Jeltema, Is GW170817 a Multimessenger Neutron Star-Primordial Black Hole Merger?, JCAP 10, 019, arXiv:2007.03686 [astro-ph.HE] .
  * Sasaki  _et al._ [2022] M. Sasaki, V. Takhistov, V. Vardanyan, and Y.-l. Zhang, Establishing the Nonprimordial Origin of Black Hole–Neutron Star Mergers, Astrophys. J. 931, 2 (2022), arXiv:2110.09509 [astro-ph.CO] .
  * Cholis  _et al._ [2016] I. Cholis, E. D. Kovetz, Y. Ali-Haïmoud, S. Bird, M. Kamionkowski, J. B. Muñoz, and A. Raccanelli, Orbital eccentricities in primordial black hole binaries, Phys. Rev. D 94, 084013 (2016), arXiv:1606.07437 [astro-ph.HE] .
  * [38] http://www.astropy.org.
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Dalal  _et al._ [2006] N. Dalal, D. E. Holz, S. A. Hughes, and B. Jain, Short grb and binary black hole standard sirens as a probe of dark energy, Phys. Rev. D 74, 063006 (2006), arXiv:astro-ph/0601275 .
  * Yagi and Seto [2011] K. Yagi and N. Seto, Detector configuration of DECIGO/BBO and identification of cosmological neutron-star binaries, Phys. Rev. D 83, 044011 (2011), [Erratum: Phys.Rev.D 95, 109901 (2017)], arXiv:1101.3940 [astro-ph.CO] .
  * Magano  _et al._ [2017] D. M. N. Magano, J. M. A. Vilas Boas, and C. J. A. P. Martins, Current and Future White Dwarf Mass-radius Constraints on Varying Fundamental Couplings and Unification Scenarios, Phys. Rev. D 96, 083012 (2017), arXiv:1710.05828 [astro-ph.CO] .
  * Creighton and Anderson [2011] J. D. E. Creighton and W. G. Anderson,  _Gravitational-wave physics and astronomy: An introduction to theory, experiment and data analysis_ (2011).
  * Maggiore [2007] M. Maggiore,  _Gravitational Waves. Vol. 1: Theory and Experiments_ , Oxford Master Series in Physics (Oxford University Press, 2007).
  * Jaranowski and Krolak [2009] P. Jaranowski and A. Krolak,  _Analysis of gravitational-wave data_ (Cambridge Univ. Press, Cambridge, 2009).
  * Paczyński [1967] B. Paczyński, Gravitational Waves and the Evolution of Close Binaries, Acta Astron 17, 287 (1967).
  * Kremer  _et al._ [2017] K. Kremer, K. Breivik, S. L. Larson, and V. Kalogera, Accreting Double white dwarf binaries: Implications for LISA, Astrophys. J. 846, 95 (2017), arXiv:1707.01104 [astro-ph.HE] .
  * Iben and Tutukov [1984] I. Iben, Jr. and A. V. Tutukov, Supernovae of type I as end products of the evolution of binaries with components of moderate initial mass (M not greater than about 9 solar masses), Astrophys. J. Suppl. 54, 335 (1984).
  * Webbink [1984] R. F. Webbink, Double white dwarfs as progenitors of R Coronae Borealis stars and type I supernovae., Astrophys. J. 277, 355 (1984).
  * Markin  _et al._ [2023] I. Markin, A. Neuweiler, A. Abac, S. V. Chaurasia, M. Ujevic, M. Bulla, and T. Dietrich, General-relativistic hydrodynamics simulation of a neutron star–sub-solar-mass black hole merger, Phys. Rev. D 108, 064025 (2023), arXiv:2304.11642 [gr-qc] .
  * Navarro  _et al._ [1996] J. F. Navarro, C. S. Frenk, and S. D. M. White, The Structure of cold dark matter halos, Astrophys. J. 462, 563 (1996), arXiv:astro-ph/9508025 .
  * Navarro  _et al._ [1997] J. F. Navarro, C. S. Frenk, and S. D. M. White, A Universal density profile from hierarchical clustering, Astrophys. J. 490, 493 (1997), arXiv:astro-ph/9611107 .
  * Ludlow  _et al._ [2016] A. D. Ludlow, S. Bose, R. E. Angulo, L. Wang, W. A. Hellwing, J. F. Navarro, S. Cole, and C. S. Frenk, The mass–concentration–redshift relation of cold and warm dark matter haloes, Mon. Not. Roy. Astron. Soc. 460, 1214 (2016), arXiv:1601.02624 [astro-ph.CO] .
  * Bird  _et al._ [2016] S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, Did LIGO detect dark matter?, Phys. Rev. Lett. 116, 201301 (2016), arXiv:1603.00464 [astro-ph.CO] .
  * Behroozi  _et al._ [2010] P. S. Behroozi, C. Conroy, and R. H. Wechsler, A Comprehensive Analysis of Uncertainties Affecting the Stellar Mass-Halo Mass Relation for 0 $<$ z $<$ 4, Astrophys. J. 717, 379 (2010), arXiv:1001.0015 [astro-ph.CO] .
  * Kroupa [2001] P. Kroupa, On the variation of the initial mass function, Mon. Not. Roy. Astron. Soc. 322, 231 (2001), arXiv:astro-ph/0009005 .
  * Mróz  _et al._ [2017] P. Mróz, A. Udalski, J. Skowron, R. Poleski, S. Kozłowski, M. K. Szymański, I. Soszyński, Ł. Wyrzykowski, P. Pietrukowicz, K. Ulaczyk, D. Skowron, and M. Pawlak, No large population of unbound or wide-orbit Jupiter-mass planets, Nature (London) 548, 183 (2017), arXiv:1707.07634 [astro-ph.EP] .
  * Williams  _et al._ [2009] K. A. Williams, M. Bolte, and D. Koester, Probing The Lower Mass Limit for Supernova Progenitors and the High-Mass End of the Initial-Final Mass Relation from White Dwarfs in the Open Cluster M35 (NGC 2168), Astrophys. J. 693, 355 (2009), arXiv:0811.1577 [astro-ph] .
  * Tinker  _et al._ [2008] J. L. Tinker, A. V. Kravtsov, A. Klypin, K. Abazajian, M. S. Warren, G. Yepes, S. Gottlober, and D. E. Holz, Toward a halo mass function for precision cosmology: The Limits of universality, Astrophys. J. 688, 709 (2008), arXiv:0803.2706 [astro-ph] .
  * [60] Binary black hole signals in LIGO open data, https://github.com/losc-tutorial/LOSC_Event_tutorial/blob/master/LOSC_Event_tutorial.ipynb.

  *[PBH]: primordial black hole
  *[PBHs]: primordial black hole
  *[GW]: gravitational wave
  *[GWs]: gravitational wave
  *[WD]: white dwarf
  *[WDs]: white dwarf
  *[LVK]: LIGO–Virgo–KAGRA
  *[NS]: neutron star
  *[SNR]: signal-to-noise ratio
  *[PSD]: power spectral density
  *[NFW]: Navarro–Frenk–White