# On an existence problem of periodic points in intervals whose images cover
themselves

Wang Yihan
<EMAIL_ADDRESS>
(April 2022)

###### Abstract

We study a system of intervals $I_{1},\ldots,I_{k}$ on the real line whose
images under a continuous map $f$ contain themselves. It’s conjectured that
there exists a periodic point of period $\leq k$ in $I_{1}\cup\ldots\cup
I_{k}$. We prove the conjecture for $k=5$ in this paper. We also propose a
discretization method in attempt to solve the problem.

##  1 Introduction

We consider a continuous mapping $f$: $\mathbb{R}\to\mathbb{R}$, and $k$
closed intervals $I_{1},I_{2},...,I_{k}$. Throughout this paper, $f^{l}$
denotes the $l$-th iteration of $f$.

######  Definition 1.1.

We call $(I_{1},I_{2},\ldots,I_{k};f)$ a $covering\ system$ if $f(I_{1}\cup
I_{2}\cup\ldots\cup I_{k})\supseteq I_{1}\cup I_{2}\cup\ldots\cup I_{k}$.

This concept was first introduced by S. A. Bogatyi and E. T. Shavgulidze in
their paper [1], when they tried to obtain an analogue of Sharkovskii’s
theorem [2] for an arbitrary tree. They proved the existence of a periodic
point of period $\leq g(n)$ in a covering system and established an upper
bound for $g(n)$ $(g(n)\leq 2(n^{2}-1)g(n-1)^{n-2})$. Moreover, they proposed
the following conjecture in the remark.

######  Conjecture 1.2.

For any covering system $(I_{1},I_{2},\ldots,I_{k};f)$, there exists a point
$x_{0}\in I_{1}\cup I_{2}\cup\ldots\cup I_{k}$ such that $f^{l}(x_{0})=x_{0}$
for some $l\leq k$.

It is proved in [1] that Conjecture 1.2 is true for $k\leq 4$. By similar but
more detailed discussions, we can extend the results to $k=5$. And we obtain
the first main theorem of this paper.

######  Theorem 1.3.

For any covering system $(I_{1},\ldots,I_{5};f)$, there exists a periodic
point of period $\leq 5$ in $I_{1}\cup\ldots\cup I_{5}$.

Inspired by the process of proof for $k=5$, we develop new methods which
relate the initial problem to a discrete one. Indeed, we propose the following
conjecture and show that Conjecture 1.2 is a corollary of it.

######  Conjecture 1.4.

For any mapping $f$: $\\{1,2,\ldots,n\\}$ $\to$ subsets of
$\\{1,2,\ldots,n\\}$ such that $\bigcup_{i=1}^{n}f(i)=\\{1,2,\ldots,n\\}$, and
for any partition

$I_{1}=\\{1,2,\ldots,i_{1}\\},\quad
I_{2}=\\{i_{1}+1,\ldots,i_{2}\\},\quad\ldots,\quad
I_{k}=\\{i_{k-1}+1,\ldots,n\\},$

there exists $j\in\\{1,2,\ldots,k\\}$ and $r,s\in I_{j}$ such that

$(convf)^{l}(\\{r,s\\})\supseteq\\{r,s\\},\quad for\ some\quad l\leq k.$

Here, $convf(A)=$ convex hull of $f(A)=\\{\min f(A),\ldots,\max f(A)\\}$ for
any finite set $A\subseteq\mathbb{N}$.

######  Theorem 1.5.

Conjecture 1.4 implies Conjecture 1.2.

Moreover, we can prove that it suffices to consider the case $f$ is a
transitive permutation in the symmetry group $S_{n}$ in Conjecture 1.4. And we
will see Conjecture 1.4 is equivalent to the following one, which concerns
only the property of a permutation.

######  Definition 1.6.

Let $f$ be a transitive permutation in the symmetry group $S_{n}$, $i.e.$, $f$
can be written as $(i_{1}i_{2}\ldots i_{n})$, and $A_{i}=\\{i,i+1\\}$
$(i=1,2,\ldots,n-1)$ be $(n-1)$ particular sets. We define the
$characteristic\ number$ $m_{i}$ of each $A_{i}$ as

$m_{i}=\min\\{m\arrowvert\ (convf)^{m}(A_{i})\supseteq A_{i}\\}.$

Here, $convf(A)=$ convex hull of $f(A)=\\{\min f(A),\ldots,\max f(A)\\}$ for
any finite set $A\subseteq\mathbb{N}$.

And we define the $characteristic\ sequence$ of $f$ to be

$m_{1}^{\prime}\leq m_{2}^{\prime}\leq\ldots\leq m_{n-1}^{\prime},$

where $\\{m_{i}^{\prime}\\}$ is a rearrangement of $\\{m_{i}\\}$.

######  Conjecture 1.7.

For any transitive $f$ in the symmetry group $S_{n}$, the characteristic
sequence $\\{m_{i}^{\prime}\\}$ of $f$ satisfies $m_{i}^{\prime}\leq i$,
$i=1,2,\ldots,n-1$.

One of the advantages of Conjecture 1.7 is the independence of the partition.
Furthermore, using the directed graph associated to a periodic orbit (see [3]
page 7 for definitions), we can finally prove:

######  Theorem 1.8.

Conjecture 1.7 is equivalent to Conjecture 1.4.

We organize the paper in the following way.

In section 2 we explain our main idea of proof for $k=5$ through specific
examples. In fact, we claim that for $k$ small (at least for $k\leq 5$), the
existence of a periodic point only depends on the ordering of the end points
of the intervals and their images under $f^{i}(i\leq k)$. Therefore, Theorem
1.3 can be proved case by case.

In section 3 we give new ideas for attempts to prove Conjecture 1.2. We
explain in detail how we can reduce the initial problem to a discrete one. And
the main results, Theorem 1.5 and Theorem 1.8, are proved in this section.

##  2 The case $k=5$

The purpose of this section is to prove Theorem 1.3. Before coming to the
proof for $k=5$, we first take a look at some basic examples which help us
understand the problem. We will always assume the mappings under consideration
are continuous from now on.

######  Example 2.1 (the case $k=1$).

Let $f:I\to\mathbb{R}$ satisfy $f(I)\supseteq I$. Then there exists a fixed
point of $f$ in $I$.

######  Remark 2.2.

This fact is an easy consequence of the intermediate value theorem in
calculus. And we will see later that for $k$ larger, most cases can be reduced
to this fundamental one.

######  Example 2.3.

(See Figure 1) In this example, each interval $I_{i}$ is mapped to some
$f(I_{i_{j}})\supseteq I_{j}$. Then $f$ induces a permutation on
${I_{1},\ldots,I_{k}}$. Thus, $\exists$ $l\leq k$, $s.t.$
$f^{l}(I_{1})\supseteq I_{1}$. And by example 2.1, we can find $x_{0}\in
I_{1}$ with $f^{l}(x_{0})=x_{0}$.This example explains why we think the period
is not bigger than $k$ in conjecture 2.1.

Figure 1: Example 2.3

We state a well-known lemma which is useful when we want to simplify the
initial problem.

######  Lemma 2.4.

If $f([a,b])\supseteq[c,d]$, then there is a subinterval
$[a^{\prime},b^{\prime}]\subseteq[a,b]$ such that
$f([a^{\prime},b^{\prime}])=[c,d]$. And furthermore, the end points of
$[a^{\prime},b^{\prime}]$ is mapped to the end points of $[c,d]$.

###### Proof.

There are points $x_{c}$ and $x_{d}$ in $[a,b]$ such that $f(x_{c})=c$ and
$f(x_{d})=d$. Without loss of generality, we can assume $x_{c}<x_{d}$. Then
the proof is completed by choosing $a^{\prime}=\sup
f^{-1}(c)\cap[x_{c},x_{d}]$ and $b^{\prime}=\inf
f^{-1}(d)\cap[a^{\prime},x_{d}]$. ∎

######  Remark 2.5.

We see from this lemma that if necessary, we can replace $I_{i}$ in a covering
system $(I_{1},\ldots,I_{k};f)$ by some subintervals $I_{i}^{\prime}$ such
that $(I_{1}^{\prime},\ldots,I_{k}^{\prime};f)$ is still a covering system,
but the end points of each $I_{i}^{\prime}$ are mapped to the end points of
$f(I_{i}^{\prime})$.

Moreover, we can assume the initial covering system $(I_{1},\ldots,I_{k};f)$
is $minimal$, that is, there doesn’t exist any covering system
$(J_{1},\ldots,J_{l};f)$ with the same mapping $f$ such that $l\leq k$ and
each $J_{j}$ is contained in some $I_{i}$. (See lemma 1 in [1] for more
details.) In particular, we can assume that the interior sets of $f(I_{i})$
and $f(I_{j})$ are disjoint. Of course, we can also assume $I_{i}$ and $I_{j}$
are disjoint in Conjecture 1.2, since if not, we can consider a new interval
$I_{i}\cup I_{j}$ instead.

Next, we introduce our methods to deal with different cases for $k=5$ through
the following three examples. All the covering systems below are assumed to be
minimal. Hence, each end point of $I_{i}$ is mapped into $I_{1}\cup\ldots\cup
I_{k}$. (See [1] lemma 4.)

Figure 2: Example 2.6

######  Example 2.6.

The behavior of $f$ is illustrated in Figure 2.

As mentioned above, we assume the end points of $I_{1}$ are mapped to the end
points of $f(I_{1})$. We can assume further that $f(a)=c,\ f(b)=d$. Otherwise,
there will exist a fixed point of $f$ in $I_{1}$. For the same reason,
$f(x)>x,\forall x\in I_{1}$. So we obtain

$f^{2}(a)=f(c)>c>a,\quad f^{2}(b)=f(d)<b,$

which shows that $f^{2}(x_{0})=x_{0}$ for some $x_{0}\in I_{1}$.

We find that $I_{1}\cap f(I_{1})\neq\emptyset$, which is much special in this
example.

Figure 3: Example 2.7

######  Example 2.7.

See Figure 3 for the information of the covering system. And recall that we
always assume end points of $I_{i}$ are mapped to end points of $f(I_{i})$.

Case (1): For $I_{2}$, assume $f(d)=a,\ f(q)=c$. In this case, we consider the
interval $I_{5}$. We see that $d,e\in f(I_{5})$, and $f(d),f(e)\in
f^{2}(I_{5})$. Since $f$ is continuous,
$f^{2}(I_{5})\supseteq[f(e),f(d)]=[f(e),a]\supseteq[b,a]\supseteq I_{3}\cup
I_{4}$. Therefore, $f^{3}(I_{5})\supseteq f(I_{3}\cup I_{4})\supseteq I_{1}$.
And finally $f^{4}(I_{5})\supseteq f(I_{1})\supseteq I_{5}$. By example 2.1,
we know there is an $x_{0}\in I_{5}$ with $f^{4}(x_{0})=x_{0}$.

Case (2): For $I_{3}$, assume $f(e)=p,\ f(h)=b$. We argue in a similar way to
deduce $f^{5}(I_{5})\supseteq I_{5}$. And again we come back to example 2.1.

Case (3): None of the conditions in the above two cases are satisfied. Then
$f(d)=c,\ f(e)=b,\ f(q)=a,\ f(h)=p$. This time we consider the point $f(c)$.
And after similar discussions, we will obtain: $f^{4}(I_{2})\supseteq I_{2},$
if $f(c)\geq q$; $f^{5}(I_{5})\supseteq I_{5}$, if $f(c)<q.$ Anyway, there is
a periodic point of period $\leq 5$ in $I_{1}\cup...\cup I_{5}$.

In example 2.7, we no longer have the property mentioned in example 2.6.
However, we can study the image of some intervals and the behavior of their
end points under the mapping $f$ or even $f^{2}$ to reduce the problem to
example 2.1. Note that we frequently use the fact: if $f(I)$ contains $x,y\
(x<y)$ where $I$ is an interval, then $f(I)\supseteq[x,y]$.

######  Example 2.8.

See Figure 4.

Of course we can check in the same way as in example 2.6 and 2.7. However, we
can simplify the problem by using the induction trick. Indeed, in this
example, we can consider another covering system $(I_{1},I_{2};f^{2})$ instead
of the given one to reduce it to the case $k=2$.

Figure 4: Example 2.8

######  Remark 2.9.

According to the ordering of all the end points of $I_{i}$ and $f(I_{i})$, we
can reduce Conjecture 1.2 to finitely many cases when $k=5$. However, since
there are about $(5!)^{2}$ cases in total, it is impossible to list them all
in this paper. We claim here without proof that all the cases can be checked
using the methods in the above three examples. Thus Theorem 1.3 holds.

##  3 Discrete problems

In this section we show that conjecture 1.2 may be proved by solving another
discrete problem related to it. The main idea comes from our attempts to find
a general method when trying to prove Conjecture 1.2 for $k=5$. And we now
explain it.

As we have seen in the examples in section 2, it will be helpful for
simplification if the image of an interval contains entirely another one or
more intervals, like in example 2.3. However, not all mappings have this nice
property. Therefore, we have to discuss where the end points of each interval
go under the mapping $f$ or even $f^{2}$, just like what we have done in
example 2.7. So this inspires us to divide the initial intervals into smaller
subintervals, such that the image of each subinterval contains the whole of
some other subintervals, not just part of them.

For example, if we cut the intervals in example 2.7 into small pieces, then
more details of the mapping $f$ will be exhibited in the graph, besides the
behavior of the end points of each interval under $f$. (See Figure 5.)
However, it may happen that no matter how you divide the intervals, there will
be an interval whose image under $f$ contains only part of another interval.
What saves us here is that, we can consider perturbation of $f$ if the
intervals are cut into enough small pieces. Indeed, we have the following
lemma:

######  Lemma 3.1.

If there is a covering system $(I_{1},I_{2},\ldots,I_{k};f)$ such that

$\forall x\in I_{1}\cup\ldots\cup I_{k},\quad and\quad\forall l\leq k,\quad
f^{l}(x)\neq x,$

then $\exists$ $\delta>0$, $s.t.$, whenever $\lVert g-f\rVert<\delta$, it
holds that $g^{l}(x)\neq x$, for $l\leq k$, $x\in I_{1}\cup\ldots\cup I_{k}$.

###### Proof.

Recall that $f$ has a periodic point of period $\leq k$ in
$I_{1}\cup\ldots\cup I_{k}$ is equivalent to that
$F(x)=(f(x)-x)\ldots(f^{k}(x)-x)$ has a zero in $I_{1}\cup\ldots\cup I_{k}$.
Therefore, the proof is completed since $f$ is uniformly continuous on a
compact set containing $f^{k}(I_{1}\cup\ldots\cup I_{k})$. ∎

Now we can introduce the discrete problem related to Conjecture 1.2. Consider
the covering system in Lemma 3.1. Let
$a_{1}<b_{1}<a_{2}<b_{2}<\ldots<a_{k}<b_{k}$ be all the end points of the
intervals $I_{1},\ldots,I_{k}$ and denote
$M_{0}=\\{a_{1},b_{1},\ldots,a_{k},b_{k}\\}$. Consider their images
$f(M_{0})=\\{f(a_{1}),\ldots,f(b_{k})\\}$ and we eliminate the points which
don’t belong to $I_{1}\cup\ldots\cup I_{k}$ to obtain a set $S_{1}\subseteq
f(M_{0})$. Put $M_{1}=M_{0}\cup S_{1}$. Inductively, in each step we consider
the image of $M_{i}$ under the mapping $f$ and we eliminate those points
dropping out of $I_{1}\cup...\cup I_{k}$ to get a set $S_{i+1}\subseteq
f(M_{i})$. Then let $M_{i+1}=M_{i}\cup S_{i+1}$.

Finally we obtain a sequence of sets $M_{0}\subseteq
M_{1}\subseteq\ldots\subseteq M_{i}\subseteq\ldots$, where each $M_{i}$
contains the images of the initial end point set $M_{0}$ under
{$id,f,f^{2},\ldots,f^{i}$} except for those dropping out of
$I_{1}\cup\ldots\cup I_{k}$.

Then we can find a sufficiently large integer $N$ such that $\forall x\in
M_{N}-M_{N-1}$, $dist(x,M_{N-1})<\delta$, where $\delta$ comes from Lemma 3.1.
Therefore, we can take small perturbation of $f$ to get a map $\tilde{f}$ with
$\lVert\tilde{f}-f\rVert<\delta$, and $\widetilde{M_{N-1}}=\widetilde{M_{N}}$
for $\tilde{f}$. Indeed, we just move each point $(x_{0},f(x_{0}))\in
M_{N-1}\times M_{N}$ on the graph of $f$ to the nearest point
$(x_{0},y_{0})\in\\{x_{0}\\}\times M_{N-1}$ without destroying the continuity.
Thus, if we divide the intervals into small pieces using the points in
$\widetilde{M_{N}}$, then $\tilde{f}$ will have the nice property that the
image of each subinterval contains entirely another one or more subintervals.
In other words, $\tilde{f}$ can be regarded as a discrete mapping from
$\\{1,2,\ldots,n\\}$ to subsets of $\\{1,2,\ldots,n\\}$. (See Conjecture 1.4.)
Here each number $j$ represents a small piece of interval $J_{j}$. And
$j^{\prime}\in f(j)$ if $f(J_{j})\supseteq J_{j^{\prime}}$. Note that the
image of some interval may be empty, since we ignore the part outside
$I_{1}\cup\ldots\cup I_{k}$. And we only care about the interval determined by
$\\{f(a),f(b)\\}$ if $a,b\in\widetilde{M_{N}}$, although the image $f([a,b])$
may be larger.

It should also be remarked here that the only useful information of $f$ on the
so called “gaps” (see [1]) between two adjacent intervals $I_{i}$ and
$I_{i+1}$ is the behavior of the end points of each gap. In other words, for
the image $f(J_{j})$ of each small subinterval $J_{j}$, we only focus on the
part $f(J_{j})\cap(I_{1}\cup\ldots\cup I_{k})$. This is the reason why we
eliminate the points outside $I_{1}\cup\ldots\cup I_{k}$ in each step. We
believe this information is enough for the proof because in thousands of
examples for $k\leq 5$ it does suffice to get a conclusion.

To sum up, once we have a mapping $f$ satisfying the conditions in Lemma 3.1,
we can find an $\tilde{f}$ close enough to $f$, which also satisfies the
conditions but can be viewed as a discrete mapping. Such $f$ and $\tilde{f}$
are counterexamples to Conjecture 1.2. Therefore, if we argue by
contradiction, we can prove Theorem 1.5, which says that Conjecture 1.2. is a
corollary of Conjecture 1.4.

###### Proof of Theorem 1.5.

Suppose the covering system $(I_{1},I_{2},\ldots,I_{k};f)$ is a counterexample
to Conjecture 1.2. Then as mentioned above, after perturbation, we can assume
$f$ satisfies the conditions of Conjecture 1.4. Therefore the conclusion of
Conjecture 1.4 gives us a subinterval of $I_{j}$ represented by
$conv\\{r,s\\}$ whose image under $f^{l}$ contains itself. Thus there must
exist a periodic point of period $\leq l\leq k$ in $I_{j}$, a contradiction. ∎

Although Conjecture 1.4 reduces Conjecture 1.2 to a discrete problem, it is
not satisfactory since it involves something inconvenient to deal with, like
convex hull and partition. In the following we will simplify the conditions in
Conjecture 1.4 and find another description of it which only depends on the
own property of a permutation.

Let $f$ be in Conjecture 1.4. We can assume further:

(1)$\quad\forall i\neq j$, $f(i)\cap f(j)=\emptyset$.

(2)$\quad\forall i$, $f(i)\neq\emptyset$.

(3)$\quad\forall i$, $f(i)$ contains exactly one element.

(4)$\quad f$ is a transitive permutation in the symmetry group $S_{n}$,
$i.e.$, $f$ can be written as $(i_{1}i_{2}\ldots i_{n})$.

Indeed, (1) can be realized since we can move out the common part of $f(i)$
and $f(j)$ from one of them without changing the property
$\bigcup_{i=1}^{n}f(i)=\\{1,2,\ldots,n\\}$.

(2) is also satisfied because we can restrict the mapping $f$ to
$\\{1,2,\ldots,n\\}-\\{i\\}$ if necessary.

(3) is a consequence of (1), (2) and the condition
$\bigcup_{i=1}^{n}f(i)=\\{1,2,\ldots,n\\}$.

For (4), note that $f\in S_{n}$ because of (3). Suppose $f$ is not transitive.
Then choose an orbit $(i_{1}\ldots i_{m})$ of $f$ ($m<n$) and restrict both
$f$ and the partition $I_{1},\ldots,I_{k}$ to $\\{i_{1},\ldots,i_{m}\\}$.

In conclusion, we obtain:

######  Proposition 3.2.

It is sufficient to consider $f\in S_{n}$ and $f$ is transitive in Conjecture
1.4.

Figure 5: Example 3.3

######  Example 3.3.

We explain in this example how to obtain a discrete mapping in Conjecture 1.4.

For the covering system in Example 2.7 (Figure 3), we cut the initial
intervals into ten pieces:

$I_{1}=\\{1,2\\},\ I_{2}=\\{3,4\\},\ \ldots,\ I_{5}=\\{8,9,10\\}.$

And the image of each piece under $f$ happens to be some other pieces. (See
Figure 5.) Therefore $f$ can be regarded as a discrete mapping as follows:

$f(1)=\\{8,9\\},\ f(2)=\\{10\\},\ f(3)=7,\ f(4)=6,\ \ldots,\ f(9)=\emptyset,\
f(10)=\\{5\\}.$

Moreover, $f$ can be reduced to a transitive permutation in $S_{9}$ as $(1\to
8\to 4\to 6\to 2\to 10\to 5\to 3\to 7\to 1)$. And $I_{5}$ reduces to
$\\{8,10\\}$, respectively. Note here that if we follow the discussion at the
beginning of this section, there will be no need for a point between 8 and 9,
since it is not the image of the initial end points under any $f^{i}$. We add
this point in order to show how we reduce the mapping $f$ to a transitive
permutation in general cases.

Now we turn to the proof of the equivalence of Conjecture 1.7 and Conjecture
1.4. (See also Definition 1.6.)

######  Proposition 3.4.

Conjecture 1.7 implies Conjecture 1.4.

###### Proof.

As mentioned above, we can assume $f\in S_{n}$ is transitive, without loss of
generality.

Suppose Conjecture 1.7 is true. Then for any partition

$I_{1}=\\{1,2,\ldots,i_{1}\\},\quad
I_{2}=\\{i_{1}+1,\ldots,i_{2}\\},\quad\ldots,\quad
I_{k}=\\{i_{k-1}+1,...,n\\},$

the $(k-1)$ characteristic numbers
$\\{m_{i_{1}},m_{i_{2}},\ldots,m_{i_{k-1}}\\}$ cannot cover
$\\{m_{1}^{\prime},\ldots,m_{k}^{\prime}\\}$. Therefore, $\exists t\neq
i_{1},\ldots,i_{k-1}$, such that $m_{t}\leq k$, since
$m_{1}^{\prime}\leq\ldots\leq m_{k}^{\prime}\leq k$ in the characteristic
sequence. In other words,

$(convf)^{m_{t}}(A_{t})\supseteq A_{t},\quad where\quad
A_{t}=\\{t,t+1\\}\subseteq some\ I_{j}.$

Therefore, we can just choose $r=t,\ s=t+1$ in Conjecture 1.4.

∎

Recall that in one dimensional dynamic systems we often consider a directed
graph associated to a periodic orbit (For more details, see [4] or chapter 1
of [3]). Namely, if $f$ is a transitive permutation in Definition 1.6, we can
construct a directed graph $\Gamma_{f}$ with $(n-1)$ vertices
$\\{A_{1},A_{2},\ldots,A_{n-1}\\}$, where $A_{i}$ is defined in 1.6. And there
is an edge from $A_{i}$ to $A_{j}$ if and only if $convf(A_{i})\supseteq
A_{j}$. Thus, we obtain:

######  Proposition 3.5.

For a transitive $f\in S_{n}$ and its associated directed graph $\Gamma_{f}$,
the characteristic number $m_{i}$ is equal to the minimal length of cycles
starting from $A_{i}$ in $\Gamma_{f}$.

Using directed graphs, we can complete the proof of Theorem 1.8.

###### Proof of Theorem 1.8.

Suppose Conjecture 1.4 is true. For a transitive $f\in S_{n}$, suppose on the
contrary that there is a $k\leq n-1$ such that $m_{k}^{\prime}\geq k+1$, and
$k$ is the minimal number with this property. Then there are only $(k-1)$
characteristic numbers $m_{i_{1}},m_{i_{2}},\ldots,m_{i_{k-1}}\leq k$. Without
loss of generality, assume $i_{1}<\ldots<i_{k-1}$. And consider the partition:

$I_{1}=\\{1,2,\ldots,i_{1}\\},\quad
I_{2}=\\{i_{1}+1,\ldots,i_{2}\\},\quad\ldots,\quad
I_{k}=\\{i_{k-1}+1,\ldots,n\\}.$

By Conjecture 1.4, we can find $j,l\leq k$ and $r,s\in I_{j}$ such that

$(convf)^{l}(\\{r,s\\})\supseteq\\{r,s\\}.$

Now if we extend $f$ to be a piecewise linear mapping from the interval
$[1,n]$ to $[1,n]$, then we will obtain $f^{l}([r,s])\supseteq[r,s]$.
Therefore, there exists $x_{0}\in[r,s]$ with $f^{l}(x_{0})=x_{0}$. Note that
since $f$ is transitive as a permutation in $S_{n}$, $x_{0}$ cannot belong to
$\\{1,2,\ldots,n\\}$. Consequently, $x_{0}\in(t,t+1)$ for some integer
$t\in[r,s]$ and it is similar for
$f(x_{0}),f^{2}(x_{0}),\ldots,f^{l-1}(x_{0})$. Since $f$ is piecewise linear,
$x_{0}\in(t,t+1)$, $f(x_{0})\in(t^{\prime},t^{\prime}+1)$ imply
$f([t,t+1])\supseteq[t^{\prime},t^{\prime}+1]$, and it is similar for other
pairs of points. So the orbit
$\\{x_{0},f(x_{0}),\ldots,f^{l-1}(x_{0}),f^{l}(x_{0})\\}$ gives a cycle of
length $l\leq k$ starting from the vertex $A_{t}$ in the associated directed
graph $\Gamma_{f}$, which implies $m_{t}\leq k$. But this $t$ together with
$i_{1},\ldots,i_{k-1}$ gives $k$ characteristic numbers $\leq k$,
contradicting $m_{k}^{\prime}\geq k+1$.

Combined with Proposition 3.4, the proof is completed. ∎

######  Remark 3.6.

(1) The condition that $f$ is transitive cannot be removed in Conjecture 1.7.
Because if $f\in S_{3}$ and $f(1)=3,\ f(2)=2,\ f(3)=1$, then the
characteristic numbers will be $m_{1}=m_{2}=2$.

(2) $\\{m\arrowvert\ (convf)^{m}(A_{i})\supseteq A_{i}\\}$ cannot be replaced
by $\\{m\arrowvert(f^{m}(i)-i)(f^{m}(i+1)-(i+1))<0\\}$ in the definition of
the characteristic number in Conjecture 1.7. As an example, consider $f$:
$(1\to 2\to 4\to 3\to 1)$. Then the characteristic numbers will become
$m_{1}=m_{3}=3,\ m_{2}=1$, if we don’t take convex hull in each step.

We summarize what we have done. We show in Theorem 1.5 that the initial
problem on existence of periodic points in a covering system follows from a
discrete one. And we only need to focus on the transitive permutations by
Proposition 3.2. Finally, Conjecture 1.7 and Theorem 1.8 tell us that it is
sufficient to study the characteristic sequence of a transitive permutation,
which is a property depending only on the permutation itself.

At the end of this section, we give some examples of periodic orbits whose
characteristic sequences are known.

######  Example 3.7.

(1) Let $f\in S_{n}$ be $(1\to 2\to\ldots\to n\to 1)$. Then the characteristic
sequence of $f$ is $1\leq 2\leq\ldots\leq n-1$.

(2) Let $f\in S_{2n+1}$ be of Stefan type (see [5] for more details), that is,

$1\to(n+1)\to(n+2)\to n\to(n+3)\to(n-1)\to\ldots\to 2n\to 2\to(2n+1)\to 1.$

Then the characteristic sequence of $f$ is

$1\leq 2\leq 2\leq 4\leq 4\leq\ldots\leq 2n-2\leq 2n-2\leq 2n.$

(3) It has been checked by computer program that Conjecture 1.7 holds for all
$n\leq 12$. And in [6], the author gives all the directed graph associated to
transitive permutations in $S_{5}$.

## Acknowledgement

I would like to thank Prof. S.A. Bogatyi for introducing the problem to me.
And I thank my supervisor, Prof. Tian Gang, for his encouragement as well as
helpful suggestions. I’m also grateful to Dr. Ye Yanan and Dr. Zhao Yikai for
writing computer programs to verify some claims.

## References

  * [1] S. A. Bogatyi and E. T. Shavgulidze.  Periodic points of a map of a system of intervals.  Mathematical Notes of the Academy of Sciences of the Ussr, 43(3):210–219, 1988.
  * [2] A. N. Sharkovskii.  Coexistence of cycles of a continuous map of the line into itself.  International Journal of Bifurcation and Chaos, 5(5):1263–1273, 1995.
  * [3] L. S. Block and W. A. Coppel.  Dynamics in One Dimension.  Springer, 1992.
  * [4] Uhland Burkart.  Interval mapping graphs and periodic points of continuous function.  Journal of Combinatorial Theory, Series B, 32(1):57–68, 1982.
  * [5] P. Štefan.  A theorem of Šarkovskii on the existence of periodic orbits of continuous endomorphisms of the real line.  Communications in Mathematical Physics, 54(3):237–248, 1977.
  * [6] Philip D. Straffin.  Periodic points of continuous functions.  Mathematics Magazine, 51(2):99–105, 1978.


* [Ser58] Jean-Pierre Serre.  Sur la topologie des variétés algebriques en caractéristique p. Symposium Internacional de Topologia Algebraica  1958\.
  * [Ser62] Jean-Pierre Serre.  Corps locaux. Hermann  1962\.
  * [Sim96] Carlos Simpson.  The Hodge filtration on nonabelian cohomology.https://arxiv.org/abs/alg-geom/9604005 1996, https://arxiv.org/abs/alg-geom/9604005.
  * [SS03] Stefan Schwede, Brooke Shipley.  Equivalences of monoidal model categories.https://arxiv.org/abs/math/0209342 2003, https://arxiv.org/abs/math/0209342.
  * [SS13] Pierre Schapira, Jean-Pierre Schneiders.  Derived category of filtered objects.https://arxiv.org/abs/1306.1359 2013, https://arxiv.org/abs/1306.1359.
  * [To06] Bertrand Toën.  Affine stacks (Champs affines).https://arxiv.org/abs/math/0012219 2006, https://arxiv.org/abs/math/0012219.
  * [To06b] Bertrand Toen  Higher and derived stacks: a global overview.https://arxiv.org/abs/math/0604504 2006, https://arxiv.org/abs/math/0604504.
  * [To13] Bertrand Toen  Simplicial presheaves and derived algebraic geometry.https://hal.archives-ouvertes.fr/hal-00772850v2/document 2013, https://hal.archives-ouvertes.fr/hal-00772850v2/document.
  * [To20] Bertrand Toën.  Classes caractéristiques des schémas feuilletés.https://arxiv.org/abs/2008.10489 2020, https://arxiv.org/abs/2008.10489.
  * [TV06] Bertrand Toen, Gabriele Vezzosi  Homotopical Algebraic Geometry II: geometric stacks and applications.https://arxiv.org/abs/math/0404373 2006, https://arxiv.org/abs/math/0404373.
  * [TV09] Bertrand Toën and Gabriele Vezzosi.  $S^{1}$-Equivariant simplicial algebras and de Rham theory.https://arxiv.org/abs/0904.3256 2009, https://arxiv.org/abs/0904.3256.
  * [TV20] Bertrand Toën and Gabriele Vezzosi.  Algebraic foliations and derived geometry: the Riemann-Hilbert correspondence.https://arxiv.org/abs/2001.05450 2020, https://arxiv.org/abs/2001.05450.
  * [TV20b] Bertrand Toën and Gabriele Vezzosi.  Algebraic foliations and derived geometry II: the Grothendieck-Riemann-Roch theorem.https://arxiv.org/pdf/2007.09251.pdf 2020, https://arxiv.org/pdf/2007.09251.pdf.


# Impact of signal clusters in wide-band searches for continuous gravitational
waves

Lorenzo Pierini1,2, Pia Astone2, Cristiano Palomba2, Aidan Nyquist3, Simone
Dall’Osso 2, Sabrina D’Antonio4, Sergio Frasca2, Iuri La Rosa1,2,5, Paola
Leaci1,2, Federico Muciaccia2, Ornella J. Piccinni6, Luca Rei7.

1Università di Roma La Sapienza, Roma, Italy;
2INFN, Sezione di Roma, Roma, Italy;
3North Park University, Chicago, Illinois, USA;
4INFN, Sezione di Roma Tor Vergata, Roma, Italy;
5University of Savoia Mont Blanc and CNRS/IN2P3 Annecy France;
6Institut de Física d’Altes Energies (IFAE), Barcelona Institute of Science
and Technology, E-08193 Barcelona, Spain;
7INFN, Sezione di Genova, Genova, Italy.

###### Abstract

In this paper we present a study of some relevant steps of the hierarchical
frequency-Hough (FH) pipeline, used within the LIGO and Virgo Collaborations
for wide-parameter space searches of continuous gravitational waves (CWs)
emitted, for instance, by spinning neutron stars (NSs). Because of their weak
expected amplitudes, CWs have not been still detected so far. These steps,
namely the spectral estimation, the peakmap construction and the procedure to
select candidates in the FH plane, are critical as they contribute to
determine the final search sensitivity. Here, we are interested in
investigating their behavior in the (presently quite) extreme case of signal
clusters, due to many and strong CW sources, emitting gravitational waves
(GWs) within a small (i.e. $<1$ Hz wide) frequency range. This could happen
for some kinds of CW sources detectable by next generation detectors, like
LISA, Einstein Telescope and Cosmic Explorer. Moreover, this possibility has
been recently raised even for current Earth-based detectors, in some scenarios
of CW emission from ultralight boson clouds around stellar mass black holes
(BHs). We quantitatively evaluate the robustness of the FH analysis procedure,
designed to minimize the loss of single CW signals, under the unusual
situation of signal clusters. Results depend mainly on how strong in amplitude
and dense in frequency the signals are, and on the range of frequency they
cover. We show that indeed a small sensitivity loss may happen in presence of
a very high mean signal density affecting a frequency range of the order of
one Hertz, while when the signal cluster covers a frequency range of one tenth
of Hertz, or less, we may actually have a sensitivity gain. Overall, we
demonstrate the FH to be robust even in presence of moderate-to-large signal
clusters.

##  I Introduction

Searches for continuous gravitational waves (CWs), emitted for instance by
rapidly rotating neutron stars (NSs), either isolated or binary systems (see
ref:lasky ; ref:CWreview  for reviews of sources and possible emission
mechanisms), are among the highest priorities of the LIGO-Virgo-KAGRA
Collaboration ref:WP . The third scientific run of the LIGO-Virgo advanced
detectors (O3) took place from April 1, 2019 to the end of March 2020. CW
signals have a duration that is longer than the typical run time of the
detectors and contain important information on the nature and dynamics of the
source. Presently, no CW source has been detected, due to the weakness of this
kind of GW emission, compared to that of compact binaries coalescence
ref:GWTC1 ; ref:GWTC2 ; ref:GWTC2.1 ; ref:GWTC3 .

CW sources we are looking for fall into five broad categories ref:WP : (1)
nonaccreting known pulsars for which timing data are available; (2) other
known or suspected isolated neutron stars, with known sky location but limited
or absent timing information; (3) unknown isolated NSs in any direction; (4)
accreting stars in known/unknown binary systems; (5) long-transients (with
duration of hours/days/months), due for example to a postmerger newborn NS, to
r-mode instabilities or to newborn highly magnetized NSs (magnetars)
ref:GW170817PM ; ref:LtransO3 . Moreover, recently a growing interest has been
focused on boson clouds forming around spinning BHs, as potential CW sources
ref:bosonTEOR . The data analysis method used is different for each of these
targets. For unknown sources, i.e. when the source parameters (namely,
position in the sky, GW frequency, frequency time variation) are not known,
fully coherent searches are computationally unfeasible. For this reason,
hierarchical procedures (which alternate coherent and incoherent analysis
steps) have been set up. Recent results on all-sky CW searches using data of
the second and third scientific run of the LIGO and Virgo detectors (O2 and
O3) have been reported in ref:allskyO2 , ref:allskyO3a , ref:allskyO3 . A
short review of the main LIGO-Virgo all-sky CW search methods is also
presented in ref:CWmethods .

One of the standard search methods for all-sky searches is the frequency-Hough
(FH) pipeline, described in detail in ref:FH . The FH pipeline is a
hierarchical procedure that aims at identifying the most significant CW
candidates in a wide parameter space. A follow-up procedure, described in
ref:allskyO2 , is then applied to these candidates, in order to confirm or
reject them. Recently, a GPU implementation of the FH was run to analyze the
highest, and computationally demanding, part of the O3 all-sky search, showing
a computational gain of more than one order of magnitude IuriGPU . The FH
hierarchical procedure strongly reduces the computational cost with respect to
a fully coherent search, but at the price of a sensitivity loss. In
particular, some thresholds and selection criteria are applied at various
stages of the analysis. Any possible signal lost during one of these steps
will not be recoverable by later steps. For this reason, each step of the
procedure is designed to reduce as much as possible the loss of possible CW
signals, which could be hidden by the presence of strong disturbances or even
by other superposing signals.

In this paper we review the procedures for spectral estimation, peakmap (PM)
construction and candidate selection applied on recent LIGO-Virgo data, with
the aim of testing their effectiveness when clusters of signals are present.
All these steps play a relevant role in determining the signal detectability
and may bring to unwanted sensitivity losses if not properly addressed. In
particular, by adding simulated signals to O2 data, we demonstrate the
robustness of the FH pipeline in presence of signal clusters, quantifying the
very qualitative claims reported in ref:papaetc .

The paper is organized as follows. In Sec. II we briefly remind the main steps
of the FH pipeline. In Sec. III we discuss the main issue this paper deals
with, i.e., the possibility that some steps of the FH analysis are negatively
affected by the presence of several concurrent large CW signals confined in a
very small frequency range. Section IV is devoted to describe in detail two
relevant analysis steps, i.e., spectral estimation and construction of the
time-frequency PM, which could be affected by the presence of dense clusters
of signals. Sections V and VI are the core of the paper and are focused on the
detailed studies we have done, by means of software simulated signals added to
O2 data, to quantify the effect on the detection efficiency of dense clusters
of signals. We end with a discussion on the results in Sec. VIII.

##  II The FH search procedure

A scheme of the FH pipeline is given in figure 1 and is briefly described in
the following.

Figure 1: Scheme of the hierarchical FH procedure, in the case of data from 2
detectors. See Sec. II for a short description of the steps of the procedure.

  * •

[A] Creation of the short fast Fourier transform (FFT) database (SFDB).

Detector calibrated data are split into chunks of duration 8192, 4096, 2048
and 1024 seconds. Each chunk is then Fourier transformed, using the FFT
algorithm, obtaining four different sets of frequency domain data covering,
respectively, the range [0-128] Hz, [128-512] Hz, [512-1024] Hz, and
[1024-2048] Hz. At each FFT chunk are also associated additional data, such as
position and velocity of the detector at the chosen reference time. There is a
work in progress to make this step in a more flexible way, making use of the
“band sampled data” (BSD) BSD  framework.

  * •

[B] Autoregressive (AR) spectral estimation.

An AR algorithm, described in ref:pia_sfdb , is used to smooth the noise
spectral density of each chunk of data in such a way to preserve possible CW
signals, which in a single FFT would be confined within one single frequency
bin, and at the same time to be able to follow slower spectral variations. The
filtering is usually done backwards through the frequency bins, that is from
higher to smaller frequencies, to better adapt recursion to the noise behavior
of the detectors, which shows much higher and rapidly varying values in the
low frequency band, i.e., below $\approx$ 30 Hz.

  * •

[C] PM creation.

For each FFT, local maxima above a given threshold, called peaks, are selected
in the ratio $\cal{R}$ among the periodogram (square modulus of the FFT) and
the AR spectrum estimation. The PM is the collection of all the peaks, each
identified by a value of the frequency and a time (the middle time of the
corresponding FFT), both measured at the detector. AR spectral estimation and
PM construction have been used since many years for these analyses and their
mathematics (and statistics) is described in ref:pia_sfdb .

  * •

[D] FH transform.

For each sky position, the peaks in the PM are shifted in frequency to remove
the Doppler effect for that position at the peak time. The shifted peaks are
fed to the Hough transform, which maps the time frequency peaks to the source
frequency and spin-down plane.

  * •

[E] Candidate selection and coincidence analysis.

Candidates are selected in the FH plane by finding, for each sky position, the
points with highest number count in each frequency and spin-down subinterval.
Candidates from the analysis of a given data set are clustered together, and
cross-checked against candidates found in the analysis of another data set (of
the same detector or of a different detector) by means of coincidences in the
signal parameter space.

  * •

[F] Candidates follow-up and verification.

A deeper analysis is done on those candidates which survive the coincidence
step, in order to check their astrophysical origin or demonstrate they are
compatible with noise or due to detector artifacts. The standard follow-up
approach consists in rerunning the search starting from longer duration FFTs,
in a limited portion of the parameter space around each candidate. In case of
no detection, upper limits (UL) on the signal strain amplitude are computed as
a function of the frequency.

##  III Boson clouds constraints from FH search on O2 data

Recent results obtained for the O2 all-sky CW search ref:UL  have also been
used to derive constraints for a class of CWs emitters different from spinning
NSs, namely “clouds” of ultralight bosons around black holes (BHs). Such
clouds are expected to emit CWs at a frequency which depends at first order on
the boson mass and at second order on the product of the boson mass and BH
mass (as explained in ref:UL  and references therein). Even if more specific
methods have been developed to search for this kind of signals, see for
example ref:nostroBOS ; ref:bosonO3 , it is possible to map the results of a
standard all-sky CW search into exclusions limits in a plane defined by the
mass of the scalar boson field and the mass of the BH, as discussed in ref:UL
. The basic assumption is that all the steps of the hierarchical FH procedure,
especially designed and tested for the weak and rare CW signals emittted by
isolated NSs, are compatible with the characteristics of the boson cloud
population. In particular, we need to verify if the sensitivity is degraded
when the emitted signals are clustered in frequency. Such scenario could arise
under optimistic assumptions about the number of stellar mass BHs present in
our Galaxy ref:BHpop1 ; ref:BHpop2 : if we assume that the boson clouds
formation mechanism affects a significant fraction of galactic BHs, we could
have a large number of concurrent CW emitters. These signals would have proper
frequency centered at a value corresponding to the boson mass (which
reasonably takes just one value), with a spread due to the BH mass
distribution. If many signals are concentrated in a small frequency range and
the signals are “strong” enough111There are presently estimations which give
quite different results, depending on the basic assumptions. See for example
ref:BHpop1 ; ref:BHpop2 ., it could happen that the stronger signals hide the
weaker ones or that there is some degradation of sensitivity due to the fact
that signals mix, making it difficult to separately identify them. Such
possible complication might happen especially in future GW detectors, like
Cosmic Explorer ref:futureCE  and the Einstein Telescope ref:futureET , which
will be the first pan-European ground-based GW antenna ref:adele . These 3G
ground-based detectors will achieve one order of magnitude better
sensitivities with respect to the existing detectors. On the other side, the
space-based detector LISA ref:futureL  will give us access to the mHz
frequency band.

In a recent publication ref:papaetc , the impact of signal clusters on the
spectral estimation used in the FH approach has been very qualitatively
discussed. Specifically, plots suggesting the AR procedure cancels signal
peaks have been shown, when the signals are concentrated in such a way to
produce bumps in the detector noise, thus reducing the search sensitivity. We
will demonstrate these conclusions are wrong in most cases.

We quantify these qualitative predictions, studying how spectral estimation,
PM construction and candidate selection behave in this situation. In fact,
these are the steps of the procedure that might suffer from problems related
to the presence of many strong signals. We inject fake CW signals into O2 data
to mimic this effect and to see if, and at which level of density and strength
of the signals, we might need to apply modifications to the procedure. For
simplicity, this study is done using one week of data, which is enough to
evaluate possible problems and sensitivity losses. In fact, by analyzing a
longer data set the discriminatory power of the Doppler modulation, see Sec.
V, would be even stronger making our robustness tests conservative.

##  IV Spectral AR estimation and PM construction

The first important step with an impact on the sensitivity of the hierarchical
FH procedure is the estimation of the average power spectral density. The
procedure is described in details in ref:pia_sfdb  and we briefly remind the
reader here its basic aspects, and give some examples using O2 data. A nearly
monochromatic wave with enough high amplitude will produce a delta spike above
the noise floor in the periodogram. That peak will be included into the PM
only if its power is greater than the corresponding point in the AR estimation
by a fixed factor. So, if the AR estimator follow that peak, its chances to be
selected in the PM decrease. So, a good spectral estimator for CW searches
should have the following properties:

  * •

If narrow peaks in the frequency domain are present, the estimator should not
be affected by that peaks. This should be as much as possible independent on
the signal-to-noise ratio (SNR) of the peak.

  * •

If the noise level varies, either slowly or rapidly, the estimator should be
able to follow the noise variations.

Let $x_{i}$ be the data samples of the FFT. The amplitude spectral density is
estimated from an AR estimation, as shown in the equations below:

$y_{i}=x_{i}+w\cdot y_{i-1}$ $w=e^{-\delta\nu/\tau_{f}}$

where $y_{i}$ are the samples of the not normalized AR mean, obtained using
$w$ as weight, with $\delta\nu$ being the FFT resolution and $\tau_{f}$ the
memory of the AR mean (with dimensions of a frequency). The values of the
normalized AR estimations are given by

$\mu_{i}=\frac{y_{i}}{Z_{i}},$

where the normalization constant is $Z_{i}=1+w\cdot Z_{i-1}$. Besides this, in
order to have a “clean” estimator, that is not affected by spectral peaks, we
define a threshold $V_{max}$ and an age $A_{max}$. While
$r=\frac{x_{i}}{\mu_{i-1}}$ is lower than (or equal to) the threshold, the new
datum $x_{i}$ is used to evaluate the actual mean and the age of the estimator
is set to zero (expressed in number of samples).
When $r=\frac{x_{i}}{\mu_{i-1}}$ is larger than the threshold, the new datum
is not used to evaluate the actual mean and the age of the estimator is
incremented by 1 bin (that is by a frequency equal to the resolution). This
eliminates or at least reduces the effect of peaks from the estimation.
If the estimator becomes too old, i.e., if the age becomes greater than the
maximum age we have set, we deduce that we are not in presence of a peak (and
thus of a possible signal), but only that the noise characteristics have
changed. We thus go back by a number of samples $n=A$, and begin a new
evaluation of the mean, restarting from zero at the sample $i-A$. This is
needed to deal with all those situations when the noise is highly non-
stationary. Figure 2 shows an example, where the AR estimation is plotted
together with the absolute value of the corresponding FFT. From the plots on
the left and from the zoom on the right it is possible to note how, when
narrow peaks are present in the FFT, they do not appear, or appear
significantly reduced, in the AR estimation, and so are not (or only slightly)
suppressed in the PM.

Figure 2: An example of one FFT (green dots, blue lines), of duration 8192 s,
done using LIGO Livingston O2 data ref:GWOSC , and the corresponding AR
amplitude spectral density estimation (red). Here, for comparison, the FFT has
been normalized to represent the amplitude spectral density and we show the
absolute value. The left plot covers the band [5-512] Hz and the second plot
is a zoom in the frequency region from 11 Hz up to 23 Hz. It is possible to
appreciate that the AR estimation follows sharp changes on noise floor, like
those around 14 Hz and 17 Hz. On the other side, narrow peaks like those at
13.15 Hz, 14.85 Hz, 17.25 Hz and 19.65 Hz do not appear in the AR estimation,
as required.

More examples are given in Sec. V, where we show the results after adding fake
CW signals to the O2 data.

##  V Signal clusters: preliminary considerations

As anticipated, we want to study the robustness of our method. In the
following, we will focus on the AR spectral estimation and PM characterization
steps, showing how they behave in presence of many CWs.

###  V.1 Qualitative effect on AR estimation

In what follows, we have generated and injected several CW signals, with
randomly chosen sky localization, nearby source frequencies and, for
simplicity, without spin-down. The range of signal amplitudes we have
considered is $2\times 10^{-25}-10^{-24}$, which is above the ULs found in O2
ref:UL  at the chosen frequencies, and corresponding to the range predicted in
ref:papaetc  for the detectable galactic boson clouds/BHs population. Please
note that the signal amplitude is well above the noise level: if the signals
were weaker, the cumulative effect on AR estimation and on the overall
procedure would be even lower.
Figure 3 shows what happens to the AR estimation inside a single FFT, in 4
different injection configurations. The FFTs are built using 4096 seconds of
LIGO Livingston O2 data ref:GWOSC , so their frequency resolution (as well as
for the AR estimation) is $\delta\nu=1/4096\rm s\simeq 2.44\cdot 10^{-4}~{}\rm
Hz$. All added signals have their proper frequencies $f_{0}$ within few
frequency bins, but because of their randomly distributed sky position they
are further spread by Doppler modulation. In details, the configurations are
in Table 1.

label | n.sig |  $f_{0}$ range [Hz] | $h_{0}$
---|---|---|---
1 | 10 | $\left[338.5,~{}333.5+\delta\nu\right]$ | $1\times 10^{-24}$
2 | 20 | $\left[338.5,~{}333.5+2\delta\nu\right]$ | $1\times 10^{-24}$
3 | 50 | $\left[338.5,~{}333.5+5\delta\nu\right]$ | $1\times 10^{-24}$
4 | 50 | $\left[338.5,~{}333.5+5\delta\nu\right]$ | $6\times 10^{-25}$
Table 1: Parameters used to simulate CWs in different density configurations,
labeled from 1 to 4.

Figure 3: Plots of the 4 different injections listed in Table 1 made on one
FFT calculated using 4096 s of LIGO Livingston O2 data ref:GWOSC . Top left,
top right, bottom left, bottom right show the configurations respectively of
label 1, 2, 3 and 4.Blue lines and black dots represent the absolute value of
one normalized FFT, in a 3 Hz band, around 338.5 Hz. The red line is the
corresponding AR amplitude spectral density estimation. The effect of the
signal peaks in the AR estimation is negligible in the first two cases, very
small in the last case, while it reduces the contribution of some of the
strongest peaks in the third case (50 strong signals, with source frequency
$f_{0}$ within 5 bins).

By looking at the AR estimation (red curve in Fig. 3), it is clear that
generally it is negligibly affected by the presence of an ensemble of signals
at low densities: the only case when it is affected is when we add 50 strong
signals with $f_{0}$ within 5 bins, which means within 1.23 mHz. In this case,
the ratio $\cal{R}$ obtained on the central peaks is reduced by $\sim$30%.
This fact obviously motivates a further investigation. However, as it will be
shown in the next subsection, this amplitude reduction does not automatically
affect the sensitivity of the search, as

  * •

it is not affecting all the peaks and, most of all, not always the same peaks
when their time evolution is considered;

  * •

the amplitude of the peaks is never used in the analysis, so what matters is
only the limit posed by the threshold used in the PM creation.

These aspects are discussed in more detail in the following subsection.

###  V.2 Peakmap characterization

As anticipated, it is important not to lose possible signal peaks at this
level, as they cannot be recovered after. We also notice that in order to
build the PM, first the ratio $\cal{R}$ of the periodogram (square modulus of
the FFT) to the AR power spectrum estimation is computed. Then, local maxima
of $\cal{R}$ above a given threshold are selected. In fact the actual value of
$\cal{R}$ for the selected peaks is not used in the analysis. This choice
allows us to reduce the impact of large noise spectral disturbances.
A very relevant point is the role of Doppler correction. This is done by
properly shifting the peaks frequency in the input PM, for each sky position
for which we are running the search. As shown in Figure 4, even in the case of
multiple signals belonging to the same frequency bin at a given time, the
Doppler correction, which is different for each sky position, permits to
clearly distinguish the various signals. This is a very powerful feature as,
moreover, it enhances the SNR of CW signals with respect to noise lines. The
first plot in Figure 4 shows a zoom of a PM with 10 injected signals. The
signals have been generated with the same intrinsic frequency $f_{0}$, that is
within the same frequency bin, but with randomly chosen sky positions, so that
their frequencies at the detector appear separated, due to the Doppler effect
which introduces a maximum frequency shift of

$\Delta f_{dop}\simeq 2f_{0}\left|\frac{\vec{v}\cdot\hat{n}}{c}\right|\simeq
2f_{0}\times 10^{-4},$ (1)

where $\vec{v}$ is the detector velocity vector and $\hat{n}$ is the versor
identifying the source position in the sky. When the Doppler effect is removed
correctly for each of these signals, we can see the injected signal as a
sequence of peaks along a straight line at the right $f_{0}$ frequency. This
is shown from the second to the fourth plots of Figure 4. It is possible to
note here that there is only one straight line, which corresponds to the
properly Doppler corrected peaks.

Figure 4: The top left plot shows a zoom on a PM of LIGO Livingston O2 data
ref:GWOSC , where different fake CW signals are added. Those signals have all
the same proper frequency, $f_{0}=338.5$ Hz, but different sky positions. They
are strong enough to produce patterns visible by eye: each one follows a
different time-frequency trajectory, having different Doppler modulations. In
the other plots (top right, bottom left and bottom right) the Doppler
correction is applied respectively for one of the signals sky direction. When
the Doppler correction matches the sky localization of a signal source, its
corresponding trajectory becomes a straight line confined on $f_{0}$. On the
other side, the other signals trajectories get further distortion.

In conclusion, these results show that, even if signals are clustered in
frequency, the presence of the Doppler modulation gives a powerful tool to
distinguish them. Thus, in all analyses where candidates are selected
separately for each sky position, as in the case of the FH pipeline, the
presence of signal clusters does not blind the search. Clearly, the resolving
power is limited by the resolution of the sky grid used, which depends on the
frequency resolution. The candidate separation effect is clearly stronger as
the observing time improves, and it is maximum for one year of data, when the
Doppler annual modulation completes one cycle.

##  VI Efficiency study: Frequency-Hough can handle multiple signals

This section is devoted to properly quantify the effect of multiple signals
and their impact on the pipeline performances. In order to do this, we need to
define an observable. We consider signals that, after the Doppler correction
are monochromatic, meaning they are represented in the PM by a straight line
at the constant emission frequency $f_{0}$. So, after the proper corrections,
we construct an histogram of the number of peaks as a function of the
frequency, using the same frequency resolution of the search. We then expect
to find an excess of counts in the histogram bin associated to the source
frequency $f_{0}$. We define the efficiency as

$\eta=\frac{n_{sig}-\bar{n}_{noise}}{n_{FFT}},$ (2)

where $n_{sig}$ is the histogram number count at the frequency $f_{0}$,
$\bar{n}_{noise}$ is the average number of counts at all the other frequencies
except for $f_{0}$, and $n_{FFT}$ is the number of FFTs used to construct the
PM. Figure 5 shows an example of efficiency versus frequency, having injected
a signal at $f_{0}=338.5$ Hz, constructed the PM and removed the Doppler
effect.

Figure 5: Example of efficiency calculation from full PM after Doppler
correction. A histogram of the peaks is built using the frequency resolution
as bin width. Here an excess of counts is clearly visible at the frequency
$338.5$ Hz, where a signal has been injected into $n_{FFT}=1000$ discrete
FFTs. The red line represents the theoretically calculated expected noise
counts, $E[\bar{n}_{noise}]=p_{0}\cdot n_{FFT}$.

We need to underline that the detection efficiency $\eta$ defined here does
not determine the full sensitivity of the FH procedure, even if the two
concepts are strictly related. The expected value of $\eta$ is the difference
between the probability of selecting a peak when a signal with normalized
spectral amplitude $\lambda$ is present, namely $p_{\lambda}$, and the
probability of selecting a peak when only Gaussian noise is present, that is
$p_{0}$ ref:FH :

$E[\eta]=p_{\lambda}-p_{0}.$ (3)

In the following, we will estimate $\eta$ in different situations. We use one
week of data222We use only one week of data in order to limit the
computational cost of the study, as the goal of the paper is a comparison of
different situations on the same data set, selected from the LIGO Livingston
O2 run ref:GWOSC , starting from GPS=1186606177. The efficiency is computed as
a function of the density of signals in given frequency intervals. We call
$\Delta f_{0}$ the interval in which the source-frame frequencies are randomly
generated. Because of the random sky localization of the sources, the
detector-frame observed frequency will be unevenly distributed into a
frequency range given by

$\Delta f=\Delta f_{0}\pm\Delta f_{dop}.$ (4)

where $\Delta f_{dop}$ is the spread of the signal frequency at the detector
frame due to the Doppler effect, given by Eq. 1. In this way we can establish
how many signals to inject in a given frequency interval by tuning $f_{0}$ and
$\Delta f_{0}$. The number of frequency bins covered by the signals is

$n_{\mathrm{bin}}=\frac{\Delta f}{\delta\nu}$ (5)

and we can refer to the resulting mean density of signals-per-frequency-bin as

$\rho_{sig}=\frac{N_{\mathrm{sig}}}{n_{\mathrm{bin}}}.$ (6)

In the following subsections we will quantify the results in terms of
estimates of detection efficiency, first adding signals with the same
amplitude and then with uniformly distributed amplitudes.

###  VI.1 Constant amplitude signals

As a first case study, we consider a rather extreme scenario. In what follows
50 signals, all with the same very large amplitude, are simulated with proper
frequency in a band respectively $5\delta\nu$ and 10 Hz wide. The injection
parameters are listed in Table 2. Figure 6 shows all the obtained efficiencies
in both cases, together with mean values, median and standard deviation of the
results. Estimations of the overall efficiency, also with the signal-per-bin
densities, are reported in Table 2.

label | $\Delta f_{0}$ | $h_{0}$ | $\rho_{sig}$ | $E\left[\eta\right]\pm\sigma_{\eta}$
---|---|---|---|---
$5\delta\nu$ | $\left[338.5,338.5+5\delta\nu\right]$ | $1\times 10^{-24}$ | $\sim 0.15$ | $0.58\pm 0.09$
10Hz | $\left[333.5,343.5\right]$ | $1\times 10^{-24}$ | $\sim 10^{-3}$ | $0.72\pm 0.09$
Table 2: Parameters used to simulate CWs in two different density
configurations, labeled as $5\delta\nu$ and 10Hz, the obtained signal
densities $\rho_{sig}$ and efficiency estimations.

The mean efficiency loss (for the case of nearby signals with respect to the
case of widely spaced signals) is $\sim 19\%$. In other words, in that extreme
situation we recover in average the $81\%$ of pixels per signal with respect
to the widely spaced one. In Sec. VII much higher densities than
$\rho_{sig}\sim 0.15$ will be explored, but very different efficiency losses
will be obtained. The reason is that the constant signal amplitudes well above
the noise floor strongly impact the AR estimation. Anyway, we note that the
presence of only strong signals seems to be so unrealistic in the context of
CWs from boson clouds. We could have a similar scenario in future signal-
dominated detectors, like LISA. In the following, we simulate signals with
randomly distributed amplitudes.

Figure 6: Efficiency $\eta$ in the two cases of: 50 injected signals with
$h_{0}=10^{-24}$ distributed in $\Delta f_{0}=5\delta\nu\sim 1.2$ mHz - left
plot - and widely spaced (over 10 Hz) - right plot. Mean, median and standard
deviation of the efficiency are also shown.

###  VI.2 Uniform amplitude signals with low densities

In this subsection we present the results of multiple injections in a more
realistic regime in which signal amplitudes are not constant, as in the
previous case, but uniformly distributed. In order to get a complete
understanding, we cover the whole frequency region where the detectors have
the best sensitivity, [70 - 512] Hz. The signals have been simulated with
source frequencies around 12 different values in that interval and with source
positions randomly chosen. For each one of that 12 central frequencies, 10
signals have been injected with proper frequencies distributed in six
different ranges $\Delta f_{0}$. The choice is done in order to search for
evidence of a decreasing efficiency when the signals are concentrated in a
smaller region. All the chosen parameters are reported below:

$h_{0}\in\left[2\times 10^{-25},1\times 10^{-24}\right]$
---

$f_{0}$ $\in$ | {70.5; | 105.5; | 120.5; | 170.5; | 205.5; | 220.5; |
---|---|---|---|---|---|---|---
| 270.5; | 330.5; | 380.5; | 420.5; | 442.5; | 492.5} | Hz
$\Delta f_{0}\in$ | {$\delta\nu$; |  $10\delta\nu$; |  $20\delta\nu$; |  $50\delta\nu$; |  $100\delta\nu$; | 10 Hz} |

Thus, the signal densities $\rho_{\mathrm{sig}}$ explored in this section are
in the range [0.001 - 0.1]. We can refer to this part as a ”low density”
regime.

Figure 7 shows bar plots of median values of the obtained efficiencies. They
are grouped according to the different $\Delta f_{0}$ (represented by
different colors), for the different injection frequencies. The efficiency
values do not show any systematic decrease when going toward the smaller
$\Delta f_{0}$, instead they seem more random. In Figure 8 the bar plots are
combined together and the average efficiencies (with respect to the injection
frequency), together with standard deviations, are shown as a function of
$\Delta f_{0}$. The highest standard deviation in the case of signals injected
in 10 Hz is due to the fact that within a so large band it is more likely to
find strong narrow disturbances that worsen the detection efficiency in that
band. Also in this plot, there is no clear effect on efficiencies due to the
clustering of source frequencies. We can conclude that in the low density
regime there is not efficiency reduction and the procedure is able to recover
all the signals.

Figure 7: Bar plots of the simulation done in Sec. VI.2. The bars show the
median efficiency computed for 10 signals with source frequency within
respectively N$\delta\nu$ bins (N=1, 10, 20, 50, 100) and 10 Hz
(40960$\delta\nu$ bins). The color code is: 10 signals in 1 (red) 10 (green)
20 (cyan) 50 (yellow) 100 (magenta) bins and 10 Hz (black). On the bottom of
each box the central source frequency of the injected signals is shown. Signal
amplitudes have been generated from a uniform distribution in the range
$h_{0}=2\times 10^{-25}-1\times 10^{-24}$. Those plots show that there is not
an evident efficiency loss when signals are concentrated in smaller $\Delta
f_{0}$. Figure 8: Error plot of the simulation done in Sec. VI.2, in the “low-
density” regime. Differently from Fig. 7, efficiencies are combined together
grouping for the different 12 frequency injections. Average efficiency $\eta$
(stars) and standard deviation (bars) are shown as function of $\Delta f_{0}.$

##  VII Simulation of signals clusters from boson clouds

In this section, we simulate situations that are expected in the most extreme
cases of CWs emitted by boson clouds around galactic BHs. In the simulations
proposed in ref:papaetc , it has been suggested that the most extreme BHs -
namely, those with spins close to the maximum allowed value - could
potentially produce up to thousands of signals with amplitude above the O2 ULs
at the detectors. Depending on galactic BH population and the boson mass,
those signals could cluster into a frequency region having a width that spans
from $\sim 0.01\mathrm{Hz}$ to $\sim 1\mathrm{Hz}$. Based on these
indications, we reproduced two scenarios, namely one with signals clustered
into a [0.04 - 0.06] Hz band and another with signals in an order of magnitude
larger $0.8$ Hz band. Also in these cases we repeated the simulations
centering on different frequencies in the range [70 , 512] Hz. Over 4,000
signals were generated with random parameters and amplitude in the range
$\left[1\,-20\right]\times 10^{-25}$, following a power law distribution such
to mimic the ratio of detectable signals below and above the UL as found in
ref:papaetc .

In Figure 9 the results of the simulation are shown. The plotted quantity is
the ratio between the detection efficiency of signals in the clustered
configuration and the efficiency of the same signals, when injected alone.
This is done in order to compare efficiency variations at different
frequencies and on signals with different amplitudes, which have different
detection efficiencies. This efficiency ratio is plotted as function of the
mean signal density. Colored areas span from minimum to maximum values of the
mean ratio obtained from injections around different frequencies, in the range
[100,400] Hz, in clean bands. The blue area represents the ratio for the whole
ensemble of injected signals, whereas the green and red belts refer to the
subsets of signals with $h_{0}$ respectively above and below O2 ULs. Looking
at the different regimes, we recognize qualitatively different behaviors.

1) In the case of signals injected into a $\sim 0.06$ Hz band (Figure 9, left
plot), for signals with amplitude above O2 ULs, the efficiency loss ratio is
of few percents at the most. After a minimum of $\approx 0.9$ reached between
$\rho_{sig}\approx 1-2$, the ratio grows up to 1. On the other side, for
signals with amplitude below O2 ULs the efficiency ratio increases for growing
signal densities, which means their detection efficiency increases - i.e. we
recover more pixels belonging to those signals. On the whole set of signals,
the result is a global efficiency gain in the range $\left[1,2.2\right]$, the
highest reached at the highest signal densities. An analogous reinforcement
effect is used in ref:5nvect ; ref:HC2F , where the contributions of multiple
subthreshold CWs are combined to enhance their detection chances in the
context of targeted searches. Consequently, in this configuration there is no
overall efficiency loss and, moreover, also weak signals have a growing
probability to be recovered. Since also stronger signals do not suffer a
significant efficiency loss, we conclude that both detection efficiency and
evaluation of ULs are not affected.

2) The case of signals injected into a $\sim 0.8$ Hz wide band (Figure 9,
right plot) corresponds to the most extreme situation depicted in ref:papaetc
. In this case the behavior of the efficiency ratio changes significantly.
First of all, the two distinct signal subsets have no opposite dynamics as
they show both an efficiency loss as the signals density increases. While the
efficiency ratio of signals with amplitude below O2 ULs stops decreasing at a
signal density $\rho_{sig}\approx 1$ (with a corresponding efficiency ratio of
$\sim$ 0.8 $\pm$ 0.05), the efficiency ratio for the subset of signals with
amplitude above the UL continues decreasing, reaching values around 0.6 for
signal densities $\rho_{sig}\approx 2$. For the whole set of signals, values
in the range 0.75 $\pm$ 0.05 are reached. Since the detection efficiency is
directly related to the individual probability $p_{\lambda}$ to select a peak
if a signal is present, and given that the overall search sensitivity is
proportional to the square root of that probability ref:FH , we conclude that
in the most extreme - and probably very unrealistic - cases the search
procedure could lose up to $\sim 15\%$ of the “optimal” search sensitivity
(i.e. that with no superposition of signals). As discussed in the following, a
mean signal density of two is much larger than predicted in ref:papaetc .

Figure 9:  Injections in the “high-density” regime. Mean values of the ratio
between detection efficiencies in the case of signals clustered together
($\mathrm{Eff}$) and in the case of each of those signals taken alone
($\mathrm{Eff}_{0}$), as function of the signals-per-bin density (bottom x
axis) and number of injected signals (top x axis). The signals are injected
around 170 Hz, 240 Hz, 338 Hz, 380 Hz: colored areas span from minimum to
maximum values obtained at the different frequencies. The blue area shows the
overall efficiency ratio for the whole signal ensemble, whereas green and red
areas refer to signals that are respectively above and below O2 ULs. Left plot
shows the result when signals are clustered in a wide 0.06 Hz frequency range.
Right plot shows the result when signals are clustered in a wide 0.8 Hz
frequency range.

Both the described dynamics have a clear explanation in the way the AR
spectral estimator works. When the signal density is such to produce a wide
excess power in frequency, the answer of the AR estimator depends on how wide
the bump is with respect to the AR memory $\tau_{f}$ defined in Section IV.
The chosen value for AR memory is typically fixed for each frequency band. In
the present case, in the frequency range [128-512] Hz, it is used
$\tau_{f}=0.02$ Hz. In case 1), the signal cluster covers a wide $0.06$ Hz
frequency range, which is the same order of $\tau_{f}$. Consequently, AR
estimate is not fast enough to adapt to the increased noise (plus signal)
level. When $\rho_{sig}\geq 1$, every frequency bin has an expected occupation
of 1 signal. After that, any added signal would superpose to signals already
present in the same bin, thus increasing its power content. The result is that
the signal power in all bins of interest increases with respect to the noise
floor. If the AR estimation is not fast enough to adapt, it will be below the
most part of bins power, so that also weaker signals get more detection
chances. This explains the observed growing efficiency in case 1). On the
other hand, in case 2) the signal cluster covers a wide $0.8$ Hz frequency
range, which is one order of magnitude bigger than $\tau_{f}$. In this case,
even if signals continue to accumulate in the same frequency bins at
$\rho_{sig}\geq 1$ and their power content increases, the AR estimation is
able to adapt to the changed level. Thus, weak signals remain below the AR
level, while strong signals are weakened with respect to the previous noise
level.

We note two important considerations. First, the loss happens when we have
actually high numbers of detectable signals (i.e., above the upper limit),
while it affect much less the estimation of ULs. Second, we plot the
efficiency loss as a function of the mean signal density. For instance, a
density of two signals per bin, means that on average in all bins of the
injection band (e.g., 0.8 Hz in 9, right plot) we have such density. Looking
at the results of ref:papaetc , in one of their “worst” situations - see
ref:papaetc  figure 40, left panel - where BH spins up to an unrealistically
large value of 1 are considered, the signal density above the UL can reach
values as large as 5 in a very few bins, while in the vast majority of bins is
below 2 (and, actually, is zero in most bins). The average density is
difficult to estimate, but it is quite likely well below 0.5. As a
consequence, the sensitivity loss is likely of the order or smaller than
5$\%$, as the impact on the AR estimation would be smaller.

##  VIII Discussion

In this paper, mainly inspired by the claims of ref:papaetc , we have studied
in detail how the basic steps of the FH algorithm for CW searches could be
affected by the presence of clusters of signals. Two regimes, of “low-density”
(with up to $\sim 0.1$ signals per frequency bin) and “high-density” signals
(with up to 10 signals per frequency bin), have been considered. In particular
we have evaluated the impact of signal clusters in the estimation of the
average spectrum, and on the construction of the PMs, finding them very
robust. We have found that the procedure used to estimate the average
spectrum, needed to normalize the FFTs and then to construct the PMs, is
robust with respect to the presence of ensemble of signals and works even in
the most extreme situations. The small signal peak amplitude loss is
compensated by the fact that signals do have a different evolution in time
(due to their position in the sky) and then through the proper Doppler
correction each signal can be properly reconstructed. In addition to this, in
the case of homogeneous high density of signals across a small frequency band,
we have shown that when the signals are spread in a frequency interval smaller
than $\sim 0.1$ Hz, the overall efficiency increases thanks to the presence of
signal clusters, as the probability to select a peak at a given time is
enhanced. When the signals are spread on a much larger frequency range of
$\sim 0.8$ Hz, we find an efficiency loss of up to $\sim 15\%$, when the mean
signal density approaches a value of two, due to the impact the signals have
in the AR estimation of the average power spectrum. On the other hand, the
“worst” cases considered in ref:papaetc  correspond to an mean signal density
significantly smaller than two, so that the corresponding sensitivity loss is
reduced to a few percent at the most.
We conclude that the frequency-Hough procedure is robust with respect to the
presence of signal ensembles, with negligible losses even in extreme cases.

The results reported here are relevant also in view of searches using the next
generation detectors. The expected improved sensitivity of Einstein Telescope
ref:futureET  with respect to current detectors, especially in the band 3-20
Hz, could make the issue of signal ensembles important. In particular, it
could play a role not only for the case of emission from boson clouds around
Kerr BHs, but especially for the early inspiral of NS binary systems, that
would produce long duration signals and could be searched adapting techniques
derived from standard CW searches. Something similar is expected for the LISA
detector ref:LISA . Our result shows that the basic elements of FrequencyHough
algorithm could be used to resolve most of these sources, exploiting their
different sky position and the algorithm robustness with respect to the
presence of signal ensembles.

## References

  * (1)

##  IX References

* (2) P. D. Lasky, PASA 32, e034 (2015)
* (3) Piccinni, O. J. (2022). Galaxies, 10(3), 72.
* (4) The White paper of the LSC and Virgo collaboration (and Executive Summaries of searches) https://dcc.ligo.org/LIGO-T1700214/public.
* (5) LIGO Scientific collaboration, Virgo collaboration Phys. Rev. X 9, 031040 (2019).
* (6) LIGO Scientific collaboration, Virgo collaboration Phys. Rev. X 11, 021053 (2021).
* (7) LIGO Scientific collaboration, Virgo collaboration arXiv:2108.01045, submitted to PRD.
* (8) LIGO Scientific collaboration, Virgo collaboration, KAGRA Collaboration arXiv:2111.03606, submitted to PRX.
* (9) LIGO Scientific collaboration, Virgo collaboration, Astrophysical Journal 875, 160 (2019).
* (10) Abbott, R., et al.Astrophysical Journal 932, 133 (2022).
* (11) Brito, R., et al.Phys. Rev. D 96, 064050 (2017)
* (12) LIGO Scientific collaboration, Virgo collaboration, Phys. Rev. D 100, 024004 (2019).
* (13) LIGO Scientific collaboration, Virgo collaboration, KAGRA Collaboration Phys. Rev. D 104, 082004 (2021).
* (14) LIGO Scientific collaboration, Virgo collaboration, arXiv:2201.00697 (2022), submitted to PRD.
* (15) Tenorio, R., Keitel, D., Sintes, A. M. (2021). Universe, 7(12), 474.
* (16) P.Astone, A.Colla, S.DAntonio, S.Frasca, C.Palomba,Phys. Rev. D 90, 042002 (2014)
* (17) Iuri La Rosa, P. Astone, S. D’ Antonio, et al., Universe 7,218 (2021)
* (18) S. Zhu, M. Baryakhtar, M.A. Papa, et al. ,Phys. Rev. D 102, 063020 (2020).
* (19) O. J. Piccinni, P. Astone, S. D’ Antonio et al., Classical Quantum Gravity 36:1, 015008 (2019).
* (20) P. Astone, S. Frasca and C. Palomba, Class Quantum Grav 22, S1197 (2005)
* (21) C. Palomba, S. D’Antonio, P. Astone, et al., Phys. Rev. Lett. 123, 171101 (2019).
* (22) S. D’Antonio, C. Palomba, P. Astone, et al. ,Phys. Rev. D 98, 103017 (2018)
* (23) LIGO Scientific collaboration, Virgo collaboration, KAGRA Collaboratio ,Phys. Rev. D 105, 102001 (2022)
* (24) A. Arvanitaki et al. ,Phys. Rev. D 81, 123530 (2010)
* (25) A. Arvanitaki et al. ,Phys. Rev. D 83, 044026 (2011)
* (26) Reitze, D. et al., (2019). arXiv preprint arXiv:1907.04833.
* (27) M. Maggiore, C. Van Den Broek, N. Bartolo et al.arXiv:1912.02622v4  2020\.
* (28) A. La Rana, _Eur. Phys. J. H_ 47, 3 (2022)
* (29) S. Babak, J. Gair, A. Sesana et al, Phys. Rev. D 95, 103012 (2017), and references therein (this is the fifth science case).
* (30) Abbott, R., et al.(2021). SoftwareX, 13, 100658.
* (31) D’Onofrio, L., et al., Phys. Rev. D 105(6), 063012 (2022)
* (32) Bennett, M. F., Melatos, et al., Astrophysical Journal 766(2), 99 (2013)
* (33) Amaro-Seoane, Pau, et al., Classical Quantum Gravity29.12 (2012): 124016.

## Acknowledgement

We thank the National Science Foundation for the international REU program,
grants No. NSF PHY-1950830 and No. NSF PHY-1460803.
This research has made use of data or software obtained from the Gravitational
Wave Open Science Center (gw-openscience.org), a service of LIGO Laboratory,
the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA. LIGO
Laboratory and Advanced LIGO are funded by the United States National Science
Foundation (NSF) as well as the Science and Technology Facilities Council
(STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of
Niedersachsen/Germany for support of the construction of Advanced LIGO and
construction and operation of the GEO600 detector. Additional support for
Advanced LIGO was provided by the Australian Research Council. Virgo is
funded, through the European Gravitational Observatory (EGO), by the French
Centre National de Recherche Scientifique (CNRS), the Italian Istituto
Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions
by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan,
Monaco, Poland, Portugal, Spain. The construction and operation of KAGRA are
funded by Ministry of Education, Culture, Sports, Science and Technology
(MEXT), and Japan Society for the Promotion of Science (JSPS), National
Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea,
Academia Sinica (AS) and the Ministry of Science and Technology (MoST) in
Taiwan.
We also thank the Amaldi Research Center, for the clusters hosted in Rome
INFN, where we have stored the LIGO/Virgo data used in this research and run
part of the present analysis. We thank the INFN-CNAF computing staff for the
resources we have used in this analysis and for their constant support.
This project has received funding from the European Union’s Horizon2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement N. 754496.